OpenFOAM的不可压缩流算法

对不可压缩流的几个经典算法及其在OpenFOAM中的实现进行了梳理。

主要针对OpenFOAM 5.x版本

SIMPLE, PISO和PIMPLE算法

简介

均质-不可压缩-无体积力-常粘性的NS方程为:

OpenFOAM和其他的CFD软件常用来解它的算法有:

  • SIMPLE: Semi-Implicit Method of-Pressure Linked Equations,用于稳态计算
  • PISO:= Pressure Implicit Split Operator,用于瞬态计算,可以使用的Courant数往往小于11
  • PIMPLE:= Merged PISO–SIMPLE,可以使用Courant数>>1

但是找文献和教科书对的时候你会发现很多时候算法和OpenFOAM里的实现对不上。比如wiki和cfd-online上给的SIMPLE算法描述就是(Wiki: SIMPLE Algorithm):

  1. 修正边界条件;
  2. 计算速度和压力梯度;
  3. 求解离散动量方程( discretized momentum equation),求得体心临时速度场;这一步是显式的
  4. 用体心临时速度场插值计算面心质量通量场;
  5. 计算压力修正方程(pressure correction equation),得到体心压力修正场;这一步是隐式的
  6. 对体心压力场进行亚松弛修正;
  7. 对压力场修正边界条件;
  8. 修正面心质量通量场;
  9. 修正体心速度场;
  10. 更新密度;

但是明显OpenFOAM中的simpleFoam并没有求解过所谓的修正方程,而是自己玩了一个Picard迭代9直接求的全量方程(你会发现OpenFOAM似乎从来就没求过修正量的Newton迭代,个人理解是因为OpenFOAM玩Newton迭代时边界条件不是很方便处理。)

于是我做了一点儿考古和分析。

SIMPLE算法

SIMPLE: Semi-Implicit Method of Pressure Linked Equations,这个方程的名字就满是问题,首先什么是半隐式方法(Semi-Implicit Method),什么又是压力联系方程(Pressure Linked Equations)?在原始文献2 中是没有用SIMPLE这个名字的。搜索文献也很难找到压力联系方程的说法。我至今没有找到这个术语最开始的出处,查到最早的是他们自己1973年的文献5 是叫SIMPLE了,但似乎大家都是以讹传讹地引用1972年的那篇文献。

根据一些参考文献3 的分类法,SIMPLE, PISO都属于一类叫压力修正格式(Pressure Correction Schemes)的算法,其特征是构造一个压力场来满足不可压缩条件,再得到一个速度场,至于这个速度场能不能满足新得到的压力场,呵呵,要是能满足就不用再继续迭代了。由于求解的动量方程已经是线性化之后的了,所以你把这个速度对应的面心通量场再去线性化一次动量方程,如此往复才能最后得到稳态解。

压力联系方程之到底是哪个方程

查来查去,发现不可压缩CFD解的主要是以下方程:

  • 压力修正方程(Pressure Correction Equation, PCE)4 : 两次迭代间压力需要的修正量的方程
  • 压力泊松方程(Pressure Poisson Equation, PPE) 6: 动量方程求散度,消去速度散度项

    稳定化的一阶欧拉时间离散时的形式。

​ 速度场保持无散时的形式。

  • 压力联系方程(Pressure Linked Equation, PLE):原始文献并没有说什么是PLE,参考另一篇文献7的说法,压力联系方程的形式如下(符号定义参考文献):

可见,所谓PLE非常诡异,它不是个微分方程,而是一个代数方程。对应到OpenFOAM里的simpleFoam则是8

同时,对比PLE和PPE,可以发现,二者都是把某种形式的动量方程带入了质量守恒方程得到的,但不同之处在于

  1. PLE是把线性化离散动量方程中得到速度表达式带入了质量守恒方程,得到$\nabla\cdot(\frac 1 a \nabla p) = \dots $。
    1. 大概是这么个过程:
    2. 线性化离散动量方程记为 $M\cdot u = \nabla p $, $M$是一个矩阵。
    3. 先LDU加法分解: $M\cdot u = (D+L+U) \cdot u = \nabla p$
    4. 然后进行Jacobi迭代:$ D \cdot u = \nabla p - (L+U) \cdot u^0$
    5. 左边: $D$记为$A$ (实际上是对角系数除以单元体积,以保持量纲一致性), 因为是对角阵,所以可以把$D^{-1}=A^{-1}$记为 $\frac 1 A$ ,OpenFOAM里实现为一个几何场而不是矩阵fvMatrix
    6. 右边:$(L+U)\cdot u^0$记为 $H$
    7. 这样$u = \frac 1 A \nabla p - \frac H A$
    8. 然后再利用$\nabla \cdot u = 0 $
    9. 得到$\nabla\cdot(\frac 1 A \nabla p ) = \nabla \cdot \frac H A$ ,这就是PLE,也是实际OpenFOAM中求解的方程。
    10. 英语里$\frac H A$ 可以读作”H by A”
    11. OpenFOAM里的$\frac 1 A$ 对应变量为rUA
  2. PPE是把动量方程直接求散度得到的,所以是$\nabla\cdot\nabla p=\dots$。

FAQ

  • 问:为啥要费劲心力去算$\mathbf{HbyA}$,再算PLE,直接离散PPE或者质量守恒方程不行么?
  • 答:我也不知道

    • 但是玩FEM的家伙们就是这么干的(当然实际上他们也有很多流派,但是他们更多的是玩函数空间,LBB条件而不是玩方程系数,玩系数太low)
    • 似乎只有玩FVM的人才玩PLE。
    • 而且投影法10那个流派的FVM也是玩的PPE,参考Chorin’s_projection_method#Chorin’s_projection_method)。

    • 而且为了玩PLE,作为使用OpenFOAM的同位网格玩家,还必须引入Rhie-Chow插值这种大坑,虽然OpenFOAM似乎轻巧地避开了这个坑。

Semi-Implicit-Method半隐式方法

我认为之所以叫半隐式方法,是因为:

  • 求解过程中,求速度预测步的方程对于速度是隐式的,但对于压力是固定的显式离散,所以称之为半隐式。
  • 求解压力方程时,压力是隐式离散的,但是此时的速度是固定不变的。
  • 所以本质上这个半隐式等价于解耦算法。
  • 文献5 摘要中提到,SIMPLE的半隐式是和SIVA(SImultaneous Variable Adjustment)全隐算法相对的提法,其特点是每次迭代要同时更新周围的速度和压力。

SIMPLE算法因此具有一些特性:

  • 解耦计算,内存占用小。
  • 用于稳态计算,所以timeStep无意义,但通常设置timeStep=1,这样解算输出时间等于时间步数;
  • 由于没有时间项来稳定化计算,需要加入松弛因子(在矩阵上加在和时间项几乎相同的位置),但是这个松弛因子在实现时搞不好会对最终稳态解产生影响。

关于consistent

对于SIMPLE-C算法16而言,比SIMPLE算法只多了一个步骤:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
//applications/solvers/incompressible/pimpleFoam/pEqn.H +16
if (pimple.consistent())
{
rAtU = 1.0/max(1.0/rAU - UEqn.H1(), 0.1/rAU); //有max!
phiHbyA +=
fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf();
HbyA -= (rAU - rAtU())*fvc::grad(p);
}


//applications/solvers/incompressible/simpleFoam/pEqn.H +10
if (simple.consistent())
{
rAtU = 1.0/(1.0/rAU - UEqn.H1()); //没有max!
phiHbyA +=
fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf();
HbyA -= (rAU - rAtU())*fvc::grad(p);
}

值得注意的是,pimpleFoam和simpleFoam的consistent选项对rAtU的计算公式是有区别的!

其中修改了的变量有:

  • rAtU

    • 1.0/rAU - UEqn.H1()的实际含义应该是$[1/(\frac{1}{D})-(-(L+U)\cdot \mathbf{1})] = (L+D+U)\cdot \mathbf{1}$,也就是矩阵系数的行和,中间最多差一个网格体积的因子。
    • 0.1/rAU的含义应该是$\frac{1}{10} D$
    • 由于系数矩阵元素的大小往往是对角项为负,非对角项为正,其行和为略小的负数(对于laplace方程,某些行和可能为0,但是对于动量方程,应该不为零)。所以为了避免除以0,pimpleFoam采用了max函数的技巧。
    • 而对于simpleFoam,只取了第一项的倒数。
  • phiHbyA

    • 按照rAtU-rAU的差值,和grad(p)对界面流量进行更新

    • 把这项代入到pEqn的表达式中,可以发现和原来是一样的

    • 1
      2
      3
      4
      5
      6
      7
      8
      9
      10
      11
      12
      13
      14
      15
      16
      17
      18
      19
      //原来的方程,如pisoFoam, icoFoam就采用了这个方程
      fvScalarMatrix pEqn
      (
      fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
      );

      //simpleFoam, pimpleFoam采用的方程,将phiHbyA改记为phiHbyA2
      fvScalarMatrix pEqn
      (
      fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA2)
      );

      //两者相减
      fvScalarMatrix pEqn
      (
      fvm::laplacian(rAtU()-rAU, p)
      == fvc::div(fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf())
      );
      //可见二者是一致的。
  • HbyA

    • 同上,按照rAtU-rAU的差值,对临时体心速度进行更新。
    • 便于后面用U = HbyA - rAtU()*fvc::grad(p);去更新真正的体心速度!
  • 至此,正确性没有问题,但是为啥rAtU要这么改?!。

    • 对于pimpleFoam中的rAtUmax()函数在分母,应该是取了较为稳定的一组系数,并避免了除以0。

    • 对于simpleFoam而言,rAU是矩阵系数中的对角部分,而rAtU则是取了矩阵每一行系数的和。

      • 参考原始文献16 ,作者就是这么干的!

      • 本质上你可以任意的分解:

      • 所以对于consistent=false的情形:

      • 对于consistent=true的情形:

      • 根据wiki上的说明,如果压力-速度耦合是解收敛的主要问题,SIMPLE-C算法可以加速收敛,SIMPLE-C比SIMPLE要快大概是120%-130%左右。

      • 如果压力-速度耦合不是收敛的主要障碍,SIMPLE-C效率与SIMPLE差不多。

关于H1和H的区别

  • 意义不一样:假设$M=(L+D+U)$,初值为$x^0$
    • 则$H_1 = H(x^0=\mathbf1)$
    • $H(x^0)=-(L+D)\cdot x^0$
    • $\mathbf 1$表示全是1的向量。
  • 都是定义在lduMatrix中的,但是H1是不带参数的,而H是带参数的。
  • 实现的地方也不一样,一个在lduMatrixATmul.C,另一个在lduMatrixTemplates.C
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
//H1()
//src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixATmul.C +298
Foam::tmp<Foam::scalarField > Foam::lduMatrix::H1() const
{
tmp<scalarField > tH1
(
new scalarField(lduAddr().size(), 0.0)
);

if (lowerPtr_ || upperPtr_)
{
scalarField& H1_ = tH1.ref();

scalar* __restrict__ H1Ptr = H1_.begin();

const label* __restrict__ uPtr = lduAddr().upperAddr().begin();
const label* __restrict__ lPtr = lduAddr().lowerAddr().begin();

const scalar* __restrict__ lowerPtr = lower().begin();
const scalar* __restrict__ upperPtr = upper().begin();

const label nFaces = upper().size();

for (label face=0; face<nFaces; face++)
{
H1Ptr[uPtr[face]] -= lowerPtr[face];
H1Ptr[lPtr[face]] -= upperPtr[face];
}
}

return tH1;
}

//H(psi)
//src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixTemplates.C +34
template<class Type>
Foam::tmp<Foam::Field<Type>> Foam::lduMatrix::H(const Field<Type>& psi) const
{
tmp<Field<Type>> tHpsi
(
new Field<Type>(lduAddr().size(), Zero)
);

if (lowerPtr_ || upperPtr_)
{
Field<Type> & Hpsi = tHpsi.ref();

Type* __restrict__ HpsiPtr = Hpsi.begin();

const Type* __restrict__ psiPtr = psi.begin();

const label* __restrict__ uPtr = lduAddr().upperAddr().begin();
const label* __restrict__ lPtr = lduAddr().lowerAddr().begin();

const scalar* __restrict__ lowerPtr = lower().begin();
const scalar* __restrict__ upperPtr = upper().begin();

const label nFaces = upper().size();

for (label face=0; face<nFaces; face++)
{
HpsiPtr[uPtr[face]] -= lowerPtr[face]*psiPtr[lPtr[face]];
HpsiPtr[lPtr[face]] -= upperPtr[face]*psiPtr[uPtr[face]];
}
}

return tHpsi;
}

template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::lduMatrix::H(const tmp<Field<Type>>& tpsi) const
{
tmp<Field<Type>> tHpsi(H(tpsi()));
tpsi.clear();
return tHpsi;
}

其他

wiki和cfd-online上的SIMPLE算法是和simpleFoam不一致的。正确的描述可以参考openfoamwiki上的描述:the SIMPLE Algorithm

simple算法在OpenFOAM中的实现也是和版本相关的,比如foam-extend 4.0的simpleFoam里Jasak就采用了一种新的算法来避免时间步和松弛因子对解的影响。OpenFOAM.org和OpenFOAM.com的似乎还没变化。

此外,SIMPLE算法还有一些衍生算法:

  • SIMPLEC: SIMPLE-Consistent
  • SIMPLER: SIMPLE-Revised

Rhie-Chow插值

Rhie-Chow插值,有人又叫压力加权插值方法(pressure-weighted interpolation method, PWIM),但是其实这个插值并不是用压力来加权的。

刚开始我以为这个插值的作用是来插值计算面心压力的,后来发现其实它是用来计算面心速度的插值方法。

本质上它等效于引入了4阶耗散项进入压力修正方程中,有利于计算的稳定化:

按照Rhie-Chow插值,均匀一维网格的面心速度表达式会含有P的三次导数成分。

带入连续性方程之后,就会含有四次导数成分,成为稳定化压力方程的高阶耗散项,从而可以平抑压力波动:

详细的说明可以看参考文献13

而且玩高阶导数的人都应该知道,压力梯度都特别大的时候,高阶项会更大,所以有时候Rhie-Chow插值也会出问题。

面心速度的作用

(注:面心速度($\mathbf{u}_f$)和界面通量(程序里常用phi表示 )仅相差一个面积,所以实际上差不多是同一回事儿。)

根据文献11 ,面心速度有三个作用:

  1. 格式系数作用:构成关于体心速度的矩阵方程UEqn时,界面通量是系数的重要组成部分,系数不同,方程就不同,解就不同;
  2. 耦合作用:耦合速度和压力,在simple算法构造关于体心压力的矩阵方程时,需要用到界面通量;
  3. 质量守恒作用:界面通量应满足质量守恒。

按照我的理解,在OpenFOAM中,这些作用在程序中分别由不同的变量承担了,并不是同一个面心速度了。

  1. 格式系数作用:由于采用的是Picard迭代,所以OpenFOAM组成速度方程UEqn的矩阵系数是采用直接按fvScheme指定格式插值得到的界面通量phifvm::div(phi,U)生成的(准确地说在simpleFoam中第一次迭代用的插值phi,后面用的pEqn.flux()phiHbyA构造的phi)。在simpleFoam中如以下代码所示。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
//UEqn.H
tmp<fvVectorMatrix> tUEqn
(
fvm::div(phi, U) //这里用的界面通量是phi
+ MRF.DDt(U)
+ turbulence->divDevReff(U)
==
fvOptions(U)
);
fvVectorMatrix& UEqn = tUEqn.ref();

//src/finiteVolume/cfdTools/incompressible/createPhi.H
surfaceScalarField phi
(
IOobject
(
"phi",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
fvc::flux(U) //调用的是fvc::flux()
);

//src/finiteVolume/finiteVolume/fvc/fvcFlux.C
Foam::tmp<Foam::surfaceScalarField> Foam::fvc::flux
(
const volVectorField& vvf
)
{
return scheme<vector> // scheme<vector>(mesh,name) 返回 tmp<surfaceInterpolationScheme>对象
(
vvf.mesh(),
"flux(" + vvf.name() + ')'
)().dotInterpolate(vvf.mesh().Sf(), vvf); //插值并点乘
}

//pEqn.H
while (simple.correctNonOrthogonal())
{
//...
if (simple.finalNonOrthogonalIter())
{
phi = phiHbyA - pEqn.flux(); //更新phi
}
}
  1. 而耦合作用反映在phiHbyA这个界面通量中,它最终构成了pEqn中的源项部分,在OpenFOAM的simpleFoam求解器中,phiHbyA是用fvc::flux(HbyA) 构造的。
    1. 但是值得注意的是,simpleFoam只构造出了不含压力梯度贡献的临时速度$u^*$ = HbyA,并且没有直接使用它,而是用它构造了phiHbyA
    2. simpleFoam没有在pEqn求解之前构造含有压力梯度贡献的U,这是在求解完压力方程后,用U = HbyA - rAtU()*fvc::grad(p) 构造的。
    3. 所以dyfluid上说在OpenFOAM中,并没有直接采用Rhie-Chow插值的原始步骤。相反的,其通过对拉普拉斯项的离散巧妙的获得一种近似Rhie-Chow插值的原理。
    4. 因为原始步骤是通过$\nabla\cdot U=0$构造pEqn,而实际上因为pEqn是通过phiHbyArUA构造的,所以跳过了原始的Rhie-Chow插值。
    5. 而因为fvm::laplacian(rUA,p)的实现已经是利用了周围相邻单元的压力,并没有原始构造方式中压力的奇偶失耦(odd-even decoupling)情况,所以认为是实现了Rhie-Chow插值的效果。
    6. 按照参考文献18 的搞法(Listing 15.18),Rhie-Chow插值似乎是用于产生phiHbyA的,但是由于那本书严格地遵守了解PCE的思路,所以还不能完全相互比较。
  2. 质量守恒作用其实是靠pEqn实现的PLE完成的,simpleFoam除了在continuityErrs.H之外不直接检查质量守恒,检查时也是用的fvc::div(phi)而不是速度U

PISO算法

PISO = Pressure Implicit Split Operator。同样的问题,什么是压力隐式(pressure implicit),什么是算子分裂(split operator)。就icoFoam来看,我的理解是:

  • 压力隐式(pressure implicit):求解的压力方程是隐式的(废话);

  • 算子分裂(split operator):把$\frac {dx}{dt}|_t^{t+\Delta t}=O(x) = (O_1+O_2)(x) $ 近似成$ \frac{dx}{dt}|_t^{t+\Delta t} = O_1(x’), \frac{dx’}{dt}|_t^{t+\Delta t}=O_2(x)$来求解

    • 原始文献15中其实也就是解耦速度和压力的意思:

    • The principle is here extended to apply to the coupling between variables, namely, the pressure and velocity, whereby operations involving different variables are split into a series of predictor-corrector steps.

    • 算子分裂的一个问题在于,它有时候不能保证到达稳态解的正确性!详见参考文献14

算法的具体描述请参考OpenFOAMwiki的PISO Algorithm (这里的描述所用的版本有点老,和最新版OpenFOAM中使用的变量名有所不同):

  1. Set the boundary conditions 施加边界条件.
  2. Solve the discretized momentum equation to compute an intermediate velocity field. 解离散动量方程得到中间阶段的速度场$U^*$
    1. 但是其实这一步可以设置momentumPredictor=false来跳过去。
    2. 因为最重要的其实是压力方程,这个速度场方程UEqn主要用来提取$A$和$H$,以构造后面压力方程要用的phi
    3. 因为有非定常项,所以这里是没有松弛的!有松弛就不对了。
  3. Compute the mass fluxes at the cells faces,构造phiHbyA.
    1. 与simple不同,这里的phiHbyA加入了fvc::ddtCorr
    2. fvc::ddtCorr项的定义参考这里
  4. Solve the pressure equation.解压力方程,其实应该是解的PLE
  5. Correct the mass fluxes at the cell faces.更新phi
  6. Correct the velocities on the basis of the new pressure field. 更新U
  7. Update the boundary conditions.更新U的边界条件
  8. Repeat from 3 for the prescribed number of times. 重复piso循环,
    1. 这个是由piso.correct()控制的,而piso.correct()是由变量label corrPISO_控制的。
    2. 重复的时候,UEqn所引用的U是有更新的。所以用UEqn构造的UEqn.A(), UEqn.H()也可能会有所更新。
  9. Increase the time step and repeat from 1. 下一个时间步,这里是由runTime.loop()控制的。

OpenFOAM中的PISO求解器其实有两个,最简单的是icoFoam,用的是PISO算法解层流。复杂一点的是pisoFoam,带有湍流模型。

其他

solutionControl及其子类的继承关系

  • pisoControl是pimpleControl的子类!
  • pimpleContrl是solutionControl的子类
  • simpleControl是solutionControl的子类。

松弛

从代码中可以发现,pisoFoam和icoFoam对于压力和动量方程是没有进行松弛的。这是因为它要进行的是时间精确的模拟。

但是如果你用的ddt格式本来就不是时间精确的话,额,还是不要用PISO了,改用PIMPLE吧。加点松弛因子可能还稳定些。

PIMPLE算法

PIMPLE算法是SIMPLE和PISO的结合。因为PISO只有一个预测步,所以PISO的稳定性限制了它只能用于瞬态计算并且时间步不能取得太大,而SIMPLE没有时间项只能用于稳态计算(或者采用伪时间步方法用于瞬态计算)。PIMPLE算法的核心是把UEqn重新生成并解了很多次!

  • pimpleFoam中的#include "UEqn.H" 语句是被包在pimple.loop()循环内的。
  • pisoFoam中的#include "UEqn.H"语句没有被包裹,每次runTime.loop()只运行一次。
  • pisoFoam中没有使用consistent关键字,所以没有SIMPLE-C修正。
  • pimpleFoam中使用了consistent关键字,所以和simpleFoam一样,对phiHbyA进行了SIMPLE-C处理。
  • simpleFoam中使用的simple.loop()而不是runTime.loop(),但是它调用了runTime.loop()

PIMPLE里的三层循环

总结起来,PIMPLE算法每个时间步(由runTime.loop()控制)有三层循环:

  • 外层pimple.loop()类似SIMPLE循环
    • 利用当前的U,生成并求解UEqn
    • 中层pimple.correct(),几乎就是PISO循环,负责压力求解
      • 利用UEqn 构造rAU, HbyA
      • HbyA 生成phiHbyA,加入fvc::ddtCorr(U,phi)修正。
      • adjustPhi()
      • 利用consistent关键字对phiHbyA进行SIMPLE-C一致性修正
      • constrainPressure()
      • 内层pimple.correctNonOrthogonal()非正交修正循环
        • 构造压力联系方程pEqn
        • pEqn.setReference()
        • 用新的p作为初始值求解pEqn,进行非正交修正。
      • 更新Uphi

而SIMPLE算法没有中层循环,PISO算法没有外层循环。所以PIMPLE算法可以简化为PISO算法,但是由于SIMPLE算法中没有包含非定常项,所以PIMPLE算法要等同于SIMPLE算法需要使ddt项的格式需要指定为steadyState

一些编程细节

此外在程序代码中还有好多稀奇古怪的乱七八糟玩意儿。

solutionControl及其子类的变量命名规范

OpenFOAM里的求解控制类是solutionControl 及其子类simpleControl, pisoControlpimpleControl。其中定义了很多变量,其命名规范是:

  • label xxx_表示迭代变量,比如corr_, corrNonOrtho_,pimpleControl中的 corrPISO_等。
  • label nXxx_表示迭代上限次数,比如nCorrNonOrtho_(实际迭代次数是nCorrNonOrtho_+1,因为第一次求解不算所谓修正),pimpleControl中的nCorrPIMPLE_, nCorrPISO_
  • bool yyy_表示标志量,比如transonic_, consistent_, momentumPredictor_
  • residualControl_ 残差控制。
  • 以pimpleControl为例,其迭代循环变量和控制变量标志如下表
含义 迭代量的变量名 迭代上限的变量名 配置文件中的关键字 默认值 判断是否执行的函数
外部迭代,又叫PIMPLE迭代 label corr_ nCorrPIMPLE_ nOuterCorrectors 1 loop()
PISO压力修正循环 label corrPISO_ nCorrPISO_ nCorrectors 1 correct()
非正交修正 label corrNonOrtho_ nCorrNonOrtho_ 这是实际迭代次数上限减1 nNonOrthogonalCorrectors 0 correctNonOrthogonal()
是否执行动量预测步 bool momentumPredictor_ momentumPredictor true momentumPredictor()
一致性 bool consistent_ consistent false consistent()
跨声速 bool transonic transonic false transonic()
是否在PISO循环外更新密度 bool SIMPLErho_ SIMPLErho false SIMPLErho()
仅在最后迭代启用湍流 bool turbOnFinalIterOnly_ turbOnFinalIterOnly true momentumPredictor()

simpleFoam中phi变量的更新。

如果仔细的人应该会发现,phi这个变量在simpleFoam中是由createFields.H调用createPhi.H创建的,然后在UEqn.H中被fvm::div(phi,U)使用,那么在下一次迭代之前它是如何更新的呢?

我们发现它是在pEqn.H中更新的。

1
2
3
4
5
//pEqn.H
if (simple.finalNonOrthogonalIter())
{
phi = phiHbyA - pEqn.flux();
}

关于adjustPhi, 和 setReference

adjustPhi和pRef针对的是压力方程为全Neuman边界的情况12,这时解的存在唯一性可能会有两个问题:

  • 唯一性:方程系数矩阵M 列秩不满时,解有无穷多(例如$x_1+x_2=1,2x_1+2x_2=2$),互之间相差一组常数(或者一组常函数,或零空间元素),所以需要setReference
  • 存在性:方程的增广矩阵秩和M的秩不一样时(例如$x_1+x_2=1,-x_1-x_2=-2$),存在相容性问题,需要adjustPhi

关于线性代数知识可以参考23

23. 教案: 线性方程组

方程多解问题:setReference

当方程为:

首先,如果$p_1$是一个解,那么$p_2=p_1+const$也是方程的解,方程解不唯一,所以离散方程会有奇异性,无法求解。
解决方法是固定一个点的压力为给定值,就是设置setReference

方程相容性问题:adjustPhi

p变量前面是直接加了微分的,所以需要对于全Neumann条件的情形还需要满足相容性条件

这就是adjustPhi去强制满足这个条件。
带入f和g的表达式

可得:

你再看看adjustPhi出现的位置,是在pEqn.H中,对phiHbyA进行的调整。其实也就是这里的$\mathbf u^*$项。必须使其满足以上相容性条件。

关于constrainPressureconstrainHbyA

从OpenFOAM的4.0版本开始,引入了fixedFluxPressure,fixedFluxExtraoplatedPressure边界条件,从而使入口的压力可以浮动,且边界压力梯度还可以根据计算情况进行调整,于是引入了constrainPressure()constrainHbyA()函数。

其说明可以参考Fumiya Nozaki’s CFD Blog的这篇blog:Recent changes in the basics of the OpenFOAM solversfixedFluxPressure 境界条件

constrainPressure用于调整压力p的边界

可压缩流中边界速度$U_f$在离散生成UEqn之后可以表示为:

所以界面法向压力梯度可以表示为:

将上述关系运用于压力边界条件,就是constrainPressure要干的事儿。

  • void constrainPressure(p,rho=1,U,phiHbyA,rhorAU,MRF=NULLMRF)
    • 输出调整的其实是p
    • 必要的输入有U,phiHbyA,rhorAU
  • p边界条件继承于fixedFluxPressureFvPatchScalarField的才会调整。

constrainHbyA(tHbyA,U,p)用于返回调整后的tHbyANew变量

  • 循环每个边界
  • U的在这个边界的边界条件是assignable()==false
  • p的在这个边界是继承自fixedFluxExtrapolatedPressureFvPatchScalarField
  • tHbyANew的这个 这边界设置为和U的这个边界相同的值。
  • 返回tHbyANew
  • 为什么要这么做:我也不知道。可能这样做了之后constrainPressure()得到的压力和HbyA在固定速度的边界与前述PLE方程就是一致的了,毕竟前面的constrainPressure()是利用了HbyA生成的phiHbyA的。
  • 如果没有用到fixedFluxExtrapolatedPressurefixedFluxPressure边界条件,这两个函数其实相当于没有作用。

关于fvc::ddtCorr()和fvc::ddtPhiCorr()

对于非稳态的求解器,如pimpleFoam, pisoFoam和icoFoam中,phiHbyA在构造时除了用到HbyA 向的插值,还会用到fvc::ddtCorr(),这一项主要是为了修正非定常和松弛因子对解的影响。而且这一项还是和OpenFOAM版本历史相关的!

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
//OpenFOAM 5.x, 4.x
//applications/solvers/incompressible/icoFoam/icoFoam.C +77
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::flux(HbyA) //这里用了fvc::flux函数,而不是原来的interpolate 和 & 算符
+ fvc::interpolate(rAU)*fvc::ddtCorr(U, phi) //用的是ddtCorr!
);

//OpenFOAM 2.3.x, 2.4.x
//applications/solvers/incompressible/icoFoam/icoFoam.C +73
surfaceScalarField phiHbyA
(
"phiHbyA",
(fvc::interpolate(HbyA) & mesh.Sf())
+ fvc::interpolate(rAU)*fvc::ddtCorr(U, phi) //用的是ddtCorr
);
//OpenFOAM 2.2.x
//applications/solvers/incompressible/icoFoam/icoFoam.C +73
surfaceScalarField phiHbyA//用了新的变量名!
(
"phiHbyA",
(fvc::interpolate(HbyA) & mesh.Sf())
+ fvc::ddtPhiCorr(rAU, U, phi) //用的是ddtPhiCorr!
);


//OpenFOAM 2.1.x
//applications/solvers/incompressible/icoFoam/icoFoam.C +72
phi = (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rAU, U, phi);

//foam-extend 4.0
//applications/solvers/incompressible/icoFoam/icoFoam.C +86
U = HUEqn.H()/aU; //注意这里的HUEqn不包含非定常项
phi = (fvc::interpolate(U) & mesh.Sf()); //所以这里也不包含。

太老的就不管了,因为根据OpenFOAM.org官网的更新说明,2.3.0版本时已经把ddtPhiCorr换成了ddtCorr

首先找OpenFOAM 5.x的代码,看看ddtCorr()到底是干啥的 :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
//src/finiteVolume/finiteVolume/fvc/fvcDdt.C +183
template<class Type>
tmp<GeometricField<typename flux<Type>::type, fvsPatchField, surfaceMesh>>
ddtCorr
(
const GeometricField<Type, fvPatchField, volMesh>& U,
const GeometricField
<
typename flux<Type>::type,
fvsPatchField,
surfaceMesh
>& phi
)
{
return fv::ddtScheme<Type>::New //根据变量名返回tmp<ddtScheme<Type>>的临时对象
(
U.mesh(),
U.mesh().ddtScheme("ddt(" + U.name() + ')')
).ref().fvcDdtPhiCorr(U, phi);
}

//这里涉及到RTS机制,
//为了简化,参考steadyState和经典的Euler格式是怎么做的

//src/finiteVolume/finiteVolume/ddtSchemes/steadyStateDdtScheme/steadyStateDdtScheme.C +305
//steadyState返回的是0,无修正,这很正常。
//返回值的量纲是:phi.dimensions()/dimTime
//对于不可压缩流应该是L^3/T^2
template<class Type>
tmp<typename steadyStateDdtScheme<Type>::fluxFieldType>
steadyStateDdtScheme<Type>::fvcDdtPhiCorr
(
const GeometricField<Type, fvPatchField, volMesh>& U,
const fluxFieldType& phi
)
{
return tmp<fluxFieldType>
(
new fluxFieldType
(
IOobject
(
"ddtCorr(" + U.name() + ',' + phi.name() + ')',
mesh().time().timeName(),
mesh()
),
mesh(),
dimensioned<typename flux<Type>::type>
(
"0",
phi.dimensions()/dimTime,
Zero
)
)
);
}

//Euler格式的修正
//src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C +534
template<class Type>
tmp<typename EulerDdtScheme<Type>::fluxFieldType>
EulerDdtScheme<Type>::fvcDdtPhiCorr
(
const GeometricField<Type, fvPatchField, volMesh>& U,
const fluxFieldType& phi
)
{
dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT(); //时间增量的倒数

fluxFieldType phiCorr
(
phi.oldTime() - fvc::dotInterpolate(mesh().Sf(), U.oldTime())
);

return tmp<fluxFieldType>
(
new fluxFieldType
(
IOobject
(
"ddtCorr(" + U.name() + ',' + phi.name() + ')',
mesh().time().timeName(),
mesh()
),
this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr)
*rDeltaT*phiCorr
)
);
}

//涉及到fvcDdtPhiCoeff系数,继续查找,发现这个函数位于ddtScheme中。
//src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.C +151
template<class Type>
tmp<surfaceScalarField> ddtScheme<Type>::fvcDdtPhiCoeff
(
const GeometricField<Type, fvPatchField, volMesh>& U,
const fluxFieldType& phi,
const fluxFieldType& phiCorr
)
{
tmp<surfaceScalarField> tddtCouplingCoeff = scalar(1)
- min
(
mag(phiCorr)
/(mag(phi) + dimensionedScalar("small", phi.dimensions(), SMALL)),
scalar(1)
);

surfaceScalarField& ddtCouplingCoeff = tddtCouplingCoeff.ref();

surfaceScalarField::Boundary& ccbf =
ddtCouplingCoeff.boundaryFieldRef();

forAll(U.boundaryField(), patchi)
{
if
(
U.boundaryField()[patchi].fixesValue()
|| isA<cyclicAMIFvPatch>(mesh().boundary()[patchi])
)
{
ccbf[patchi] = 0.0;
}
}

if (debug > 1)
{
InfoInFunction
<< "ddtCouplingCoeff mean max min = "
<< gAverage(ddtCouplingCoeff.primitiveField())
<< " " << gMax(ddtCouplingCoeff.primitiveField())
<< " " << gMin(ddtCouplingCoeff.primitiveField())
<< endl;
}

return tddtCouplingCoeff;
}

所以OpenFOAM的ddtCorr的算法如下:

  1. 首先用过去时间的通量$\phi^{n-1}$和速度$U^{n-1}$估算通量的修正$\phi^{c}$。
  1. 计算ddt耦合系数$K_c$:

    同时,应该注意到,ddt耦合系数$K_c$是一个大于0的无量纲界面张量场!它的边界条件也被特殊地处理了,对于固定速度边界或cyclicAMI边界,它的值是0.

  2. 输出修正量$\phi^k = K_c\phi^c/\Delta t$

但是为什么这么计算耦合系数$K_c$和通量$\phi^c$,我也不知道!

  • cfd-online上一个叫eugene的id提到ddtCorr实质上是文献19中的所谓Choi Correction,对应的原始文献是20 ,但发帖者认为fvcDdtPhiCoeff = 1时是完整的Choi Correction,但此时OpenFOAM会不稳定,所以用了一个经验性的fvcDdtPhiCoeff()函数来限制这个Choi Correction。

Jasak的搞法

Jasak似乎玩了另外的花样。在foam-extend 4.0中,Jasak的icoFoam和OpenFOAM.org以及OpenFOAMplus的icoFoam就不一样了。

  • Jasak的icoFoam把UEqn分成了HUEqnddtUEqn两部分:

    • HUEqn: 对流和扩散项
    • ddtUEqn:非定常项;
  • 动量预测步:HUEqn + ddtUEqn == -grad(p)

  • PISO循环中:HHUEqn.H()中产生!与时间步完全无关。

  • 速度更新:

  • 1
    2
    3
    4
    5
    6
    7
    //https://github.com/Unofficial-Extend-Project-Mirror/foam-extend-foam-extend-4.0/blob/master/applications/solvers/incompressible/icoFoam/icoFoam.C +113
    // Note: cannot call H(U) here because the velocity is not complete
    // HJ, 22/Jan/2016
    U = 1.0/(aU + ddtUEqn.A())*
    (
    U*aU - fvc::grad(p) + ddtUEqn.H()
    );

Jasak 认为这使得代码更加清晰了(时间导数算子中完全去掉了ddtPhiCorrddtCorr)。而且完全消除了时间步和松弛因子对结果的影响。结果对比参考文献17

投影法

另外一种选择是用投影法,参考文献21。 投影法在同位网格上可以避免Rhie-Chow插值,从而避免这个ddtCorr()的问题。 他们也做了数值试验来观察这个质量通量修正项(mass flux correction term)对PISO算法的效果。发现扔掉这一项,PISO算法的耗散大大降低。

OpenFOAMplus的搞法

OpenFOAMplus在v1706中的ddtScheme加入了一个scalar ddtPhiCoeff_的定义。从而使ddtCorr()的实现更加复杂化。

代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
//https://develop.openfoam.com/Development/OpenFOAM-plus/blob/master/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.H
template<class Type>
class ddtScheme
:
public tmp<ddtScheme<Type>>::refCount
{

protected:

// Protected data

const fvMesh& mesh_;

//- Input for fvcDdtPhiCoeff (-1 default)
scalar ddtPhiCoeff_; //ebf654f2这个集成rhoPimpleAdiabaticFoam的提交中才出现的。

//...

而此时更改为了:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
// https://develop.openfoam.com/Development/OpenFOAM-plus/blob/master/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.C +151
template<class Type>
tmp<surfaceScalarField> ddtScheme<Type>::fvcDdtPhiCoeff
(
const GeometricField<Type, fvPatchField, volMesh>& U,
const fluxFieldType& phi,
const fluxFieldType& phiCorr
)
{
tmp<surfaceScalarField> tddtCouplingCoeff
(
new surfaceScalarField
(
IOobject
(
"ddtCouplingCoeff",
U.mesh().time().timeName(),
U.mesh()
),
U.mesh(),
dimensionedScalar("one", dimless, 1.0)
)
);

surfaceScalarField& ddtCouplingCoeff = tddtCouplingCoeff.ref();

if (ddtPhiCoeff_ < 0)
{
ddtCouplingCoeff -= min
(
mag(phiCorr)
/(mag(phi) + dimensionedScalar("small", phi.dimensions(), SMALL)),
scalar(1)
);
}
else
{
ddtCouplingCoeff =
dimensionedScalar("ddtPhiCoeff", dimless, ddtPhiCoeff_);
}

surfaceScalarField::Boundary& ccbf = ddtCouplingCoeff.boundaryFieldRef();

forAll(U.boundaryField(), patchi)
{
if
(
U.boundaryField()[patchi].fixesValue()
|| isA<cyclicAMIFvPatch>(mesh().boundary()[patchi])
)
{
ccbf[patchi] = 0.0;
}
}

if (debug > 1)
{
InfoInFunction
<< "ddtCouplingCoeff mean max min = "
<< gAverage(ddtCouplingCoeff.primitiveField())
<< " " << gMax(ddtCouplingCoeff.primitiveField())
<< " " << gMin(ddtCouplingCoeff.primitiveField())
<< endl;
}

return tddtCouplingCoeff;
}

这样使得可以这个系数可以手工调整。

他们所引用的参考文献22 暂时没找到全文

关于p,pFinal的选择

仔细的哥们儿会发现simpleFoam里用的是pEqn.solve(),到了icoFoam, pisoFoam和pimpleFoam里用的就是pEqn.solve(mesh.solver(p.select(piso/pimple.finalInnerIter())))

这里的区别,还是看代码吧:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
//src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H +92
//最后一步内迭代,也就是最后一次用PISO求解压力且是非正交修正的最后一步时返回true,否则返回false
inline bool Foam::pimpleControl::finalInnerIter() const
{
return
corrPISO_ == nCorrPISO_
&& corrNonOrtho_ == nNonOrthCorr_ + 1;
}

//src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C +985
//true,也就是最后一次迭代时返回名字`pFinal`,否则返回名字`p`
template<class Type, template<class> class PatchField, class GeoMesh>
Foam::word Foam::GeometricField<Type, PatchField, GeoMesh>::select
(
bool final
) const
{
if (final)
{
return this->name() + "Final";
}
else
{
return this->name();
}
}

// src/OpenFOAM/matrices/solution/solution.C +353
// 返回name对应的solver controls dictionary
const Foam::dictionary& Foam::solution::solver(const word& name) const
{
if (debug)
{
Info<< "Lookup solver for " << name << endl;
}

return solvers_.subDict(name);
}

//src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C +55
//按这个solverControls来解
template<class Type>
Foam::SolverPerformance<Type> Foam::fvMatrix<Type>::solve
(
const dictionary& solverControls
)

结论:simpleFoam是稳态求解器,所以没有pFinal一说,所以用不着这套复杂的机制。

关于收敛判据和求解控制

对于不可压缩流动而言, 收敛判据只取决于两个问题,

  1. 连续性方程是否得到满足,在多大程度上得到满足?
    1. OpenFOAM里是用#include "continuityErrs.H"来计算、记录和输出这一项的。
      1. 程序开始时还需要用 include "initContinuityErrs.H" 初始化一下。
      2. 它用的是phi而不是U
      3. 它会记录累计误差cumulativeContErr,全场流量平衡误差globalContErr,以及局部误差的和sumLocalContErr(最后一项几乎等价于速度场散度的L1模:$|\nabla\cdot U|_{L1}$)
    2. 事实上这点比较简单,因为每次计算压力的时候都是力图使连续性方程得到满足的。
    3. 而且这一项与非定常项无关。
  2. 动量方程在多大程度上得到满足:
    1. 这点在OpenFOAM里是通过对动量方程的求解来满足的。
    2. 其实OpenFOAM并没有对其进行严格地度量(指L1, L2模这样的数学度量)。
    3. 对于稳态算法SIMPLE,看最终一步的动量方程的的初始残差
    4. 对于非稳态算法(PISO, PIMPLE等)输出中的每个时间步最后一次外迭代时,动量方程的初始残差可以认为是一个度量,但是它并不完全一致和严格。
      1. 每个时间步 是因为非稳态算法一般是需要进行时间精确模拟的。如果这个时间步没有收敛,这个时间步就不对,并且会影响到下一个时间步。
      2. 最后一次 是因为前面的算是中间过程,最后并不保存下来;
      3. 外迭代是因为内迭代是解耦计算,有时候看上去残差很小,可能并没有卵用,一耦合残差又超级大。
        1. 对SIMPLE对应着simple.loop()
        2. 对PIMPLE对应着pimple.loop()
      4. 初始残差是因为OpenFOAM输出的残差都是线性求解器输出的,而其中线性求解器输出的只有初始残差是和非线性残差是一致的。
        1. 经过线性求解器的迭代,线性残差可能会很小,甚至到机器零;
        2. 但是这个解只有对初始解产生的矩阵系数是对的,如果把矩阵系数按照得到的解重新计算,得到的残差才算是非线性残差(严格说来也不完全一样,但是差不多啦);
        3. 这个非线性残差小才算是收敛 ;
        4. 另外,即使非线性残差小不代表误差小;
        5. 时间步很大时,网格很粗时,即使残差为0,误差也可能很大;
      5. 实际得到的解因为是最后一次外迭代完成之后的解,一般来说其非线性残差可能会比输出中的初始残差更小一些。但是这个数据并没有被计算和输出。
    5. 参考OpenFOAM的残差定义

参考资料

1. MATHEMATICS, NUMERICS, DERIVATIONS AND OPENFOAM®
2. Patankar, S. V. and Spalding, D.B. (1972), “A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows”, Int. J. of Heat and Mass Transfer, Volume 15, Issue 10, October 1972, Pages 1787-1806
3. Pressure Correction Scheme for Incompressible Fluid Flow
4. S.V. Patankar. Numerical Heat Transfer and Fluid Flow. Hemisphere, 1980
5. Caretto L.S., Gosman A.D., Patankar S.V., Spalding D.B. (1973) Two calculation procedures for steady, three-dimensional flows with recirculation. In: Cabannes H., Temam R. (eds) Proceedings of the Third International Conference on Numerical Methods in Fluid Mechanics. Lecture Notes in Physics, vol 19. Springer, Berlin, Heidelberg
6. http://thevisualroom.com/poisson_for_pressure.html
7. Lee S L, Tzong R Y. Artificial pressure for pressure-linked equation[J]. International journal of heat and mass transfer, 1992, 35(10): 2705-2716.
8. http://dyfluid.com/RhieChow.html
9. Nonlinear Systems: Picard and Newton methods
10. 林建忠,《流体力学》
11. DISCUSSION ON MOMENTUM INTERPOLATION METHOD FOR COLLOCATED GRIDS OF INCOMPRESSIBLE FLOW. Bo Yu, Wen-Quan Tao, Jin-Jia Wei, Yasuo Kawaguchi, Toshio Tagawa & Hiroyuki Ozoe. Numerical Heat Transfer, Part B: Fundamentals Vol. 42 , Iss. 2,2002
12. CS205b/CME306 Lecture 16
13. Collocated grids
14. Operator Splitting
15. Issa, R. I. (1986). Solution of the implicitly discretised fluid flow equations by operator-splitting. Journal of computational physics, 62(1), 40-65.
16. Van Doormaal, J. P., & Raithby, G. D. (1984). Enhancements of the SIMPLE method for predicting incompressible fluid flows. Numerical heat transfer, 7(2), 147-163.
17. Numerics Improvements and Validation Results: FOAM-Extend Update on Work-in-Progress
18. Moukalled, F., Mangani, L., & Darwish, M. (2016). The finite volume method in computational fluid dynamics.
19. Yu, B., Tao, W. Q., Wei, J. J., Kawaguchi, Y., Tagawa, T., & Ozoe, H. (2002). Discussion on momentum interpolation method for collocated grids of incompressible flow. Numerical Heat Transfer: Part B: Fundamentals, 42(2), 141-166.
20. S. K. Choi, Note on the Use of Momentum Interpolation Method for Unsteady Flows, Numer. Heat Transfer A, vol. 36, pp. 545-550, 1999.
21. Vuorinen, V., Keskinen, J. P., Duwig, C., & Boersma, B. J. (2014). On the implementation of low-dissipative Runge–Kutta projection methods for time dependent flows using OpenFOAM®. Computers & Fluids, 93, 153-163.
22. Knacke, T. (2013). Potential effects of Rhie & Chow type interpolations in airframe noise simulations. In: Schram, C., Dénos, R., Lecomte E. (ed): Accurate and efficient aeroacoustic prediction approaches for airframe noise, VKI LS 2013-03.
如果觉得此文不错,欢迎大爷/姐姐打赏