验证OpenFOAM中interFoam求解器Word格式.docx
《验证OpenFOAM中interFoam求解器Word格式.docx》由会员分享，可在线阅读，更多相关《验证OpenFOAM中interFoam求解器Word格式.docx（16页珍藏版）》请在冰豆网上搜索。
在德国，许多研究院校已经开始用OpenFOAM来工作，他们几乎每半年就组织一
次OpenFOAM的用户交流研讨会；
挪威的Risavika燃气能源公司开展了数值模拟传热多相
流研究二氧化碳的捕捉和封存，该研究机构已经发表了他们的研究成果并给出了含传热过程
的可压缩VOF两相流求解器（compressibleInterFoam）的算例；
另外，我们也可以从CFD-Online论坛上看到OpenFOAM板块的帖子数仅次于ANSYS（FluentandCFX），其用户交流活跃程
度非同一般。
当然，基于OpenFOAM开展工作也面临着一定程度的挑战和困难，相比较商
业CFD软件，OpenFOAM缺乏详尽的用户使用手册，相关学习资料也很有限，它对研究工
作者提出了更高的要求，需要研究者更加积极主动地去探索和持续不断地深入学习。
CFD-Online的OpenFOAM板块是很好学习交流平台，OpenFOAMwiki是很好的资料库，对审视OpenFOAM的所有相关内容很有帮助，Jasak的克罗地亚网站提供了很多演讲报告和
OpenFOAM相关的博士论文，是深入学习的好资料，官方发布的UserGuide和ProgrammerGuide，尽管写得远远不够详尽，初学者会觉得对读程序没有多大帮助，但是只
要对OpenFOAM自定义化的语言一定程度的熟悉之后，再回过头来看，会对OpenFOAM的认识有一个整体上的提升。
2.interFOAM
2.1（twofluidsystem）
interFOAM是OpenFOAM中用来解决two-fluidsystems一类问题的求解器，属于theflowoftwoimmisciblefluids的现象有：
水波、溃坝、水中气泡等，two-fluidsystem按照界面的拓扑结构可分为三类：
分离（segregated）、混合流（mixed）、分散流（dispersed）。
比如一个盛有部分水的封闭容器，当低频小幅度晃动时，界面完好，此时为第一种情况；
当晃动频率振幅增
加到一定时，界面变得不稳定，界面部分破碎，致使液体包进了气泡，此属第二类情况；
当
容器剧烈晃动时，大量空气以悬浮的气泡形式存在液体中，为第三类情况。
interFOAM中没有引入湍流模型，因此它主要针对第一二类问题，所有的方法能够模拟界面发生较大变形，
破裂或融合等现象，但也不排除第三类情况。
这类问题的关键技术在于要模拟出两种流体的
界面interface，由于在整个体系的运动过程中，界面往往是不断变化的，整个流体域的运动
相互耦合，从而需要在求解过程中不断生成界面，要解决好这个问题，需要考虑以下几点：
1.在离散网格中表示出界面
2.如何处理分布了两种流体的计算单元
3.界面的运动
4.界面条件与运动方程的耦合
总的来说，interFOAM中使用了VOF（volumeoffluid）来处理这类问题。
自由面或流
体界面的模拟方法总的可归为两类：
面法（假想界面上有粒子，高度函数，LevelSet法，贴面法等）和体积法（假想流体中有粒子，体积分数法等）。
interFOAM使用了体积法（volumemethod）中的体积分数法（volumefraction）来表示计算域中两种流体的分布，但是仅仅知道
计算单元（cell）中的两种流体体积比还是不够的（图1），还要能进一步表示出每个计算单元
中的流体分界面，这需要良好的interfacecapturingmethod，以期望模拟出真实的界面。
目前
出现的技术有：
线技术（linetechniques），施主受主法（donor-acceptermethod），高阶差分法（higherorderdifferencingschemes），interFOAM使用了基于施主受主法的VOF法。
图1.离散网格中的体积分数
2.2
对于two-fluidsystem这类问题，首先引入体积分数函数,，定义表示一种流体,,1
的存在，,表示另一种流体的存在，表示两种流体的界面区域，因此,,00,,,1
（）表示了整个流场中两种流体的分别，在建立运动方程时将两种流体当作一种0,,,1
流体来处理，该流体的物理属性表示为：
,,,，（1,,）,
（1）12
,,,，（1,,）,
（2）12
OpenFOAM采用FVM法离散计算区域，,函数在界面区域经历非连续变化（astepfunction），这导致成为分段连续函数，为了能够使用连续介质力学方法建立数学模型,，,
并计算界面曲率，需要对,函数进行光滑处理，定义流体的过渡区域厚度为，则有：
1forthepoint（x,t）insidefluid1,
（3）
（,）0forthepoint（x,t）insidefluid2xt,
,01forthepoint（x,t）insidethetransitionalarea,
满足的输运方程：
D,,,（4）,，u,,,,0Dt,t
若用表示任何守恒物理量，它满足控制体形式的守恒方程为：
,dV，F,dS,F,dS,Q,dV，Q,dS（5）,,,,,cDVs,tVVVV,,,V
边界处的扩散项源项边界处的源项边界处的对流项方程（5）的微分形式为：
,（6），,,F,,,F,Q，,,QcDVs,t
取，表示单位体积质量，由于我们研究的系统中无化学反应或相变，因此无源无,,,
扩散，方程简化为：
,,（7），,,,u,0,t
利用方程
（1）（4）可进一步简化为
（8）,,u,0
动量守恒时，取F,0,,,u，表示单位体积动量，认为系统中无动量扩散，因此，D
源项由内部力和外部力组成，流体微元内部力相互抵消，微元表面存在应力，切向表现为粘
性力，法向表现为压力，把流体视为局部热平衡的牛顿流，则应力张量表示为：
,,2,,T，，，，T,,p，,,,uI，,,,u，,,u（9）,,3,,
另外，液体在自由面上还受到表面张力的影响，,是张力系数，表示自由面每增加单
位面积所需做的功，用,f表示界面曲率，表示界面的单位法向量，用表示表面张力沿n,界面法向的一个分量，此力的作用是平衡界面两边的压力差，此力只在界面处存在，在非界
面处其值为零。
f,,,（x）n（10）,
其中
,,（11）n,,,
（12）,（x）,,,n
因此，得到动量方程：
,,,,u，，，,,,u,u,T,,g，f,,t
由方程（8）（9）进一步可化为：
,,,,,u（14），,,（,uu）,,,,,u,,g,,,p,f,,t
在interFOAM中，为了解决阶段函数,的迁移问题，Weller引入了额外的人工压缩项：
,,,，，，，，，（15），,,,u，,,,1,,u,0,,t
其中u,是一个适于压缩界面的速度场，这一项只对界面区域产生影响。
函数在界面,
区域经历非连续变化，这导致密度粘性成为分段连续函数，为了能够使用连续介质力学方法
建立数学模型并计算界面曲率，还需要采取对,函数进行光滑处理等相关技术。
2.3
OpenFOAM采用基于任意非结构化网格的有限体积法，将流体区域分成有限数量的互
不重叠的小控制体，然后将微分形式的控制方程在每个小控制体积分，保证方程局部和整体
上的守恒性，然后用分段之和近似积分，并提供各种差分差值方案（Eulerimplicitorexplicit,Crank-Nicolson,centred,upwind,TVD,NVD等）来估算各种微分运算（面法向梯度，法向梯
度，拉普拉斯项，散度项，时间项，flux计算），来近似表现流体的迁移或扩散等行为，使
用时要就具体问题加以选择。
在求解代数方程组时，也有多种求解控制器可供选择（LinearSolverControl,PISOandSIMPLE等）。
在用interfacecapturing方法处理界面的情况中，容易出现数值扩散，导致界面受污，为了克服此问题，常采用界面压缩技术，在OpenFOAM中求解
方程时，Weller引入了他自己的一种技术叫做interfaceCompression，来解决传统的VOF界面压缩法存在的基本问题。
3interFOAM
利用interFOAM可以求解很大一部分属于two-fluidsystems的问题，所有的方法上一章已经介绍过，各种具体的问题的差异性只在于体积分数，压力，速度的初始条件和边界条
件的不同。
3.1
（1）.初始条件和边界条件
图2初始时刻的水柱
溃坝的初始时刻如图2所示，左右和底部是墙（无滑移），此处速度始终为零；
因为模
拟的是二维效应，所以模型的前后面设置成滑移条件（OpenFOAM本是三维模拟）；
顶部模拟的是开边界，因此在此处设置空气压力为定值，速度法向梯度为零，两种流体可通过此
边界离开计算区域，但只有空气可经此流入计算域，因此速度指向外时，体积分数为零梯度，
速度向内时，体积分数为定值；
初始时刻，速度为零，边界处压力梯度为零，体积分数如图2在水柱区域为1，其他区域为零。
为了便于后面与实验结果比较，这里取得几何模型与实
验中的一致。
（2）．在OpenFOAM中求解
要用OpenFOAM模拟溃坝，需要模型的几何信息，输运性质（密度粘性等），环境信
息（重力场），这些信息分别保存在“constant”文件夹中；
需要整个计算域初始时刻的体积
分数，速度和压力信息，它们分别保存在“0”文件夹中；
需要控制方程各项的离散方案，离
散后的代数方程的迭代方法和容许误差等信息，以及求解时需要的时间设置和读写控制等信
息，它们分别保存在“system”文件夹中。
建立了上面所述的输入信息之后，在Linux下的终端依次执行：
“blockMesh”（网格生成）“setFields”（设置流场初始分布）“interFOAM”（求解）“paraFOAM”（后处理），最后可在paraView中得到可视化的结果，计算出的各个时
刻的压力速度体积分数这些数据保存在以时间命名的文件夹里（正如初始时刻数据放在“0”文件夹中一样）。
针对Weller的界面压缩方法，需要对数值方案和压缩系数进行多种尝试
和比较，从而找出最优方案。
（3）．结果分析
溃坝实验是用于验证自由面流或两相流数学模型合理性的经典实验。
这里将所得结果
（图4）与经典实验结果（图3）做定性的比较。
图3到图6都是反映溃坝的时间演化过
程的图片。
在此过程中，重力加速度作用于左边的水柱向下向右运动以寻求可能的最低势能。
初始阶段由惯性力支配，随后粘性效应迅速增大。
下面捕捉几个主要现象加以比较和分析：
1）t=0.2s：
底部约75%的部分被水覆盖，与实验基本一致；
2）t=0.4s:
：
运动在最前端的水流爬上对面的墙上，底部的水平自由面如图3中一样略微倾斜
3）t=0.6s：
底部的水面已经平行于底部，爬上对墙的水流因为重力开始下滑
4）t=0.8s：
下滑后的水流形成一股反向运动的波并包进了部分空气，但是与图4相比，可以观察到空气泡体积偏大，另外返回的波原本在水面上溅起的大量水珠，在数值模拟中未
能捕捉到，以至原本分散到空气中的这部分还保留在波头中，使波头看起来比实验照片中要
低。
这说明模型还不能很好的模拟出更小尺度的破碎（breakup）和扩散（dispersed）现象。
5）t=1.0s：
返回波爬升到左边墙上，由于波头在水面上运动的跳跃性，在爬上墙时包进
了大量空气，空气泡中和水面上的扩散（dispersed）现象同样没有捕捉到，而其他情况基本
一致。
图3溃坝的实验结果
图4溃坝的数值模拟结果
3.2
图5.有障碍物的溃坝初始几何模型
更有趣的溃坝现象是，在水柱不远处增加了一个障碍物（如图5），相比上例此例只在几何建模时略有差异，其他信息完全不变。
模拟出的结果（如图7）：
1）t=0.1s：
实验中（图6）此刻底部流体前段已经触到障碍物，而数值模拟结果显示此
刻还未到达障碍物；
2）t=0.2s：
波头被障碍物挡住，一部分绕过障碍物从上方继续前进；
3）t=0.3s：
波头继续朝对面墙运动，并已盖住了障碍物的上边缘；
4）t=0.4s：
绕过障碍物的波头撞倒对面墙上，将下方的空气困在里面，情况跟实验吻合；
5）t=0.5s：
此时障碍物上方的部分水流由于重力作用开始下坠，撞到墙上的水流也开始
下滑，但下方困住的空气的阻碍作用，可以看到墙上下坠的水流和流过障碍物上边缘的水流，
在它们下方将会形成一个气泡（如图7中黄色小框标记出的）；
6）t=0.6s：
被困住的大气泡中的空气已经冲破上面水层，还可以看到上一时刻中所标示
的已经形成完全被水所包裹的气泡，另外绕到障碍物后方的水流在撞到底部之后，又反溅起
来（图中已标示）。
这些现象跟实验吻合。
同样，整个过程中脱离水面扩散到空气中的大量
水珠未能被模拟出来。
图6.溃坝撞击障碍物的实验结果（Koshizuka（1995））
图7.溃坝撞击障碍物的数值结果
3.3
上面溃坝的数值模拟结果已经在一定程度上显示出了interFOAM模拟这类two-fluidsystems问题的有效范围，interFOAM对自由面变形破裂和融合的模拟效果还是不错的，但
对根据细微的水和空气微粒相互扩散现象的模拟目前难以做到，需要更加细致的网格或模
型。
当然在OpenFOAM中，由于所用方法的不同，解决这类问题还可以选用其他求解器，
比如rasInterFOAM，该求解器中引入了各种湍流模型和DPE模型（DispersedPhaseElement），应该可以模拟现象更加复杂的问题。
4.interFOAM
4.1
（1）数值造波理论X,、角频率为的活塞式造波机，其推板0做简谐运动的速度为：
由线性造波理论，对于平衡位置在原点、冲程为
X,0U（t）,cos,t（16）2
'
在水深为,的波浪水槽中距造波板x处的波面为d
22,,X4sin,d4sinhkd,,x'
0nx（17）,cos（kx,t）esin,t,,，,,,22kdsinh（2kd）2,dsin2,d，，nn,,
式中满足方程k
2kgtanh（kd）,,,0（18）
而,为下面方程的第n个根n
2，，,gtan,d，,,0（19）nn
式（17）中Σ项为造波板产生的衰减立波，在离开造波板一定距离后很快就衰减掉，而第
一项则是造波板所产生的波数为,，频率为的行进波。
k
令式（17）中的x=0，则得到造波前的波面为
WLT’,U（t），U（t,）,（20）,,4其中，
24sinhkdW,2kd，sinh（2kd）
24sin,dnL,,2,d，sin2,dnn
U,为普通造波机产生所需要的余弦波面的推板速度，则有00
若令,,0（21）U（t）,0W
（2）在OpenFOAM中如何实现？
静网格法，通过设置入口处的速度、压力、gamma值，达到造波的目的，通过在waveMakerVelocityFvPatchVectorField.C文件中加入#include“waveMakerVelocity.H”，来调用
waveMakerVelocity.H文件。
waveMakerVelocity.H文件内定义了多种类型的波
相当于定义了新的速度边界条件waveMakerVelocity，在case的定义中直接调用自定义
的边界条件
有OpenFOAM自带的边界条件进行改写，这里根据
OpenFOAM-1.6/src/finiteVolume/fields/fvPatchFields/inletOutlet进行改写。
（3）线性波的数值模拟结果
建立一个长45米，高1.2米，水深为0.6米的二维水池模型，用速度入口式生产一个波高为
0.02米，周期为1s的线性波。
水池中生成的波面如图8、图9。
生成的波面图
图8线性波的波面图
图9距离造波边界处1个、2个、4个、6个波长处的波面时历曲线
4.2如何在OpenFOAM求解器中加海绵层，实现消波的目的为了更有效地解决数值水槽末端的波浪反射问题，在VOF模型中采用了海绵层阻尼消波的方法。
其原理是将某种形式的衰减系数加入到控制方程之中，使得在海绵层内部的流体满足
修正后的时均运动方程：
222,u,u,u,u1,p,u,u,u，u，v，w,g,，v（，，）,,（x）,ux222,t,x,y,z,,x,x,y,z
222,v,v,v,v1,p,v,v,v（22），u，v，w,g,，v（，，）,,（x）,vy222,t,x,y,z,,y,x,y,z
222,w,w,w,w1,p,w,w,w，u，v，w,g,，v（，，）,,（x）,wz222,t,x,y,z,,z,x,y,z模型中的海绵层采用线性衰减形式
x,x当0x,x,（x）,,xx时，，为海绵层的起点，为水池的长度；
001x,x10
当x,x,（x）,0时，。
0
OpenFOAM对方程组（22）写成如下代码进行计算：
fvVectorMatrixUEqn
•（
•fvm:
:
ddt（rho,U）
•+fvm:
div（rhoPhi,U）
•-fvm:
laplacian（muf,U）
•-（fvc:
grad（U）&
fvc:
grad（muf））
•+U*rho*thetaField//spongelayer
可见，在OpenFOAM中实现消波非常的简单，只要修改求解方程的外层代码即可。
增加了海绵层之和，在不同时刻的波面图，如图10。
图10T=30s/40s/50s/60s时的水池的波面图
参考文献：
1.OpenFOAM:
UserGuide,ProgrammerGuide
2.
3.HassanHemida.OpenFOAMtutorial:
FreesurfacetutorialusinginterFOAMand
rasInterFOAM,2008
4.OnnoUbbink.Numericalpredictionoftwofluidsystemswithsharpinterfaces.DoctoralthesisforUniversityofLondon,ImperialCollege,1997
5.HenrikRusche.ComputationalFluidDynamicsofDispersedTwo-PhaseFlowsatHighPhaseFractions.DoctoralthesisforUniversityofLondon,ImperialCollege,2002
升级会员 