基于OpenFOAM和Python的流场动态模态分解:从数据提取到POD-DMD分析

本文涉及的产品
实时计算 Flink 版,5000CU*H 3个月
智能开放搜索 OpenSearch行业算法版,1GB 20LCU 1个月
检索分析服务 Elasticsearch 版,2核4GB开发者规格 1个月
简介: 本文介绍了如何利用Python脚本结合动态模态分解(DMD)技术,分析从OpenFOAM模拟中提取的二维切片数据,以深入理解流体动力学现象。通过PyVista库处理VTK格式的模拟数据,进行POD和DMD分析,揭示流场中的主要能量结构及动态特征。此方法为研究复杂流动系统提供了有力工具。

本文探讨了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

目录
相关文章
|
6天前
|
机器学习/深度学习 数据采集 TensorFlow
使用Python实现智能食品消费模式分析的深度学习模型
使用Python实现智能食品消费模式分析的深度学习模型
92 70
|
26天前
|
数据采集 缓存 定位技术
网络延迟对Python爬虫速度的影响分析
网络延迟对Python爬虫速度的影响分析
|
27天前
|
算法 Python
在Python编程中,分治法、贪心算法和动态规划是三种重要的算法。分治法通过将大问题分解为小问题,递归解决后合并结果
在Python编程中,分治法、贪心算法和动态规划是三种重要的算法。分治法通过将大问题分解为小问题,递归解决后合并结果;贪心算法在每一步选择局部最优解,追求全局最优;动态规划通过保存子问题的解,避免重复计算,确保全局最优。这三种算法各具特色,适用于不同类型的问题,合理选择能显著提升编程效率。
39 2
|
28天前
|
数据采集 存储 JSON
Python爬虫开发中的分析与方案制定
Python爬虫开发中的分析与方案制定
|
28天前
|
图形学 Python
SciPy 空间数据2
凸包(Convex Hull)是计算几何中的概念,指包含给定点集的所有凸集的交集。可以通过 `ConvexHull()` 方法创建凸包。示例代码展示了如何使用 `scipy` 库和 `matplotlib` 绘制给定点集的凸包。
27 1
|
29天前
|
JSON 数据格式 索引
Python中序列化/反序列化JSON格式的数据
【11月更文挑战第4天】本文介绍了 Python 中使用 `json` 模块进行序列化和反序列化的操作。序列化是指将 Python 对象(如字典、列表)转换为 JSON 字符串,主要使用 `json.dumps` 方法。示例包括基本的字典和列表序列化,以及自定义类的序列化。反序列化则是将 JSON 字符串转换回 Python 对象,使用 `json.loads` 方法。文中还提供了具体的代码示例,展示了如何处理不同类型的 Python 对象。
|
29天前
|
数据采集 Web App开发 iOS开发
如何使用 Python 语言的正则表达式进行网页数据的爬取?
使用 Python 进行网页数据爬取的步骤包括:1. 安装必要库(requests、re、bs4);2. 发送 HTTP 请求获取网页内容;3. 使用正则表达式提取数据;4. 数据清洗和处理;5. 循环遍历多个页面。通过这些步骤,可以高效地从网页中提取所需信息。
|
1月前
|
数据可视化 开发者 Python
Python GUI开发:Tkinter与PyQt的实战应用与对比分析
【10月更文挑战第26天】本文介绍了Python中两种常用的GUI工具包——Tkinter和PyQt。Tkinter内置于Python标准库,适合初学者快速上手,提供基本的GUI组件和方法。PyQt基于Qt库,功能强大且灵活,适用于创建复杂的GUI应用程序。通过实战示例和对比分析,帮助开发者选择合适的工具包以满足项目需求。
92 7
|
1月前
|
存储 数据处理 Python
Python科学计算:NumPy与SciPy的高效数据处理与分析
【10月更文挑战第27天】在科学计算和数据分析领域,Python凭借简洁的语法和强大的库支持广受欢迎。NumPy和SciPy作为Python科学计算的两大基石,提供了高效的数据处理和分析工具。NumPy的核心功能是N维数组对象(ndarray),支持高效的大型数据集操作;SciPy则在此基础上提供了线性代数、信号处理、优化和统计分析等多种科学计算工具。结合使用NumPy和SciPy,可以显著提升数据处理和分析的效率,使Python成为科学计算和数据分析的首选语言。
36 3
|
1月前
|
存储 机器学习/深度学习 算法
Python科学计算:NumPy与SciPy的高效数据处理与分析
【10月更文挑战第26天】NumPy和SciPy是Python科学计算领域的两大核心库。NumPy提供高效的多维数组对象和丰富的数学函数,而SciPy则在此基础上提供了更多高级的科学计算功能,如数值积分、优化和统计等。两者结合使Python在科学计算中具有极高的效率和广泛的应用。
52 2