版本对应:

OpenFOAM-11中的incompressibleFluid模块,对应OpenFOAM-10之前的simpleFoam与pimpleFoam。OpenFOAM-10中的稳态瞬态求解器单独处理,OpenFOAM-11中都整合在了incompressibleFluid模块。

CFD: 不可压 + 瞬态

Issa在1986年提出的PISO算法,最重要的特点在于其为一个非迭代的瞬态压力速度耦合算法。PISO算法适用于不可压缩流以及可压缩流。OpenFOAM中的瞬态PISO算法与原始PISO算法大同小异。本文讨论的算法主要取自OpenFOAM-11中的incompressibleFluid模块。需要注意的是,incompressibleFluid模块可以同时处理瞬态、稳态,本文仅仅关注瞬态算法。

在老版本的OpenFOAM中存在3个瞬态不可压缩单相流求解器:icoFoam、pimpleFoam、pisoFoam。icoFoam为一个demo式的求解器,主要便于初学者理解代码,其不可调用湍流模型。pisoFoam为一个纯正的瞬态不可压缩单相流求解器。pimpleFoam可以调用pimple算法(将在后续介绍)。在OpenFOAM-11中,三者不再进行区分。incompressibleFluid模块可以涵盖所有的瞬态不可压缩压力速度耦合算法。

另外需要注意的是,OpenFOAM-11中的incompressibleFluid模块可以进行DNS模拟。其可用于模拟层流流动或流场剪切力引致的湍流准直接模拟 (quasi-DNS),也被称之为准DNS。在这里需要注意的是,得益于谱方法以及有限差分法高阶精度的易用性,DNS计算通常采用谱方法或有限差分法。有限体积法由于本身在计算通量时引入的中点插值法则,大大限制了有限体积法的精度。然而Komen等通过研究圆管内的湍流流动,认为在网格分辨率足够的情况下,使用有限体积法同样可以获得高质量的平均速度以及脉动速度流场 。Komen等进一步的使用有限体积法进行准DNS计算,并与有限体积法大涡模拟进行标定,认为大涡模拟方法会引致一定的数值耗散。Axtmann and Rist也使用OpenFOAM进行了准DNS直接模拟,并研究了编译器与并行计算的特性。Liu等通过icoFoam对开口顶盖驱动流进行了准直接模拟 ,进一步验证了使用OpenFOAM进行准直接模拟的可行性。

控制方程与算法

给定瞬态的N-S方程:

(1)\[ \nabla\cdot\bfU=0, \]
(2)\[ \frac{\p\bfU}{\p t}+\nabla \cdot (\mathbf{U}\mathbf{U})-\nabla \cdot(\nu \nabla \mathbf{U})=-\nabla p, \]

在这里我们不考虑湍流粘度的影响。首先对方程(2)中的时间项进行对速度\(\bfU\)关于时间\(t\)的欧拉全隐离散有 :

(3)\[ \int \int {\frac{{\partial \mathbf{U}}}{{\partial t}}\mathrm{d}V\mathrm{d}t = \left( {\mathbf{U}_\mathrm{P}^* - \mathbf{U}_\mathrm{P}^t} \right)\;} {V_\mathrm{P}}, \]

对流项隐性离散:

(4)\[ \int \int {\nabla \cdot \left( {\mathbf{U}\mathbf{U}} \right)\mathrm{d}V\mathrm{d}t = \int \int {\mathbf{U}\mathbf{U}\cdot\mathrm{d}\bfS\mathrm{d}t} = } \sum {{{\left( {{\mathbf{U}^*}{\mathbf{U}^t}} \right)}_f}} \cdot\bfS_f\Delta t = \Delta t\sum {{F_f^t \mathbf{U}_f^*} } , \]

拉普拉斯项隐性离散:

(5)\[\begin{split} \begin{equation} \begin{split} {\int \int \nabla \cdot \left( {\nu \nabla \mathbf{U}} \right)\mathrm{d}V\mathrm{d}t = \int \int {\nu \nabla \mathbf{U}\cdot\mathrm{d}\bfS\mathrm{d}t = \sum {\left( {\nu \nabla {\mathbf{U}^*}} \right)} } _f}&\cdot\bfS_f\Delta t \\ &= \Delta t\sum { \nu (\nabla\mathbf{U}^*)_f\cdot\bfS_f } , \end{split} \end{equation} \end{split}\]

压力项显性离散:

(6)\[ \int \int \nabla p \mathrm{d}V\mathrm{d}t=\int\int p \mathrm{d}\bfS\mathrm{d}t=\Delta t\sum \left(p_f^t\bfS_f\right), \]

其中上标\(^t\)表示为当前的时间步(已知),上标\(^*\)表示预测步(待求),下标\(_f\)表示网格单元面上的值,\(V_\rP\)表示网格单元体积,\(\bfS_f\)表示网格单元的各个面的面矢量,\(\Delta t\)表示时间步长,\(F_f^t\)为通量,\(\nu\)为动力粘度。方程(5)\((\nabla\mathbf{U}^*)_f\cdot\bfS_f \)需要进一步处理,其中\(\nabla\mathbf{U}^*\)为定义在体心的量,即网格体心的速度梯度。\((\nabla\mathbf{U}^*)_f\)为定义在面心的量,即网格面心的速度梯度。

考虑方程(3)-(6),有:

(7)\[ \frac{{\mathbf{U}_\mathrm{P}^* - \mathbf{U}_\mathrm{P}^t}}{{\Delta t}}{V_\mathrm{P}} + \sum {{F_f^t \mathbf{U}_f^*} } = \sum { \nu (\nabla\mathbf{U}^*)_f\cdot\bfS_f } -\sum p_f^t\bfS_f. \]

需要注意的是,方程(2)中的左边第二项(对流项)是非线性的。在求解的过程中,要么选择非线性求解器,要么将对流项线性化。由于非线性求解器非常复杂,因此在OpenFOAM均采用线性化处理。具体的,在对方程(7)中的左边第二项对流项离散中,其中的通量\(F_f^t\)采用当前已知时间步的速度\(\mathbf{U}^t\)来计算,同时保留\(\mathbf{U}^*_\mathrm{P}\)为未知量。这种将高阶非线性项降为一阶线性项的过程即为线性化操作。线性化的一个问题即为通量速度的信息有所滞后。

方程(7)\(\mathbf{U}_f^*\)被定义为面上的预测速度。面速度需要从体心速度进行插值来获得,在此步可以引入各种插值格式。假设在均一网格上使用中心线性格式:

(8)\[ \mathbf{U}_f^* = \frac{{\mathbf{U}_\mathrm{P}^* + \mathbf{U}_\mathrm{N}^*}}{2}. \]
(9)\[ p_f^t = \frac{{p_\mathrm{P}^t + p_\mathrm{N}^t}}{2}. \]

同时,为使用紧致基架点防止数值震荡,\((\nabla\mathbf{U}^*)_f\cdot\bfS_f\)通常表示为面法向梯度与面矢量的模的乘积:

(10)\[ (\nabla\mathbf{U}^*)_f \cdot\mathbf{S}_f= \left((\nabla\mathbf{U}^*)_f\cdot\frac{\mathbf{S}_f}{|\mathbf{S}_f|} \right ) \cdot|\mathbf{S}_f| \]

其中\((\nabla\mathbf{U}^*)_f\cdot\frac{\mathbf{S}_f}{|\mathbf{S}_f|} \)表示面法向梯度,有

(11)\[ (\nabla\mathbf{U}^*)_f\cdot\frac{\mathbf{S}_f}{|\mathbf{S}_f|} = \frac{{\mathbf{U}_\mathrm{N}^* - \mathbf{U}_\mathrm{P}^*}}{{\left| \mathbf{d} \right|}}, \]

其中\(\mathbf{d}\)表示网格单元体心之间矢量差(距离矢量),下标\(_\mathrm{N}\)表示相邻网格单元的速度,下标\(_\mathrm{P}\)表示当前网格单元的速度。速度的面法向梯度在三维情况下为一个矢量(变量的面法向梯度与变量的阶数相同)。在这里需要注意的是,方程(11)只对矩形网格精准,网格非正交性增加,需要特定的数值处理提高其计算精准性。

注意:

紧致基架点与拓展基架点的英文被称之为compact stencil与extended stencil,有关基架点的详细描述可参考无痛苦NS方程笔记中有关laplacian(ϕ) = div(grad(ϕ))?的讨论。

将方程(10)(8)(9)代入到方程(7)

(12)\[ \frac{{\mathbf{U}_\mathrm{P}^* - \mathbf{U}_\mathrm{P}^t}}{{\Delta t}}{V_\mathrm{P}} + \sum { {F_f^t \frac{{\mathbf{U}_\mathrm{P}^* + \mathbf{U}_\mathrm{N}^*}}{2}} } = \sum { {\nu \left| \bfS_f \right|\frac{{\mathbf{U}_\mathrm{N}^* - \mathbf{U}_\mathrm{P}^*}}{{\left| \mathbf{d} \right|}}} } -\sum \frac{{p_\mathrm{P}^t + p_\mathrm{N}^t}}{2}\bfS_f, \]
(13)\[\begin{split} \left( {\frac{{{1}}}{{\Delta t}} + \frac{1}{V_\mathrm{P}}\sum { {\frac{{F_f^t}}{2}} + \frac{1}{V_\mathrm{P}}\sum { {\nu \frac{{\left| \bfS_f \right|}}{{\left| \mathbf{d} \right|}}} } } } \right)&\mathbf{U}_\mathrm{P}^* = \\ - \sum \frac{1}{V_\mathrm{P}}{ {\left( {\frac{{F_f^t}}{2} - \nu \frac{{\left| \bfS_f \right|}}{{\left| \mathbf{d} \right|}}} \right)\mathbf{U}_\mathrm{N}^*} } & + \frac{{{1}}}{{\Delta t}}\mathbf{U}_\mathrm{P}^t-\frac{1}{V_\mathrm{P}}\sum \frac{{p_\mathrm{P}^t + p_\mathrm{N}^t}}{2}\bfS_f. \end{split}\]

将上式简写为

(14)\[ {A_\mathrm{P}^t}\mathbf{U}_\mathrm{P}^* = -\sum {A_\mathrm{N}^t\mathbf{U}_\mathrm{N}^*} +S_\mathrm{P}^t-\frac{1}{V_\mathrm{P}}\sum \frac{{p_\mathrm{P}^t + p_\mathrm{N}^t}}{2}\bfS_f, \]

其中

(15)\[ A_\mathrm{P}^t= {\frac{{{1}}}{{\Delta t}} + \frac{1}{V_\mathrm{P}}\sum { {\frac{{F_f^t}}{2}} + \frac{1}{V_\mathrm{P}}\sum { {\nu \frac{{\left| \bfS_f \right|}}{{\left| \mathbf{d} \right|}}} } } } , \]
(16)\[ A_\mathrm{N}^t= \frac{1}{V_\mathrm{P}}\left({\frac{{F_f^t}}{2} - \nu \frac{{\left| \bfS_f \right|}}{{\left| \mathbf{d} \right|}}}\right) , \]
(17)\[ S_\mathrm{P}^t=\frac{{{1}}}{{\Delta t}}\mathbf{U}_\mathrm{P}^t. \]

注意:

如果单单看扩散项的系数,有\(A_\mathrm{P}^t=\frac{1}{V_\mathrm{P}}\sum { {\nu \frac{{\left| \bfS_f \right|}}{{\left| \mathbf{d} \right|}}} }\)\(-\sum A_\mathrm{N}^t=\sum \frac{1}{V_\mathrm{P}} { {\nu \frac{{\left| \bfS_f \right|}}{{\left| \mathbf{d} \right|}}} }\)。很明显\(A_\mathrm{P}^t=-\sum A_\mathrm{N}^t\)。此对应关系大量的出现在CFD教材中。对流项也有类似的特征,只不过不是相等的关系,而是相减的关系。

注意:

另外,上面的离散都可以被认为是内部面的离散,对于边界面,可能会导致不同分量离散出来的矩阵系数可能不同。单OpenFOAM里面的UEqn.A()只有一个。这方面的介绍请参考 这个链接。边界面的离散不影响理解整个算法流程。

PISO算法

可以看出在某个时间步内\(A_\mathrm{P}^t\)\(A_\mathrm{N}^t\)、和\(S_\mathrm{P}^t\)保持不变(注意\(A_\mathrm{N}^t\)\(\sum A_\mathrm{N}^t\bfU_\rN\)\(-\sum A_\mathrm{N}^t\bfU_\rN+S_\rP^t\)的区别)。求解方程(14)即可获得预测速度\(\mathbf{U}^*_\mathrm{P}\)。方程(14)即为OpenFOAM中的动量预测方程。需要注意的是,动量预测步骤并不是必须的,这主要和速度压力耦合求解策略有关。如果不调用动量预测步骤,则\(\mathbf{U}^*_\mathrm{P}=\mathbf{U}^t_\mathrm{P}\)

N-S方程求解的关键问题之一在于并没有压力的方程出现。仔细分析得知上文中求解动量方程获得了速度,但连续性方程并不是关于压力的方程。压力泊松方程的构建正式基于连续性方程而来,其通过对动量方程继续做散度并相加获得。考虑最终收敛的情况下,对连续性方程\(\nabla\cdot\bfU=0\)进行离散后的形式为

(18)\[ \sum \left(\mathbf{U}_{\mathrm{P},f}^{t+\dt} \cdot \bfS_f\right)=0. \]

如果能用压力表示方程(18)中的\(\mathbf{U}_{\mathrm{P},f}^{t+\dt}\),是不是就可以得出压力泊松方程了呢?答案是肯定的。首先,方程(14)可以写为:

(19)\[ {A_\mathrm{P}^t}\mathbf{U}_\mathrm{P}^{*}{\rm{ + }}\sum {A_\mathrm{N}^t\mathbf{U}_\mathrm{N}^{*}} =S_\mathrm{P}^t-\frac{1}{V_\rP}\sum p_f^{t}\bfS_f, \]

在当前时间步收敛情况下为:

(20)\[ {A_\mathrm{P}}^{t+\dt}\mathbf{U}_\mathrm{P}^{t+\dt}{\rm{ + }}\sum {A_\mathrm{N}^{t+\dt}\mathbf{U}_\mathrm{N}^{t+\dt}} = S_\mathrm{P}^{t+\dt}-\frac{1}{V_\rP}\sum p_f^{t+\dt}\bfS_f, \]

方程(20)为一个非线性方程组。使其线性化后可以写为:

(21)\[ {A_\mathrm{P}}^{t}\mathbf{U}_\mathrm{P}^{t+\dt}{\rm{ + }}\sum {A_\mathrm{N}^{t}\mathbf{U}_\mathrm{N}^{t+\dt}} = S_\mathrm{P}^{t}-\frac{1}{V_\rP}\sum p_f^{t+\dt}\bfS_f, \]

方程(20)相对于方程(21)存在一定的滞后。在\(\dt\)较小的情况下这被认为是较小的。若将方程(21)减去(19)有:

(22)\[ {A_\mathrm{P}}\mathbf{U}_\mathrm{P}'+\sum {A_\mathrm{N}\mathbf{U}_\mathrm{N}'} = -\frac{1}{V_\rP}\sum p_f '\bfS_f, \]

其中\(\mathbf{U}_\mathrm{P}'=\mathbf{U}_\mathrm{P}^{t+\dt}-\mathbf{U}_\mathrm{P}^{*}\)(其他同理)。在这里未调用任何略去临点的假定。对方程(21)移项有

\[ \mathbf{U}_\mathrm{P}^{t+\dt} = \mathbf{HbyA}^{t+\dt}_\mathrm{P} - \frac{1}{{{A_\mathrm{P}^t}}} \frac{1}{V_\rP}\sum p_f^{t+\dt}\bfS_f. \]
(23)\[ \mathbf{HbyA}^{t+\dt}_\mathrm{P} = \frac{1}{{{A_\mathrm{P}^t}}}\left( { - \sum {{A_\mathrm{N}}\mathbf{U}_\mathrm{N}^{t+\dt}} + S_\mathrm{P}^n} \right). \]

类似的,面上的速度可以通过下面的公式获得:

(24)\[ \mathbf{U}_{\mathrm{P},f}^{t+\dt} = \mathbf{HbyA}_f^{t+\dt} - \frac{1}{{{A_{\mathrm{P},f}^t}}} \left(\frac{1}{V_\rP}\sum p_f^{t+\dt}\bfS_f\right)_f. \]

将方程(24)代入到方程(18)有:

(25)\[ \sum \left( \mathbf{HbyA}_f^{t+\dt} - \frac{1}{{{A_{\mathrm{P},f}^t}}}\left(\frac{1}{V_\rP}\sum p_f^{t+\dt}\bfS_f\right)_f \right) \cdot \bfS_f =0, \]
(26)\[ \sum { {\mathbf{HbyA}_f^{t+\dt}} \cdot \bfS_f } = \sum { {\frac{1}{{A_{\mathrm{P},f}^t}}\left(\frac{1}{V_\rP}\sum p_f^{t+\dt}\bfS_f\right)_f} \cdot \bfS_f }, \]

方程(26)即下述方程的离散形式:

(27)\[ \nabla \cdot \mathbf{HbyA}^{t+\dt} = \nabla \cdot\left(\frac{1}{A} \nabla p^{t+\dt}\right). \]

方程(27)即为压力泊松方程。若有\(\mathbf{HbyA}^{t+\dt}\),则可求得收敛的压力。需要注意的是,\(\nabla \cdot\left(\frac{1}{A^t} \nabla p^{t+\dt}\right)\)可以通过先求梯度再求散度进行,也可以直接通过拉普拉斯积分进行。方程(26)右端即先求梯度再求散度的方法,这种方法会引起数值震荡。若通过拉普拉斯并进行高斯积分进行离散,则会调用紧致格式防止数值震荡 (Li, 2019)

可以看出,最后求解的速度\(\mathbf{U}_\mathrm{P}^{t+\dt}\)和压力\(p^{t+\dt}\)应该符合方程(14)(27),分别对应动量方程和连续性方程。然而目前,我们只有通过动量预测步骤求出来的\(\mathbf{U}_\mathrm{P}^*\)\(\mathbf{HbyA}_\mathrm{P}^{t+\dt}\)\(p^{t+\dt}\)均为未知的。压力泊松方程不能求解。1986年Issa提出了PISO算法可以解决这个问题 (Issa, 1986)。PISO算法在提出时主要针对不可压缩非稳态计算,其为一种在时间步内非迭代的算法,在随后也被拓展用于稳态问题以及可压问题中。PISO算法通过对当前时间步内的多次修正来获得最终的\(\mathbf{U}_\mathrm{P}^{t+\dt}\)\(p^{t+\dt}\)。依据PISO算法,在方程(22)的基础上引入略去临点影响的假定有:

(28)\[ {A_\mathrm{P}}\mathbf{U}_\mathrm{P}' = -\frac{1}{V_\rP}\sum p_f '\bfS_f, \]

将方程(28)(19)相加有

(29)\[ {A_\mathrm{P}^t}\mathbf{U}_\mathrm{P}^{**}{\rm{ + }}\sum {A_\mathrm{N}^t\mathbf{U}_\mathrm{N}^{*}} = S_\mathrm{P}^t-\frac{1}{V_\rP}\sum p_f^*\bfS_f, \]

其中\(\mathbf{U}_\mathrm{P}^{**}=\mathbf{U}_\mathrm{P}^{*}+\mathbf{U}_\mathrm{P}'\)\(p_f^{*}=p_f^{t}+p_f'\)。相比方程(21),方程(29)由于忽略了临点的影响并不是精准的。

注意:

如果不进行动量预测步计算,则\(\mathbf{U}_\mathrm{P}^{**}=\mathbf{U}_\mathrm{P}^{n}+\mathbf{U}_\mathrm{P}'\)

对方程(29)进行移项有

(30)\[ \mathbf{U}_\mathrm{P}^{**} = \mathbf{HbyA}^{*}_\mathrm{P} - \frac{1}{{{A_\mathrm{P}^t}}}\frac{1}{V_\rP}\sum p_f^*\bfS_f. \]
(31)\[ \mathbf{HbyA}^{*}_\mathrm{P} = \frac{1}{{{A_\mathrm{P}^t}}}\left( { - \sum {{A_\mathrm{N}^t}\mathbf{U}_\mathrm{N}^{*}} + S_\mathrm{P}^t} \right). \]
(32)\[ \mathbf{U}_{\mathrm{P},f}^{**} = \mathbf{HbyA}_f^{*} - \frac{1}{{{A_{\mathrm{P},f}^t}}}\left(\frac{1}{V_\rP}\sum p_f^*\bfS_f\right)_f. \]

方程(32)中可参考方程(25)-(27)的步骤组建压力泊松方程,即

(33)\[ \nabla \cdot\left(\frac{1}{A^t} \nabla p^{*}\right)=\nabla \cdot \mathbf{HbyA}^{*} \]

对方程(27)求解后有\(p^{*}\)。将\(p^*\)回代到方程(30)的时候,有\(\mathbf{U}_{\mathrm{P}}^{**}\)。这里的\(\mathbf{U}_{\mathrm{P}}^{**}\)满足连续性方程。但由于方程(29)忽略了临点的影响,\(\mathbf{U}_{\mathrm{P}}^{**}\)\(p^{*}\)并不严格满足动量方程。因此在实际操作中,需要以计算后的\(\mathbf{U}_{\mathrm{P}}^{**}\)继续组建\(\bfHbyA^{**}\),然后继续求解压力方程获得新的压力场\(p^{**}\),然后进一步获得新的速度场\(\mathbf{U}_{\mathrm{P}}^{***}\)。这即为时间步内的第二次迭代。当然也可以在时间步内继续进行迭代。

因此对于PISO算法,每个时间步内的迭代可以表示为:

  1. 依据初始条件求解预测速度\(\mathbf{U}^*\)。通过方程(32)组建\(\mathbf{HbyA}^*\)

  2. 求解方程(33)求解压力\(p^*\)

  3. 通过方程(32)依据压力\(p^*\)以及\(\mathbf{HbyA}^*\)更新速度有\(\mathbf{U}^{**}\)

  4. 回到第二步迭代求解几次。

PISO算法在时间步内的迭代一般只需要3步即可。相对于稳态SIMPLE算法的若干次迭代,瞬态PISO算法的区别主要在于:

  1. 瞬态PISO算法的时间步长要保证库朗数足够小,因此每一步的流场变动相对来说并不是很激进。因此,在求解方程(33)与方程(29)的时候,可以将残差设定为非常小的值。保证动量方程与连续性方程都是收敛的。这样,误差仅仅来自于略去临点的误差。大大增加每个时间步内的收敛速度。

  2. 瞬态PISO算法不需要对压力场进行松弛。但可以对速度、湍流变量的矩阵进行少量的松弛,过量的松弛会导致时间步进不精准。因此,PISO算法时间步内的算法可以理解为没有松弛的SIMPLE算法。

  3. 理论上,瞬态PISO算法也可以与SIMPLEC相结合,形成单时间步内迭代更快的算法。但由于瞬态求解器的库朗数的严格限制。使用SIMPLE算法就可以在每时间步内保证3次以内收敛,因此实践上,并没有使用SIMPLEC算法的实际意义。

PIMPLE算法

PIMPLE算法是OpenFOAM基金会提出特有的计算方法。并没有见过在文献里面当做专门的算法而发表的文章。PIMPLE算法更像是一种算法的融合,而非重新开发。PIMPLE算法有一些非常吸引人的特性:

  • 调用库郎数远大于1的时间步长且同时进行瞬态步进;

  • 在时间精度不重要的时候,获得一种拟瞬态的稳态解;

PIMPLE算法将PISO算法与SIMPLE算法融合,获得大时间步长步进的计算能力。在瞬态PISO算法中(例如上文讨论的icoFoam),需要保证时间步长要保证库朗数足够小,这样才能不需要加松弛,且保证每个时间步内计算的结果是收敛的。在PIMPLE算法中,可以调用很大的时间步长。在这种情况下,每个时间步长内的流场或许会变化的很激进容易发散,因此在PIMPLE算法中,可以在每个时间步内添加松弛。可以把PIMPLE理解为,在较大时间步长下,每个时间步内去添加松弛来获得收敛解,然后向后进行推进。PIMPLE算法将低松弛去掉,并关闭PIMPLE循环,可以完全的演变为瞬态PISO算法。

若使用PIMPLE算法,建议使用比较高的时间精度格式来减少高库郎数导致的误差。然而在全局库郎数比较大的情况下,时间历史步进的不精准是难以避免的,同时松弛因子对结果的影响也会加大。因此计算结果如果需要非常高精度的时间步进历史,需要使用PISO算法。若对时间步进要求存在要求且网格大小分布非常不均一,可以使用PIMPLE算法增大库郎数限制,但要关注全局库郎数的变化。若全局库郎数的变化特别大。不可否认的会牺牲流场变化的时间精度。若对时间步进的历史过程不关心,仅仅关注最后稳态的情况,且同时算例存在收敛解。完全可以使用SIMPLE算法来进行。也可以通过PIMPLE算法来增加到足够的库郎数来进行时间推进。也可以通过局部时间步Local Time Step的方式来进行。

下面是一个典型的PIMPLE设置。其中的nOuterCorrectors用来控制PIMPLE循环次数,也被称之为外循环。表示在每个时间步长内进行100次迭代(可以理解为100次SIMPLE迭代)。nCorrectors表示内循环。表示每次外循环内的压力方程求解次数,可以理解为PISO的迭代次数。在PIMPLE算法中,一般内循环设置为1次即可,即不需要多次PISO内循环。

因此,nOuterCorrectors设置为1的时候,PIMPLE完全的变为了PISO算法。outerCorrectorResidualControl用来控制PIMPLE外循环的收敛标准(类似SIMPLE算法的收敛标准)。比如本算例设置100次外循环,但如果满足相应的收敛标准,则停止迭代。对于一个存在近期时间稳态的算例,在求解后期,PIMPLE外循环可以降低为10次以下甚至更少,且保持非常高的库郎数递进。relaxationFactors用来设置外循环的松弛因子。其设置完全可以参考SIMPLE算法。

PIMPLE
{
    momentumPredictor no;
    nOuterCorrectors 100;
    nCorrectors     1;
    nNonOrthogonalCorrectors 0;

    outerCorrectorResidualControl
    {
        "(h|p)"
        {
            relTol          0;
            tolerance       0.0001;
        }
    }
}

relaxationFactors
{
    fields
    {
        p      0.3;
        pFinal   1;
    }
    equations
    {
        "U|k|epsilon|h"     0.7;
        "(U|k|epsilon|h)Final"   1;
    }
}

关键代码

incompressibleFluid模块中的速度方程通过下述代码来进行组建:

fvVectorMatrix UEqn
(
    fvm::ddt(U)
  + fvm::div(phi, U)
  - fvm::laplacian(nu, U)
);

其中fvm::ddt(U)对应方程(3),这里的fvm表示隐形离散,反之,若写为fvc则表示显性离散。fvm::div(phi, U)对应方程(4),注意其中的phi,在OpenFOAM中,phi大部分情况下表示为通量。进一步的,在不可压缩流体中,phi表示定义在面上的体积通量\(\bfU_f\cdot\bfS_f\),在可压缩流体中,phi表示定义在面上的质量通量\(\rho_f\bfU_f\cdot\bfS_f\)fvm::laplacian(nu, U)则对应方程(5)

方程(15)中定义的\(A_\rP^t\)为一个volScalarField,其可以通过UEqn.A()来表示。\(A_\rP\)的倒数,在代码中被表示为rAU,其代码如下:

volScalarField rAU(1.0/UEqn.A());

在方程(23)中,定义了\(\mathbf{HbyA}\),其可以理解为进行一次迭代后获得的一个速度,更激进的,可以把HbyA理解为速度中间量。方程(23)右侧括号内的\({ - \sum {{A_\mathrm{N}^t}\mathbf{U}_\mathrm{N}^{t+\dt}} + S_\mathrm{P}^t}\)可以用代码UEqn.H()来表示。因此进一步的,\(\mathbf{HbyA}\)可以表示为

volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();

其中第一行仅仅为声明场,第二行为对其进行赋值。在获得\(\mathbf{HbyA}\)之后,方程(27)可以用来求解压力,相应的代码如下:

fvScalarMatrix pEqn
(
    fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
);

如果进一步的看foamRun代码,

{
    solver.moveMesh();//网格变形
    solver.fvModels().correct();//源项添加
    solver.prePredictor();//关键步1
    solver.momentumPredictor();//关键步2
    solver.thermophysicalPredictor();//关键步3
    solver.pressureCorrector();//关键步4
    solver.postCorrector();//后处理
}

其中的:

  1. 关键步1不执行;

  2. 关键步2初始速度矩阵,若需要动量预测,则进行动量预测;

  3. 关键步3不执行;

  4. 关键步4进行压力求解;