OpenFoam中的Simple系列算法_openfoam 中simple.consistent

作者：笔触狂放9 | 2024-05-03 11:15:31

踩

openfoam 中simple.consistent

简介

Simple系列算法是求解不可压缩流的主要工具。OpenFoam中提供了SimpleFoam求解器、PisoFoam求解器以完成不可压缩流的求解。

但是OpenFoam中的代码与流体力学教材中的讲述并不一致，导致初学者难以理解代码的含义。笔者推荐一篇Martínez撰写的文章《Influence of momentum interpolation methods on the accuracy and convergence of pressure–velocity coupling algorithms in OpenFOAM》，该文章详细地推导了Simple算法、SimpleC算法和Piso算法，更重要的是其最终形式就是OpenFoam所采用的形式。

本文将直接使用Martínez的推导结果，将结果与OpenFoam中的代码对应，方便初学者理解代码含义。

Simple算法

文中公式（40）表示使用Rhie-Chow插值得到的面心速度
$u_{e}^{m}=\frac{\alpha_{u}\left(\sum_{N B} A_{N B} u_{N B}^{*}+B_{P}\right)_{e}}{\left(A_{P}\right)_{e}}-\frac{\alpha_{u} \Delta y\left(p_{E}^{m}-p_{P}^{m}\right)}{\left(A_{P}\right)_{e}}$ 面心速度是用来构造连续方程的。考虑不可压缩流，对于任意一个单元，进出各面心的流量之和必然为0。所以针对一个单元，对上式求和（面积加权），即可得到压力方程
$\sum {{S_f}{{{\alpha _u}\Delta y\left( {p_E^m - p_P^m} \right)} \over {{{\left( {{A_P}} \right)}_f}}}} = \sum {{S_f}{{{\alpha _u}{{\left( {\sum\limits_{NB} {{A_{NB}}} u_{NB}^* + {B_P}} \right)}_f}} \over {{{\left( {{A_P}} \right)}_f}}}}$
现在看看SimpleFoam求解器中的代码

fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA)
1

代码与公式是严格对应的，rAtU是动量方程对角矩阵Ap的倒数。phiHbyA是H矩阵与Ap之商在面心插值后的结果。总之，代码与公式是严格对应的。

SimpleC算法

文中公式（48）是使用SimpleC算法和Rhie-Chow插值得到的面心速度
$u_{e}^{m}=\frac{\alpha_{u}\left(\sum_{N B} A_{N B} u_{N B}^{*}+B_{P}\right)_{e}}{\left(A_{P}\right)_{e}}-\alpha_{u}\left(\frac{1}{\left(A_{P}\right)_{e}}-\frac{1}{\left(A_{P}-\sum_{N B} A_{N B}\right)_{e}}\right) \Delta y\left(p_{E}^{m-1}-p_{P}^{m-1}\right)-\frac{\alpha_{u} \Delta y\left(p_{E}^{m}-p_{P}^{m}\right)}{\left(A_{P}-\sum_{N B} A_{N B}\right)_{e}}$
同样对面求和之后有
$\sum {{S_f}{{{\alpha _u}\Delta y\left( {p_F^m - p_P^m} \right)} \over {{{\left( {{A_P} - \sum\limits_{NB} {{A_{NB}}} } \right)}_f}}}} = \sum {{S_f}\left( {{{{\alpha _u}{{\left( {\sum\limits_{NB} {{A_{NB}}} u_{NB}^* + {B_P}} \right)}_f}} \over {{{\left( {{A_P}} \right)}_f}}} - {\alpha _u}\left( {{1 \over {{{\left( {{A_P}} \right)}_f}}} - {1 \over {{{\left( {{A_P} - \sum\limits_{NB} {{A_{NB}}} } \right)}_f}}}} \right)\Delta y\left( {p_F^{m - 1} - p_P^{m - 1}} \right)} \right)}$
现在看看SimpleFoam求解器中的代码
SimpleC在Simple算法的基础上对一些系数做了修改，如下

if (simple.consistent())
{
    rAtU = 1.0/(1.0/rAU - UEqn.H1());
    phiHbyA +=
        fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf();
    HbyA -= (rAU - rAtU())*fvc::grad(p);
}
1
2
3
4
5
6
7

rAtU 变成了 $\over {{{\left( {{A_P} - \sum\limits_{NB} {{A_{NB}}} } \right)}_f}}}$ ，phiHbyA变成了更复杂的形式。
压力泊松方程依然没有变化

fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA)
1

代码与公式是严格对应的。

Piso算法

文中公式（55）是使用Piso算法和Rhie-Chow插值得到的面心速度
$u_{e}^{m}=\frac{\alpha_{u}\left(\sum_{N B} A_{N B} u_{N B}^{* *}+B_{P}\right)_{e}}{\left(A_{P}\right)_{e}}-\frac{\alpha_{u} \Delta y\left(p_{E}^{m}-p_{P}^{m}\right)}{\left(A_{P}\right)_{e}}$
代码与公式之间是严格对应的。

事实上，Piso算法与Simple算法的压力修正方程并无区别。在OpenFoam中，Piso算法是为了算瞬态，而Simple算法只能算稳态。二者之间的区别主要在压力修正的计算次序上

Piso算法的流程如下：

while (runTime.loop())
     {
         Info<< "Time = " << runTime.userTimeName() << nl << endl;
 
         #include "CourantNo.H"
 
         // Pressure-velocity PISO corrector
         {
             fvModels.correct();
 
             #include "UEqn.H"
 
             // --- PISO loop
             while (piso.correct())
             {
                 #include "pEqn.H"
             }
         }
 
         viscosity->correct();
         turbulence->correct();
 
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
             << nl << endl;
     }
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

代码中的runTime是物理时间，在每个物理时间步下完成速度预测与若干次压力修正。收敛之后，迈向下一个时间步。

Simple算法的流程如下：

while (simple.loop(runTime))
     {
         Info<< "Time = " << runTime.userTimeName() << nl << endl;
 
         fvModels.correct();
 
         // --- Pressure-velocity SIMPLE corrector
         {
             #include "UEqn.H"
             #include "pEqn.H"
         }
 
         viscosity->correct();
         turbulence->correct();
 
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
             << nl << endl;
     }
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

代码中的runTime不是物理时间，而是一个递增的整数，代表迭代次数。每一步迭代都需要求解速度预测方程与压力修正方程。收敛后得到的是稳态解。

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/笔触狂放9/article/detail/529352