本文探討了Python指令碼與動態模態分解(DMD)的結合應用。我們將利用Python對從OpenFOAM模擬中提取的二維切片資料進行DMD計算。這種方法能夠有效地提取隱藏的流動模式,深化對流體動力學現象的理解。
使用開源CFD軟體OpenFOAM,有兩種方法可以對CFD資料進行DMD計算。第一種方法是直接將OpenFOAM的場資料讀入Python;第二種方法則是從OpenFOAM中提取二維切片,然後對這些資料進行DMD計算。
本文將重點介紹第二種方法,即利用Python的強大庫直接分析從OpenFOAM提取的二維切片資料,執行DMD並視覺化提取的模態。
OpenFOAM案例模態分解準備指南
本研究的起點是雷諾數為100的方形圓柱周圍完全發展的、統計穩定的流動。在此基礎上,我們將模擬時間延長至100個渦脫落週期。在每個脫落週期內,我們從資料中提取16次二維切片。二維切片的提取是透過OpenFOAM中的
surfaces函式物件實現的,具體配置如下:
surfaces
{
type surfaces;
libs ("libsampling.so");
writeControl timeStep;
writeInterval 142;
surfaceFormat vtk;
fields (p U);
interpolationScheme cellPoint;
surfaces
{
zNormal
{
type cuttingPlane;
point (0 0 0.05);
normal (0 0 1);
interpolate true;
}
};
};
// ************************************************************************* //模擬完成後,在案例的
postProcessing目錄中會生成一個名為
surfaces的子目錄,其中包含所有提取的表面資料。目錄結構如下:
surfaces/
├── 4771.2000000577236
│ └── zNormal.vtp
├── 4772.6200000577546
│ └── zNormal.vtp
├── 4774.0400000577856
│ └── zNormal.vtp
├── 4775.4600000578166
│ └── zNormal.vtp
.
.
.在進行後續分析之前,請確保案例模擬已完成且表面資料已成功提取。
表面資料提取
為了從OpenFOAM生成的VTK檔案中提取資料,我們將使用PyVista庫。PyVista是視覺化工具包(VTK)的Python介面,透過NumPy包裝VTK庫,提供了直接訪問陣列的方法和類。它為VTK的強大視覺化後端提供了一個文件完善的Pythonic介面,便於快速原型設計、分析和空間參考資料集的視覺化整合。
PyVista在科學計算視覺化中具有重要價值,尤其適用於演示和研究論文的圖形生成。同時它也作為其他依賴3D網格渲染的Python模組的支援庫。
匯入必要的模組,包括PyVista:
importmatplotlib.colors
importmatplotlib.pyplotasplt
importnumpyasnp
importpandasaspd
importfluidfoamasfl
importscipyassp
importos
importmatplotlib.animationasanimation
importpyvistaaspv
importimageio
importio
%matplotlibinline
plt.rcParams.update({'font.size' : 18, 'font.family' : 'Times New Roman', "text.usetex": True})接下來設定路徑變數和常量:
### 常量
d=0.1
Ub=0.015
### 路徑
Path='E:/deephub/Sq_Cyl_Surfaces/surfaces/'
save_path='E:/deephub/SquareCylinderData/'
Files=os.listdir(Path)現在可以嘗試讀取第一個快照表面:
Data=pv.read(Path+Files[0] +'/zNormal.vtp')
grid=Data.points
x=grid[:,0]
y=grid[:,1]
z=grid[:,2]
rows, columns=np.shape(grid)
print('rows = ', rows, 'columns = ', columns)
print(Data.array_names)輸出:
['TimeValue', 'p', 'U']從輸出可以看出,我們的二維切片包含了時間值、壓力場和速度場。利用PyVista,可以為每個快照提取渦量場,並將結果資料組織成一個大型矩陣,以便進行後續的POD計算。具體實現如下:
Data=pv.read(Path+Files[0] +'/zNormal.vtp')
grid=Data.points
x=grid[:,0]
y=grid[:,1]
z=grid[:,2]
rows, columns=np.shape(grid)
print('rows = ', rows, 'columns = ', columns)
### 對U場進行處理
Snaps=len(Files) # 快照數量
data_Vort=np.zeros((rows,Snaps-1))
foriinnp.arange(0,Snaps-1):
data=pv.read(Path+Files[i] +'/zNormal.vtp')
gradData=data.compute_derivative('U', vorticity=True)
grad_pyvis=gradData.point_data['vorticity']
data_Vort[:,i:i+1] =np.reshape(grad_pyvis[:,2], (rows,1), order='F')
np.save(save_path+'VortZ.npy', data_Vort)讓我們檢查一下生成的
data_Vort陣列的維度:
data_Vort.shape
### 輸出
### (96624, 1600)此外,我們可以視覺化渦量場的一個快照:
這個視覺化結果展示了方形圓柱周圍的渦量分佈,為我們提供了流場結構的直觀認識。
正交分解(POD)
為了確定動態模態分解(DMD)的最佳近似秩，我們可以對渦量場資料進行正交分解(POD)分析。POD是一種強大的降維技術，能夠捕捉流場中的主要能量結構。
以下是POD分析的Python實現：
### POD分析
# 構建資料矩陣
X=data_Vort
# 計算並去除平均場
X_mean=np.mean(X, axis=1)
Y=X-X_mean[:,np.newaxis]
# 計算協方差矩陣
C=np.dot(Y.T, Y)/(Y.shape[1]-1)
# 對協方差矩陣進行奇異值分解
U, S, V=np.linalg.svd(C)
# 計算POD模態
Phi_POD=np.dot(Y, U)
# 計算時間係數
a=np.dot(Phi_POD.T, Y)接下來可以分析POD特徵值以評估各模態的能量貢獻：
Energy=np.zeros((len(S),1))
foriinnp.arange(0,len(S)):
Energy[i] =S[i]/np.sum(S)
X_Axis=np.arange(Energy.shape[0])
heights=Energy[:,0]
fig, axes=plt.subplots(1, 2, figsize= (12,4))
ax=axes[0]
ax.plot(Energy, marker='o', markerfacecolor='none', markeredgecolor='k', ls='-', color='k')
ax.set_xlim(0, 20)
ax.set_xlabel('Modes')
ax.set_ylabel('Energy Content')
ax=axes[1]
cumulative=np.cumsum(S)/np.sum(S)
ax.plot(cumulative, marker='o', markerfacecolor='none', markeredgecolor='k', ls='-', color='k')
ax.set_xlabel('Modes')
ax.set_ylabel('Cumulative Energy')
ax.set_xlim(0, 20)
plt.show()
分析結果顯示，前21個POD模態捕捉了約99.9%的總能量。這一發現為我們後面選擇DMD的近似秩提供了重要依據，表明使用21階近似進行DMD分析是合理的。
以下是前幾個POD模態的視覺化結果，用於參考：
這些模態圖展示了流場中的主要結構，為我們理解流動特性提供了直觀的洞察。
動態模態分解(DMD)
動態模態分解是一種強大的技術，能夠提取流場中的動態特徵。以下是DMD演算法的Python實現：
defDMD(X1, X2, r, dt):
# 對X1進行奇異值分解
U, s, Vh=np.linalg.svd(X1, full_matrices=False)
# 截斷SVD矩陣
Ur=U[:, :r]
Sr=np.diag(s[:r])
Vr=Vh.conj().T[:, :r]
# 構建Atilde矩陣並計算其特徵值和特徵向量
Atilde=Ur.conj().T@X2@Vr@np.linalg.inv(Sr)
Lambda, W=np.linalg.eig(Atilde)
# 計算DMD模態
Phi=X2@Vr@np.linalg.inv(Sr) @W
# 計算連續時間特徵值
omega=np.log(Lambda)/dt
# 計算DMD模態振幅
alpha1=np.linalg.lstsq(Phi, X1[:, 0], rcond=None)[0]
b=np.linalg.lstsq(Phi, X2[:, 0], rcond=None)[0]
# DMD重構
time_dynamics=None
foriinrange(X1.shape[1]):
v=np.array(alpha1)[:,0]*np.exp( np.array(omega)*(i+1)*dt)
iftime_dynamicsisNone:
time_dynamics=v
else:
time_dynamics=np.vstack((time_dynamics, v))
X_dmd=np.dot(np.array(Phi), time_dynamics.T)
returnPhi, omega, Lambda, alpha1, b, X_dmd為了應用這個DMD函式，我們首先需要準備時間偏移的資料矩陣：
# 獲取資料矩陣的兩個時間步長偏移檢視
X1=np.matrix(X[:, 0:-1])
X2=np.matrix(X[:, 1:])然後，我們定義近似秩和時間步長：
r=21 # 根據POD分析結果選擇
dt=0.01*142接下來，我們執行DMD計算：
Phi, omega, Lambda, alpha1, b, X_dmd=DMD(X1, X2, r, dt)在進行視覺化之前，我們首先分析DMD特徵值的分佈。這有助於我們理解所識別的DMD模態的動態特性。我們將實部和虛部特徵值繪製在單位圓上：
theta=np.linspace(0, 2*np.pi, 150)
radius=1
a=radius*np.cos(theta)
b=radius*np.sin(theta)
fig, ax=plt.subplots()
ax.scatter(np.real(Lambda), np.imag(Lambda), color='r', marker='o', s=100)
ax.plot(a, b, color='k', ls='--')
ax.set_xlabel(r'$\Lambda_r$')
ax.set_ylabel(r'$\Lambda_i$')
ax.set_aspect('equal')
plt.show()
這個圖顯示所有特徵值都位於單位圓上，表明DMD模態既不增長也不衰減，呈現穩定的特性。
為了視覺化DMD模態，我們首先需要將DMD模態矩陣轉換為陣列：
A=np.squeeze(np.asarray(Phi))然後可以使用Matplotlib繪製DMD模態：
Rect1=plt.Rectangle((-0.5, -0.5), 1, 1, ec='k', color='white', zorder=2)
Mode=11
fig, ax=plt.subplots(figsize=(11, 4))
p=ax.tricontourf(x/0.1, y/0.1, np.real(A[:,Mode]), levels=1001, vmin=-0.005, vmax=0.005, cmap=cmap)
ax.add_patch(Rect1)
ax.xaxis.set_tick_params(direction='in', which='both')
ax.yaxis.set_tick_params(direction='in', which='both')
ax.xaxis.set_ticks_position('both')
ax.yaxis.set_ticks_position('both')
ax.set_xlim(-1, 20)
ax.set_ylim(-5, 5)
ax.set_aspect('equal')
ax.set_xlabel(r'$\bf x/d$')
ax.set_ylabel(r'$\bf y/d$')
ax.text(0, 4, r'$f_i ='+str(np.imag(Lambda[Mode])) +'$', fontsize=25, color='black')
plt.show()
這個圖展示了第11個DMD模態的空間結構。類似地，我們可以繪製前6個DMD模態：
這些DMD模態圖揭示了流場中的關鍵動態結構，為我們深入理解方形圓柱周圍的流動特性提供了重要依據。
透過結合POD和DMD分析，我們不僅捕捉了流場的主要能量結構，還揭示了這些結構隨時間的演化特性。這種綜合分析方法為複雜流動系統的研究提供了強大的工具，能夠幫助我們更深入地理解流體動力學現象。
總結
本文詳細介紹了一種基於OpenFOAM和Python的流場動態分析方法。我們從OpenFOAM模擬資料的提取和處理開始，利用PyVista庫高效地處理二維切片資料。透過正交分解(POD)成功捕捉了流場的主要能量結構，為動態模態分解(DMD)的應用奠定了基礎。DMD分析進一步揭示了流場的動態特徵，使我們能夠深入理解方形圓柱周圍的複雜流動現象。
這種結合OpenFOAM、POD和DMD的綜合分析方法，不僅提高了對複雜流體系統的認識，還為流體動力學研究提供了強大的工具。Python的靈活性和效率在整個分析過程中發揮了關鍵作用，展示了其在科學計算和資料視覺化方面的優勢。
https://avoid.overfit.cn/post/7d6faa4f21244df0ac7ed62f9833acd2
作者：Shubham Goswami