icoFoam 解析

Submitted by rheoworks on Thu, 06/01/2017 - 21:45

icoFoam 是 OpenFOAM User Guide 里面的算例用的一个计算程序。它的源程序位于 $FOAM_SOLVERS/incompressible/icoFoam, 进入到这个目录里，用 tree 会输出以下内容（注：#及其后属于注释，并不是使用该命令的输出结果）：

.                        # 当前目录
|-- createFields.H       # 创建变量（场）
|-- icoFoam.C            # 主程序
`-- Make                 # 编译设置的文件夹
    |-- files            # 可执行文件的输出
    `-- options          # 编译选项、库添加

icoFoam 用于求解不可压缩层流 N-S 方程。其中，连续性方程为 $$ \nabla\mathbf{U}=0\tag{1} $$ 动量方程为 $$ \frac{\partial \mathbf{U}}{\partial t}+\nabla \cdot (\mathbf{U} \otimes\mathbf{U})=-\nabla \frac{p}{\rho}+\nabla \cdot(\nu \nabla \mathbf{U}) \tag{2} $$ 式中，$\mathbf U$ 为速度矢量，$p$ 为压力，$\rho $ 为密度。由于所求解的流场为不可压缩流，因此，在 icoFoam 程序中，$p/\rho$ 就由变量（或场） p 表示。

对动量方程进行全隐半离散化，方程左边第一项 $$ \iiint_V \int \frac{\partial \mathbf U}{\partial t}\,\mathrm{d}V\mathrm{d}t = \left( \mathbf U_\mathrm{P}^r - \mathbf U_\mathrm{P}^n \right)\, V_\mathrm{P} \tag{3} $$

方程左边第二项 $$ \begin{equation} \begin{split} \iiint_V \int \nabla \cdot \left( \mathbf{U}\otimes\mathbf{U} \right)\mathrm{d}V\,\mathrm{d}t = &\oint_S \int \mathbf{U}\otimes\mathbf{U}\,\mathrm{d}\mathbf{S}\,\mathrm{d}t \ = &\sum {\left( \mathbf{U}^n\otimes\mathbf{U}^r \right)}_f S_f\Delta t\ = &\sum \left( {F_f^n \mathbf{U}_f^r} \right) \Delta t \end{split} \end{equation} \tag{4} $$

方程右边第二项 $$ \begin{equation} \begin{split} \iiint_V \int \nabla \cdot \left( \nu \nabla \mathbf U \right)\mathrm{d}V\,\mathrm{d}t = &\oint_S\int \nu \nabla \mathbf U\,\mathrm{d} \mathbf{S}\,\mathrm{d}t\ = &\sum {\left( \nu \nabla \mathbf U^r \right)}_f S_f\Delta t\ = &\sum \left( \nu \left( \nabla_f^\bot \mathbf U^r \right) \right) \Delta t \end{split} \end{equation} \tag{5} $$

式中，上标 $n$ 表示为当前的时间步（已知），上标 $r$ 表示预测时间步（待求），下标 $f$ 表示网格单元面上的值，$V$ 表示网格单元体积，$\mathrm P$ 表示当前网格，$S$ 表示网格单元的各个面的面矢量，$\Delta t$ 表示时间步长，$F_f$ 为面通量。$\nu$ 为动力粘度。$\nabla_f^\bot{\mathbf{U}^r}$ 表示面的法向梯度。

$$ \begin{equation} \nabla_f^\bot \mathbf U^r = {( \nabla \mathbf U^r )}_f S_f = |S_f| X \end{equation} \tag{6} $$

其中 $d$ 表示网格单元体心之间的距离。下标 $N$ 表示相邻网格，下标 $P$ 表示当前网格。$\nabla_f^\bot \mathbf U^r$ 在三维情况下是一矢量, $X=\frac{\mathbf U_\mathrm{N}^n - \mathbf U_\mathrm{P}^r} {|d|}$。将上述三个离散项代入原方程，有

$$ \frac{\mathbf U_\mathrm{P}^r - \mathbf U_\mathrm{P}^n}{\Delta t} V_\mathrm{P} + \sum \left( F_f^n \mathbf U_f^r \right) = \sum \left( \nu \nabla_f^\bot \mathbf{U}^r \right) \tag{7} $$

在对流项的离散中，通量 $F_f$ 采用当前已知时间步的速度 $\mathbf U^n$ 来计算，并保留另一个预测速度 $\mathbf U_\mathrm{P}^r$ 为未知量。这种把其中一个速度设为已知，并保留另外一个未知速度的过程就是线性化：使用当前时间步的速度对高阶非线性对流项进行离散。由于 N-S 方程是非线性的，要么使用非线性求解器，要么将对流项线性化。由于非线性求解器非常复杂，因此在 OpenFOAM 中均采用线性化处理。线性化带来的一个问题就是通量的信息有所滞后。

动量预测

由于 $\mathbf U_f^r$ 为面上的预测速度，然而，在计算中求得的均为体心的速度。因此，需要从体心速度进行插值得来面心速度，此处可以引入各种格式。假设我们使用中心差分：

$$ \mathbf U_f^r = \frac{\mathbf U_\mathrm{P}^r + \mathbf U_\mathrm{N}^n}{2} \tag{8} $$

将 $(8)$ 代入离散方程 $(7)$，有

$$ \begin{equation} \frac{\mathbf U_\mathrm{P}^r - \mathbf U_\mathrm{P}^n}{\Delta t} V_\mathrm{P} + \sum \left( F_f^n \frac{ \mathbf U_\mathrm{P}^r + \mathbf U_\mathrm{N}^n}{2} \right) = \sum \left( \nu \left| S_f \right|\frac{\mathbf U_\mathrm{N}^n - \mathbf U_\mathrm{P}^r}{\left| d \right|} \right) \end{equation} \tag{9} $$

整理后，可得

$$ \begin{equation} \begin{split} &\left( \frac{V_\mathrm{P}}{\Delta t} + \sum \frac{F_f^n}{2} + \sum \nu \frac{\left| S_f \right|}{\left| d \right|} \right) \mathbf U_\mathrm{P}^r \ =& - \sum \left( \left( \frac{F_f^n}{2} - \nu \frac{\left| S_f \right|}{\left| d \right|} \right)\mathbf U_\mathrm{N}^n \right) + \frac{V_p}{\Delta t}\mathbf U_\mathrm{P}^n \end{split} \end{equation} \tag{10} $$

上式可简化为

$$ \begin{equation} A_\mathrm{P} \mathbf U_\mathrm{P}^r + \sum A_\mathrm{N} \mathbf U_\mathrm{N}^n - E_\mathrm{P}^n=0 \end{equation} \tag{11} $$

式中，

$$ \begin{equation} A_\mathrm{P}=\left( \frac{V_p}{\Delta t} + \sum \frac{F_f^n}{2} + \sum \nu \frac{\left| \mathbf S_f \right|}{\left| d \right|} \right) \end{equation} $$

$$ \begin{equation} A_\mathrm{N}=\left( {\frac{{F_f^n}}{2} - \nu \frac{{\left| \mathbf S_f \right|}}{{\left| d \right|}}} \right) \end{equation} $$

$$ \begin{equation} E_\mathrm{P}=\frac{V_p}{\Delta t}\mathbf U_\mathrm{P}^n \end{equation} $$

在公式 $(11)$ 中，就可加入上个时间步已知的压力梯度项。由于是不可压缩流，因此，在 icoFoam 中，使用 p 来代替 $p/\rho$，后面的不可压缩流中的 $p$ 均指 $p/\rho$

$$ \begin{equation} A_\mathrm{P} \mathbf U_\mathrm{P}^r + \sum A_\mathrm{N}\mathbf U_\mathrm{N}^n - E_\mathrm{P}^n=-\nabla p^n \end{equation} \tag{12} $$

其中 $p^n$ 为前一时间步已知压力。此方程中，除 $\mathbf U_\mathrm{P}^r$ 外，均为已知。此方程就是在 OpenFOAM 中的动量预测方程。求解此方程则可得到预测的速度 $\mathbf U_\mathrm{P}^r$。

求得预测速度 $\mathbf U_\mathrm{P}^r$ 之后，将离散方程去掉压力项并进行移项，得到 $\mathbf{U}^*$ 来表示预测的无压力项 $\mathbf U_\mathrm{P}^r$.

$$ \begin{equation} \mathbf{U}^* = \mathbf U_\mathrm{P}^r = \frac{1}{A_\mathrm{P}}\left( - \sum A_\mathrm{N} \mathbf U_\mathrm{N}^r + E_\mathrm{P}^r \right) \tag{13} \end{equation} $$

注意，预测速度 $\mathbf U_\mathrm{P}^r$ 已经从动量预测方程 $(12)$ 求出，还未涉及到连续性方程。此时预测的速度场 $\mathbf U_\mathrm{P}^r$ 并不符合连续性方程。如果我们定义 $\mathbf U_\mathrm{P}^{n+1}$ 为 PISO 循环后最终收敛的速度，那么，进行修正后的 $\mathbf U_\mathrm{P}^{n+1}$ 就应当满足动量方程和连续性方程。

速度 $\mathbf U_\mathrm{P}^{n+1}$ 符合动量方程

$$ \begin{equation} A_\mathrm{P} \mathbf U_\mathrm{P}^{n+1} + \sum A_\mathrm{N}\mathbf U_\mathrm{N}^{n+1} - E_\mathrm{P}^{n+1} = - \nabla p^{n+1} \tag{14} \end{equation} $$

其中 $p^{n+1}$ 为 PISO 循环内最后求得的压力。将公式 $(13)$ 代入到公式 $(14)$，有：

$$ \begin{equation} \mathbf U_\mathrm{P}^{n+1} = \mathbf{U}^{*,n+1} - \frac{1}{A_\mathrm{P}}\nabla p^{n+1} \end{equation}\tag{15} $$

面上的 $\mathbf U_{\mathrm{P,}f}^{n+1}$ 并不是通过体心的速度插值来获取，Rhie-Chow 插值建议这样进行获取面上的 $\mathbf U_{\mathrm{P,}f}^{n+1}$：

$$ \begin{equation} \mathbf U_{\mathrm{P},f}^{n+1} = \mathbf{U}^{*,n+1}_f - \frac{1}{A_{\mathrm{P},f}} \nabla_f p^{n+1} \end{equation}\tag{16} $$

速度 $\mathbf U_\mathrm{P}^{n+1}$ 符合连续性方程 $$ \begin{equation} \nabla \cdot \mathbf U_\mathrm{P}^{n+1}=0 \end{equation}\tag{17} $$

离散后得

$$ \begin{equation} \sum (\mathbf U_f^{n+1} \cdot \mathbf S_f)=0 \end{equation}\tag{18} $$

将获取的面心上的速度 $\mathbf U_{\mathrm{P,}f}^{n+1}$ 代入上式，

$$ \begin{equation} \sum \left[ \left( \mathbf{U}^{*,n+1}_f - \frac{1}{A_{\mathrm{P},f}}\nabla_f p^{n+1} \right) \cdot \mathbf S_f \right] = 0 \end{equation}\tag{19} $$

可将其改写成张量形式

$$ \begin{equation} \nabla \cdot (\mathbf{U}^{*,n+1}) = \nabla \cdot\left(\frac{1}{A_{\mathrm{P},f}} \nabla p^{n+1}\right) \end{equation} \tag{20} $$

PISO 循环

方程的解 $\mathbf U_\mathrm{P}^{n+1}$ 和压力 $p^{n+1}$ 应当满足连续性方程和动量方程。而到目前为止，仅有通过动量预测方程求出来的 $\mathbf U_\mathrm{P}^r$。方程 $(14)$ 和 $(19)$ 无法自行求解得到。Issa 于 1986 年提出 PISO (Pressure Implicit with Splitting of Operators) 技术来解决这个问题。其通过对当前时间步内的多次修正来获得最终的 $\mathbf U_\mathrm{P}^{n+1}$ 和 $p^{n+1}$ 。

根据 PISO 循环

icoFoam 使用 PISO 算法解 N-S 方程。其基本程序为

根据初始条件（压力场、速度场等）求解预测速度场 $\mathbf U^r$
根据预测速度求无压力梯度项的速度场 $\mathbf U^1$
根据无压力梯度项的速度场 $\mathbf U^1$，求解压力（泊松）方程，得到压力场 $p^1$
根据压力场 $p^1$，修正预测速度 $\mathbf U^r$ 以满足连续性方程
接第 2 步，直到满足收敛要求

现在，我们来探索一下源程序。（注：C++，以下源程序版本为 foam-extend 3.1 ）

在看主程序之前，我们需要了解在各求解器中需要计算的变量（或场），这个变量（或场）的创建工作由 createFields.H 来完成。请看 icoFoam 求解器的 createFields.H 头文件。

    // 在终端输出 “Reading field T”。在 OpenFOAM 中，endl 为空行，Info 表示在屏幕输出
    Info<< "Reading transportProperties\n" << endl;

    IOdictionary transportProperties // IOdictionary 意为创建字典文件
    (
        IOobject
        (
            "transportProperties",
            runTime.constant(),  // 在 constant 文件夹存储，
            mesh,                //
            IOobject::MUST_READ, // 必读文件
            IOobject::NO_WRITE 　// 从不输出，字典文件不需要输出
        )
    );

    dimensionedScalar nu // 创建变量 nu，dimensionedScalar 表示具有量纲的标量
    (
        transportProperties.lookup("nu") // 在 transportProperties 文件中读取 nu
    );

    Info<< "Reading field p\n" << endl;  // 输出 "Reading field p" 到终端，
    volScalarField p                     // 创建 volScalarField 变量 p，
    (                                    // volScalarField 是类，表示体心标量场
        IOobject  // 采用 IOobject 注册对象，变量（场）大多采用此类定义
        (
            "p", // 变量名
            runTime.timeName(),  // 在运行时间存储
            mesh,                // 注册于网格
            IOobject::MUST_READ, // 必读，如果某个场为计算得出，可以不进行读取操作
            IOobject::AUTO_WRITE // 自动输出，按 case 文件夹中　system/controlDict
        ),                       // 的定义来确定什么时候输出
        mesh 　　// 需要读取的场可以写mesh，需要计算的场则可以改写，例如：alpha1*rho1
    );


    Info<< "Reading field U\n" << endl;
    volVectorField U         // 创建变量 U，volVectorField 是类，表示的是体心向量场
    (
        IOobject
        (
            "U",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );


#   include "createPhi.H" // 包含头文件 createPhi.H 表面流率场，该文件位于
                          // src/finiteVolume/cfdTools/incompressible

    label pRefCell = 0;
    scalar pRefValue = 0.0;
    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);

以下是 icoFoam 主程序文件。

#include "fvCFD.H" 　　　　// 所有求解器必备头文件，它涉及到构建时间、组建矩阵、
                          // 有限体积离散、组建网格、量纲等大量内容。

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
#   include "setRootCase.H"        // 设置 case 的根目录，必备，
                                   // 该头文件位于 src/OpenFOAM/include
#   include "createTime.H"         // 创建时间对象 runTime
#   include "createMesh.H"         // 创建网格对象
#   include "createFields.H"       // 包含 createFields.H 头文件
#   include "initContinuityErrs.H" // 初始连续性误差

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    Info<< "\nStarting time loop\n" << endl;  // 输出“Starting time loop”
    // 计算主程序
    while (runTime.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl; // 输出时间

#       include "readPISOControls.H" // PISO 控制文件
#       include "CourantNo.H"        // Courant 数计算，文件位于
                                     // src/finiteVolume/cfdTools/incompressible
        fvVectorMatrix UEqn       // 构造速度方程
        (
            fvm::ddt(U)           // 非稳态项，隐式离散
          + fvm::div(phi, U)      // 对流项，隐式离散
          - fvm::laplacian(nu, U) // 扩散项，隐式离散
        );

        solve(UEqn == -fvc::grad(p)); // 求解速度方程

        // --- PISO loop 计算

        for (int corr = 0; corr < nCorr; corr++)
        {
            volScalarField rUA = 1.0/UEqn.A();

            U = rUA*UEqn.H();                        // 无压力项之速度值 U*
            phi = (fvc::interpolate(U) & mesh.Sf()) // U* 使用线性插值，需要在算例中设定interpolate(U)的格式
                + fvc::ddtPhiCorr(rUA, U, phi);          // 在 ext-4.0 版本中，这一项已经去掉

            adjustPhi(phi, U, p);            // 修正压力边界（梯度为0）上的 phi，使其满足连续性方程
           //网格非正交循环
            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
            {
                fvScalarMatrix pEqn
                (
                    fvm::laplacian(rUA, p) == fvc::div(phi)  // 公式 (20)，phi = Uf*Sf，此处间接使用 Rhie-Chow 插值原理
                );

                pEqn.setReference(pRefCell, pRefValue);
                pEqn.solve();               // 求解公式 (20)

                if (nonOrth == nNonOrthCorr)
                {
                    phi -= pEqn.flux();   // 使用求解的压力场修正通量场，即在最后一次修正的时候使通量守恒，
                }    // pEqn.flux()返回公式 (16) 方程右边第二项，为 fvc::interpolate(rUA)*fvc::snGrad(p)*mag(mesh.Sf())
            }    // 某些可压缩求解器中, pEqn.flux() 项可能为 '+' 号，即 phi += pEqn.flux()。这是因为 pEqn 中的 laplacian 项为 '-'号

#           include "continuityErrs.H"

            U -= rUA*fvc::grad(p);   // 公式 (16) 
            U.correctBoundaryConditions();
        }

        runTime.write();     //求解结果输出，所有 AUTO_WRITE 声明的变量，都会输出

        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }

    Info<< "End\n" << endl;    // 提示程序结束

    return 0;                  // 返回 0
}

参考文献

李东岳. icoFoam 解析. http://dyfluid.com/icoFoam.html, retrieved 2017-06-22
H Jasak. Consider a case which has zero. https://www.cfd-online.com/Forums/openfoam-solving/59668-two-fundamental-questions-about-icofoam-while-updating-velocities-pressure.html#post190879 retrieved on 2017-07-23

Main menu

Error message

icoFoam 解析

动量预测

PISO 循环

Main menu

Error message

You are here

icoFoam 解析

动量预测

PISO 循环