【原创教程】自带 solver :基础求解器(三)icoFoam

我们从最简单的 icoFoam 开始，在安装完 OpenFOAM 后，我们常用 icoFoam 来检验一下有没有安装好。路径：

OpenFOAM/applications/solvers/Incompressible/icoFoam

包含有三个代码文件和两个 Make 文件夹中的文件：

• icoFoam.C——求解器顶层源代码

• createFields.H——头文件

• Make/files——按行存储所有源代码文件名，最后一行指定目标代码 EXE 的名称和存放位置

• Make/options——设定查找头文件和库的路径 EXE_INC 和需要链接的库 EXE_LIBS

3. icoFoam

3.1 概览

老样子，先来看看 OpenFOAM 是如何描述这个求解器的：

Application
    icoFoam
Group
    grpIncompressibleSolvers
Description
    Transient solver for incompressible, laminar flow of Newtonian fluids.

向右滑动显示完整内容

一个用于不可压缩的层流牛顿流体的求解器。一般来说不可压缩牛顿流体的控制方程可以写为：

图2 计算节点

显式方法的优点是方程建立及编程相对简单。

然而也有无法避免的缺点，每一个差分方程中的时间增量和位移增量是相互影响的，即一旦位移增量取定，时间增量就不是独立变量了，是要受到稳定性条件的限制。

即其取值必须不能大于稳定性条件的限制，否则，时间推进的过程将变得非常不稳定，计算程序也会因为数字趋于无穷大或对负数开平方的一系列原因而终止运行。

（2）隐式差分
隐式格式的方程并非彼此独立，而是在网格节点上形成相互耦合的方程组，相应的未知数在方程组中同时求解出来。

\heading Solver details
The solver uses the PISO algorithm to solve the continuity equation:
        \f[
            \div \vec{U} = 0  //PISO算法
        \f]
and momentum equation:   //动量（预测）方程
        \f[
            \ddt{\vec{U}}
          + \div \left( \vec{U} \vec{U} \right)
          - \div \left(\nu \grad \vec{U} \right)
          = - \grad p
        \f]
Where:
    \vartable  //变量表
        \vec{U}  | Velocity
        P       | Pressure
    \endvartable

    \heading Required fields //初始化
    \plaintable
        U       | Velocity [m/s]
        p       | Kinematic pressure, p/rho [m2/s2]
    \endplaintable

这是对 icoFoam 控制方程的一个解释，更好地理解代码含义，在将来的介绍中，应该不需要再刻意提醒大家了。

int main(int argc, char *argv[])
{
    argList::addNote
    (
        "Transient solver for incompressible, laminar flow"
        " of Newtonian fluids."
    );

求解器的 main 方程开始，并作出提醒：一个应用于瞬态层流不可压缩牛顿流体的求解器。

    Info<< "\nStarting time loop\n" << endl;
    while (runTime.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

开始时间迭代，并且通过 Info 语句输出模拟时间。

        fvVectorMatrix UEqn
        (
           fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
        );
        if (piso.momentumPredictor())
        {
            solve(UEqn == -fvc::grad(p));
        }

动量方程的代码：

• fvm::ddt(U)
  对方程
  

       中的时间项进行对速度 U 关于时间 t 的欧拉全隐离散

• fvm::div(phi, U)

  对流项的隐性离散

• fvm::laplacian(nu, U)
  拉普拉斯项的隐性离散

• fvc::grad(p)
  压力项的显性离散
  

// --- PISO loop
        while (piso.correct())
        {
            volScalarField rAU(1.0/UEqn.A());
            volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
            surfaceScalarField phiHbyA
            (
                "phiHbyA",
                fvc::flux(HbyA)
              + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
            );
adjustPhi(phiHbyA, U, p);

开始 PISO 算法的 while 循环，声明了一个变量 phiHbyA ，很明显是一个叫做 HbyA 的变量的通量，其中变量 HbyA 方程表达式可以写为：

while (piso.correctNonOrthogonal())
            {
                // Pressure corrector
                fvScalarMatrix pEqn
                (
                    fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
                );

                pEqn.setReference(pRefCell, pRefValue);
                pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));

                if (piso.finalNonOrthogonalIter())
                {
phi = phiHbyA - pEqn.flux(); 
}
}

压力预测方程，或者说是压力控制方程，再或者说叫做压力泊松方程 :

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

即：

U = HbyA - rAU*fvc::grad(p);
U.correctBoundaryConditions();
}
runTime.write(); 
runTime.printExecutionTime(Info);
}

更新速度，作为下一次迭代的初始值来进行下一次计算。

整体来说，icoFoam的一次迭代过程包含有：

• 依据初始条件求解预测速度；

• 通过速度组建HbyA；

• 求解压力方程；

• 根据求解的压力以及HbyA更新速度；

• 回到第二步进行下一次迭代。

参考：
[1]http://openfoamwiki.net/
[2]http://www.thevisualroom.com/explicit_implicit.html
[3]http://dyfluid.com/icofoam.html