C++(Python)肥皂泡沫普拉托邊界膜曲面模型算法

🎯要點

🎯肥皂泡二維流體模擬 | 🎯泡沫普拉托邊界膜曲面模型算法演化厚度變化 | 🎯螺旋曲面三周期最小結構生成

📜皂膜用例：Python計算物理粒子及拉格朗日和哈密頓動力學 | Python和MATLAB粘性力接觸力動態模型半隱式歐拉算法

在這里插入圖片描述

🍇Python(MATLAB)三維曲面圖

普拉托定律描述了肥皂膜的結構。這些定律是由比利時物理學家普拉托在 19 世紀根據他的實驗觀察制定的。自然界中的許多圖案都是基于遵守這些定律的泡沫。此定律描述了皂膜的形狀和構造如下：

肥皂膜由整個（完整的）光滑表面制成
肥皂膜的一部分的平均曲率在同一塊肥皂膜上的任何點處處處恒定
肥皂膜總是沿著稱為普拉托邊界的邊緣成三塊相交，并且它們以 arccos(? 1 / 2 ) = 120°‘
這些普拉托邊界在一個頂點處以四面體角 arccos(? 1 / 3 ) ≈ 109.47°。

除普拉托定律之外的結構是不穩定的，并且薄膜將很快傾向于重新排列自身以符合這些定律。美國數學家吉恩·艾倫·泰勒使用幾何測量理論在數學上證明了這些定律適用于最小曲面。

肥皂膜是被空氣包圍的薄層液體（通常為水基）。例如，如果兩個肥皂泡接觸，它們會合并并在其間形成一層薄膜。因此，泡沫由通過普拉托邊界連接的薄膜網絡組成。肥皂膜可用作極小曲面的模型系統，極小曲面在數學中被廣泛使用。

從數學角度來看，肥皂膜是最小表面。表面張力是單位面積產生表面所需的能量。薄膜——與任何物體或結構一樣——傾向于以最小勢能狀態存在。為了最小化其能量，自由空間中的液滴自然呈現球形，對于給定的體積，其表面積最小。水坑和薄膜可以在其他力的存在下存在，例如重力和對基質原子的分子間吸引力。后一種現象稱為潤濕：基質原子和薄膜原子之間的結合力會導致總能量降低。在這種情況下，物體的最低能量配置是盡可能多的薄膜原子盡可能靠近基質。這將導致無限薄的薄膜，無限廣泛地分布在基質上。實際上，粘附潤濕效應（導致表面最大化）和表面張力效應（導致表面最小化）會相互平衡：穩定的結構可以是液滴、水坑或薄膜，具體取決于作用于身體的力。

💦Python繪制三維曲面：

Matplotlib 的 mpl_toolkits.mplot3d 工具包中的 axis3d 提供了用于創建三維曲面圖的必要函數。曲面圖是通過使用 ax.plot_surface() 函數創建的。

語法：

ax.plot_surface(X, Y, Z)

其中 X 和 Y 是 x 和 y 點的二維數組，而 Z 是高度的二維數組。

示例：

# Import libraries
from mpl_toolkits import mplot3d
import numpy as np
import matplotlib.pyplot as pltx = np.outer(np.linspace(-3, 3, 32), np.ones(32))
y = x.copy().T # transpose
z = (np.sin(x **2) + np.cos(y **2) )fig = plt.figure(figsize =(14, 9))
ax = plt.axes(projection ='3d')ax.plot_surface(x, y, z)
plt.show()

梯度曲面圖是三維曲面圖與二維輪廓圖的組合。在此圖中，三維曲面的顏色與二維輪廓圖相同。曲面較高的部分與曲面較低的部分顏色不同。

語法：

surf = ax.plot_surface(X, Y, Z, cmap=, linewidth=0, antialiased=False)

屬性 cmap= 設置表面的顏色。還可以通過調用Fig.colorbar來添加顏色條。下面的代碼創建一個梯度曲面圖：

from mpl_toolkits import mplot3d
import numpy as np
import matplotlib.pyplot as pltx = np.outer(np.linspace(-3, 3, 32), np.ones(32))
y = x.copy().T # transpose
z = (np.sin(x **2) + np.cos(y **2) )fig = plt.figure(figsize =(14, 9))
ax = plt.axes(projection ='3d')
my_cmap = plt.get_cmap('hot')surf = ax.plot_surface(x, y, z,cmap = my_cmap,edgecolor ='none')fig.colorbar(surf, ax = ax,shrink = 0.5, aspect = 5)ax.set_title('Surface plot')
plt.show()

使用 Matplotlib 繪制的三維曲面圖可以投影到二維曲面上。下面的代碼創建三維圖并可視化其在二維等高線圖上的投影：

from mpl_toolkits import mplot3d
import numpy as np
import matplotlib.pyplot as pltx = np.outer(np.linspace(-3, 3, 32), np.ones(32))
y = x.copy().T # transpose
z = (np.sin(x **2) + np.cos(y **2) )fig = plt.figure(figsize =(14, 9))
ax = plt.axes(projection ='3d')my_cmap = plt.get_cmap('hot')surf = ax.plot_surface(x, y, z, rstride = 8,cstride = 8,alpha = 0.8,cmap = my_cmap)
cset = ax.contourf(x, y, z,zdir ='z',offset = np.min(z),cmap = my_cmap)
cset = ax.contourf(x, y, z,zdir ='x',offset =-5,cmap = my_cmap)
cset = ax.contourf(x, y, z, zdir ='y',offset = 5,cmap = my_cmap)
fig.colorbar(surf, ax = ax, shrink = 0.5,aspect = 5)ax.set_xlabel('X-axis')
ax.set_xlim(-5, 5)
ax.set_ylabel('Y-axis')
ax.set_ylim(-5, 5)
ax.set_zlabel('Z-axis')
ax.set_zlim(np.min(z), np.max(z))
ax.set_title('3D surface having 2D contour plot projections')plt.show()

💦MATLAB最小曲面

構造最小曲面

n=35;
[X,Y,Z]=meshgrid(linspace(-pi,pi,n));
[F,V] = isosurface(X,Y,Z,S,0);

曲面可視化

cFigure;
hold on;
title('Schwarz P-surface','FontSize',fontSize);
gpatch(F,V,'kw','k',1);
axisGeom;
camlight headlight;
drawnow;

可視化曲面所有變體：

cFigure;subplot(2,3,1);
title('Schwarz P-surface','FontSize',fontSize);
hold on;
gpatch(F,V,pColors(1,:),'none',1);
axisGeom;
camlight headlight; lighting gouraud;
view(-50,30);[X,Y,Z]=meshgrid(linspace(-pi,pi,n));
[F,V] = isosurface(X,Y,Z,S,0.1);subplot(2,3,2);
title('d','FontSize',fontSize);
hold on;
gpatch(F,V,pColors(2,:),'none',1);
axisGeom;
camlight headlight; lighting gouraud;
view(-50,30);[X,Y,Z]=meshgrid(linspace(-2*pi,2*pi,n));
[F,V] = isosurface(X,Y,Z,S,0.6);subplot(2,3,3);
title('Gyroid','FontSize',fontSize);
hold on;
gpatch(F,V,pColors(3,:),'none',1);
axisGeom;
camlight headlight; lighting gouraud;
view(-50,30);[X,Y,Z]=meshgrid(linspace(-pi,pi,n));
[F,V] = isosurface(X,Y,Z,S,0);subplot(2,3,4);
title('Neovius','FontSize',fontSize);
hold on;
gpatch(F,V,pColors(4,:),'none',1);
axisGeom;
camlight headlight; lighting gouraud;
view(-50,30);[X,Y,Z]=meshgrid(linspace(-pi,pi,n));
[F,V] = isosurface(X,Y,Z,S,0);subplot(2,3,5);
title('w','FontSize',fontSize);
hold on;
gpatch(F,V,pColors(5,:),'none',1);
axisGeom;
camlight headlight; lighting gouraud;
view(-50,30);[X,Y,Z]=meshgrid(linspace(-pi,pi,n));
[F,V] = isosurface(X,Y,Z,S,0.5);subplot(2,3,6);
title('pw','FontSize',fontSize);
hold on;
gpatch(F,V,pColors(6,:),'none',1);
axisGeom;
camlight headlight; lighting gouraud;
view(-50,30);
drawnow;