验证OpenFOAM中interFoam求解器

验证OpenFOAM中interFOAM求解器
1.OpenFOAM简介
OpenFOAM是Open Field Operation and Manipulation的简称，顾名思义，它实现的是(张量)场的运算和操作，实质上是一个应用于计算连续介质力学的C++类库。

这一开源软件起源于八十年代后期英国帝国理工大学的Gosman教授的团队。

当时的CFD代码普遍采用FORTRAN语言编写，为了寻求一种更为强大灵活通用的模拟平台，在以后的若干年里，他们利用上了C++语言的高级特性，采用更为有效的方式重新编写了很多代码。

其中最初的一些C++类出自Charlie Hill的博士论文，到1993年由Henry Weller和Hrvoje Jasak以及后来的帝国理工博士生们持续不断的开发。

这一软件起初叫做FOAM，直到2004年由Henry创立的OpenCFD公司以开放源代码的形式公之于众，更名为OpenFOAM，目前最新版本是OpenFOAM-1.7，另外还有一个Jasak引导的克罗地亚版本，目前的版本是OpenFOAM-1.5-dev，两个版本在底层类库和上层应用的涵盖范围上有许多差异。

OpenFOAM 的优越性表现在以下几个方面：
1.是最早利用C++语言编写而成的科学软件包之一（其它的主流CFD软件公司已经发
布或正在开发新一代C++代码）；
2.利用C++的运算符重载功能使得顶层代码在对偏微分方程的描述上相对简单了许
多，且可读性强，这使得OpenFOAM看上去就像一种非常自然的适用于模拟物理问题的编程语言；
3.是最早采用多面体单元网格的通用CFD软件包，而这个功能得以实现是源于对模拟
对象采用分层描述的自然结果；
4.是目前发布于开源许可下的最强大的通用CFD软件包。

OpenFOAM为CFD领域研究工作者们提供了一个强大的开源研究平台，基于OpenFOAM展开的学习和研究目前正在欧美国家的研究院校和技术公司如火如荼的进行：Chalmer University of Technology从2006年开始开设了针对OpenFOAM的博士课程，他领导的涡轮机小组开发出用于模拟空泡和回旋流的代码，在瑞典还有另外几所大学也开设了OpenFOAM相关课程，Niklas Nordin和他的研究小组开发出强大的模拟燃烧的求解器(dieselFoam)，瑞典的几家公司开始将它们应用于工业；2008年丹麦的工业界和研究院校发起了提升国家开源CFD软件的倡议，八家公司和两所院校(Aalborg University和Technical University of Denmark)联合开发丹麦的开源CFD软件平台，他们选用OpenFOAM开展了诸多项目；在德国，许多研究院校已经开始用OpenFOAM来工作，他们几乎每半年就组织一次OpenFOAM的用户交流研讨会；挪威的Risavika燃气能源公司开展了数值模拟传热多相流研究二氧化碳的捕捉和封存，该研究机构已经发表了他们的研究成果并给出了含传热过程的可压缩VOF两相流求解器(compressibleInterFoam)的算例；另外，我们也可以从CFD-Online 论坛上看到OpenFOAM板块的帖子数仅次于ANSYS(Fluent and CFX)，其用户交流活跃程
度非同一般。

当然，基于OpenFOAM开展工作也面临着一定程度的挑战和困难，相比较商业CFD软件，OpenFOAM缺乏详尽的用户使用手册，相关学习资料也很有限，它对研究工作者提出了更高的要求，需要研究者更加积极主动地去探索和持续不断地深入学习。

CFD-Online的OpenFOAM板块是很好学习交流平台，OpenFOAMwiki是很好的资料库，对审视OpenFOAM的所有相关内容很有帮助，Jasak的克罗地亚网站提供了很多演讲报告和OpenFOAM相关的博士论文，是深入学习的好资料，官方发布的UserGuide和ProgrammerGuide，尽管写得远远不够详尽，初学者会觉得对读程序没有多大帮助，但是只要对OpenFOAM自定义化的语言一定程度的熟悉之后，再回过头来看，会对OpenFOAM的认识有一个整体上的提升。

2.基本两相流求解器interFOAM简介
2.1简介两相流(two fluid system)问题
interFOAM是OpenFOAM中用来解决two-fluid systems一类问题的求解器，属于the flow of two immiscible fluids的现象有：水波、溃坝、水中气泡等，two-fluid system按照界面的拓扑结构可分为三类：分离(segregated)、混合流(mixed)、分散流(dispersed)。

比如一个盛有部分水的封闭容器，当低频小幅度晃动时，界面完好，此时为第一种情况；当晃动频率振幅增加到一定时，界面变得不稳定，界面部分破碎，致使液体包进了气泡，此属第二类情况；当容器剧烈晃动时，大量空气以悬浮的气泡形式存在液体中，为第三类情况。

interFOAM中没有引入湍流模型，因此它主要针对第一二类问题，所有的方法能够模拟界面发生较大变形，破裂或融合等现象，但也不排除第三类情况。

这类问题的关键技术在于要模拟出两种流体的界面interface，由于在整个体系的运动过程中，界面往往是不断变化的，整个流体域的运动相互耦合，从而需要在求解过程中不断生成界面，要解决好这个问题，需要考虑以下几点：
1.在离散网格中表示出界面
2.如何处理分布了两种流体的计算单元
3.界面的运动
4.界面条件与运动方程的耦合
总的来说，interFOAM中使用了VOF（volume of fluid）来处理这类问题。

自由面或流体界面的模拟方法总的可归为两类：面法（假想界面上有粒子，高度函数，Level Set 法，贴面法等）和体积法（假想流体中有粒子，体积分数法等）。

interFOAM 使用了体积法(volume method)中的体积分数法（volume fraction）来表示计算域中两种流体的分布，但是仅仅知道计算单元(cell)中的两种流体体积比还是不够的（图1），还要能进一步表示出每个计算单元中的流体分界面，这需要良好的interface capturing method，以期望模拟出真实的界面。

目前出现的技术有：线技术(line techniques)，施主受主法(donor-accepter method)，高阶差分法(higher order differencing schemes)，interFOAM使用了基于施主受主法的VOF法。

图1. 离散网格中的体积分数
2.2 数学建模
对于two-fluid system 这类问题，首先引入体积分数函数α，定义1=α表示一种流体的存在，0=α表示另一种流体的存在，10<<α表示两种流体的界面区域，因此α（10≤≤α）表示了整个流场中两种流体的分别，在建立运动方程时将两种流体当作一种流体来处理，该流体的物理属性表示为：
211ρααρρ）（-+= （1）
211μααμμ）（-+= （2）
OpenFOAM 采用FVM 法离散计算区域，α函数在界面区域经历非连续变化(a step function)，这导致μρ，成为分段连续函数，为了能够使用连续介质力学方法建立数学模型并计算界面曲率，需要对α函数进行光滑处理，定义流体的过渡区域厚度为δ，则有：
1 f o r t h e p o i n t (x ,t ) i n s i d e f l u i d 1(,)0 f o r t h e p o i n t (x ,t ) i n s i d e f l u i d 201 f o r t h e p o i n t (x ,t ) i n s i d e t h e t r a n s i t i o n a l a r e a x t δαα
⎧⎪=⎨⎪<<⎩
α满足的输运方程：
（3）
0=∇⋅+∂∂=αααu t
Dt D （4） 若用ϕ表示任何守恒物理量，它满足控制体形式的守恒方程为：
⎰⎰⎰⋅+⋅=⋅-⋅+
∂s V dS Q dV Q dS F dS F dV ϕ （5） 方程（5）的微分形式为： s V D c Q Q F F t
⋅∇+=⋅∇-⋅∇+∂∂ϕ （6） 取ρϕ=，表示单位体积质量，由于我们研究的系统中无化学反应或相变，因此无源无扩散，方程简化为：
0=⋅∇+∂∂u t
ρρ （7） 利用方程（1）（4）可进一步简化为
0=⋅∇u （8）
动量守恒时，取u
ρϕ=，表示单位体积动量，认为系统中无动量扩散，因此0=D F ，源项由内部力和外部力组成，流体微元内部力相互抵消，微元表面存在应力，切向表现为粘性力，法向表现为压力，把流体视为局部热平衡的牛顿流，则应力张量表示为： ()()T u u I u p T ⊗∇+⊗∇+⎪⎭
⎫ ⎝⎛⋅∇+-=μμ32 （9） 另外，液体在自由面上还受到表面张力的影响，σ是张力系数，表示自由面每增加单位面积所需做的功，用κ表示界面曲率，n
表示界面的单位法向量，用σf 表示表面张力沿界面法向的一个分量，此力的作用是平衡界面两边的压力差，此力只在界面处存在，在非界面处其值为零。

n x f )(σκσ= （10）
其中
α
α∇∇=n （11） n x
⋅∇=)(κ （12）
因此，得到动量方程： ()σρρρf g T u u t
u +=-⊗⋅∇+∂∂ 由方程（8）（9）进一步可化为：
σρμρρf p g u u u t
u --∇=-∇⋅∇-⋅∇+∂∂ )( （14） 在interFOAM 中，为了解决阶段函数α的迁移问题，Weller 引入了额外的人工压缩项：
()()()01=-⋅∇+⋅∇+∂∂σααααu u t
（15） 其中σu 是一个适于压缩界面的速度场，这一项只对界面区域产生影响。

α函数在界面区域经历非连续变化，这导致密度粘性成为分段连续函数，为了能够使用连续介质力学方法建立数学模型并计算界面曲率，还需要采取对α函数进行光滑处理等相关技术。

2.3 数值方法
OpenFOAM 采用基于任意非结构化网格的有限体积法，将流体区域分成有限数量的互不重叠的小控制体，然后将微分形式的控制方程在每个小控制体积分，保证方程局部和整体上的守恒性，然后用分段之和近似积分，并提供各种差分差值方案（Euler implicit or explicit, Crank-Nicolson, centred, upwind, TVD, NVD 等）来估算各种微分运算（面法向梯度，法向梯度，拉普拉斯项，散度项，时间项，flux 计算），来近似表现流体的迁移或扩散等行为，使用时要就具体问题加以选择。

在求解代数方程组时，也有多种求解控制器可供选择（Linear Solver Control, PISO and SIMPLE 等）。

在用interface capturing 方法处理界面的情况中，容易出现数值扩散，导致界面受污，为了克服此问题，常采用界面压缩技术，在OpenFOAM 中求解α方程时，Weller 引入了他自己的一种技术叫做interfaceCompression ，来解决传统的VOF 界面压缩法存在的基本问题。

3．溃坝实例验证interFOAM
利用interFOAM 可以求解很大一部分属于two-fluid systems 的问题，所有的方法上一章已经介绍过，各种具体的问题的差异性只在于体积分数，压力，速度的初始条件和边界条件的不同。

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 相比，可以观察到空气泡体积偏大，另外返回的波原本在水面上溅起的大量水珠，在数值模拟中未能捕捉到，以至原本分散到空气中的这部分还保留在波头中，使波头看起来比实验照片中要低。

这说明模型还不能很好的模拟出更小尺度的破碎（break up）和扩散(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-fluid systems 问题的有效范围，interFOAM对自由面变形破裂和融合的模拟效果还是不错的，但对根据细微的水和空气微粒相互扩散现象的模拟目前难以做到，需要更加细致的网格或模型。

当然在OpenFOAM中，由于所用方法的不同，解决这类问题还可以选用其他求解器，比如rasInterFOAM，该求解器中引入了各种湍流模型和DPE模型(Dispersed Phase Element)，应该可以模拟现象更加复杂的问题。

4.造波验证interFOAM 求解器 4.1 使用速度入口式造波 (1) 数值造波理论
由线性造波理论，对于平衡位置在原点、冲程为0X 、角频率为ω的活塞式造波机，其推板做简谐运动的速度为：
t X t U ωω
cos 2
)(0=
（16）
在水深为d 的波浪水槽中距造波板x 处的波面'
η为
⎭
⎬⎫⎩⎨⎧++-+=∑-t e d d d t kx kd kd kd X x
n n n x ωμμμωημsin 2sin 2sin 4)cos()2sinh(2sinh 42
220
'
（17） 式中k 满足方程
0)tanh(2=-ωkd kg （18）
而n μ为下面方程的第n 个根
()0tan 2=+ωμμd g n n （19）
式（17）中Σ项为造波板产生的衰减立波，在离开造波板一定距离后很快就衰减掉，而第一项则是造波板所产生的波数为k ，频率为ω的行进波。

令式（17）中的x=0，则得到造波前的波面为
)4
()(T
t U L
t U W
-+
=
ωω
η’ （20）
其中，
)
2sinh(2sinh 42kd kd kd
W +=
∑
+=d
d d
L n n n μμμ2sin 2sin 42
若令0U 为普通造波机产生所需要的余弦波面0η的推板速度，则有
W
t U ω
η00)(=
(21)
(2) 在OpenFOAM 中如何实现？
静网格法，通过设置入口处的速度、压力、gamma 值，达到造波的目的，通过在waveMakerVelocityFvPatchVectorField.C 文件中加入# include “waveMakerVelocity.H ”，来调用waveMakerVelocity.H 文件。

waveMakerVelocity.H 文件内定义了多种类型的波
● 相当于定义了新的速度边界条件waveMakerV elocity ，在case 的定义中直接调用自定义
的边界条件 ● 有
OpenFOAM
自带的边界条件进行改写，这里根据
OpenFOAM-1.6/src/finiteV olume/fields/fvPatchFields/inletOutlet 进行改写。

(3) 线性波的数值模拟结果
建立一个长45米，高1.2米，水深为0.6米的二维水池模型，用速度入口式生产一个波高为0.02米，周期为1s 的线性波。

水池中生成的波面如图8、图9。

生成的波面图
图8 线性波的波面图
图9 距离造波边界处1个、2个、4个、6个波长处的波面时历曲线
4.2 如何在OpenFOAM 求解器中加海绵层，实现消波的目的
为了更有效地解决数值水槽末端的波浪反射问题，在VOF 模型中采用了海绵层阻尼消波的方法。

其原理是将某种形式的衰减系数加入到控制方程之中，使得在海绵层内部的流体满足修正后的时均运动方程：
u x z
u
y u x u v x p g z u w y u v x u u t u x ⋅-∂∂+∂∂+∂∂+∂∂-=∂∂+∂∂+∂∂+∂∂)()(1222222μρ v x z v y v x v v y p g z v w y v v x v u t v y ⋅-∂∂+∂∂+∂∂+∂∂-=∂∂+∂∂+∂∂+∂∂)()(1222222μρ （22） w x z
w y w x w v z p g z w w y w v x w u t w z ⋅-∂∂+∂∂+∂∂+∂∂-=∂∂+∂∂+∂∂+∂∂)()(1222222μρ 模型中的海绵层采用线性衰减形式
当0x x ≥时，0
10
)(x x x x x --=θ
μ，0x 为海绵层的起点，1x 为水池的长度；
当0x x ≤时，0)(=x μ。

OpenFOAM 对方程组（22）写成如下代码进行计算：
fvVectorMatrix UEqn
• (
• fvm::ddt(rho, U)
• + fvm::div(rhoPhi, U)
• - fvm::laplacian(muf, U)
• - (fvc::grad(U) & fvc::grad(muf))
• + U*rho*thetaField // sponge layer 海绵层
可见，在OpenFOAM中实现消波非常的简单，只要修改求解方程的外层代码即可。

增加了海绵层之和，在不同时刻的波面图，如图10。

图10 T=30s/40s/50s/60s时的水池的波面图
参考文献：
1. OpenFOAM:UserGuide,ProgrammerGuide
2. /forum/messages
3. Hassan Hemida.OpenFOAM tutorial:Free surface tutorial using interFOAM and rasInterFOAM,2008
4. Onno Ubbink. Numerical prediction of two fluid systems with sharp interfaces. Doctoral thesis for University of London, Imperial College,1997
5. Henrik Rusche. Computational Fluid Dynamics of Dispersed Two-Phase Flows at High Phase Fractions. Doctoral thesis for University of London,Imperial College,2002。