return 0;
wmake

icoFoam解析

如果打不开图像,请右键在新标签页打开图像后刷新几次 如果打不开图像,请右键在新标签页打开图像后刷新几次


引言

icoFoam为OpenFOAM中的一个求解器,其求解的为不可压缩牛顿流体的NS方程,可用于层流模拟或流场剪切力引致的湍流quasi-DNS直接模拟(Komen et al., 2014; Axtmann and Rist, 2016; Liu et al., 2016)。在这个求解器中,压力和速度的耦合通过瞬态PISO算法进行计算(Issa, 1986),且没有考虑其他体积力(如重力等)。相比较OpenFOAM中的其他求解器,除去三个基本求解器(laplacianFoam, potentialFoam, scalarTransportFoam),icoFoam为一个最基本的求解器范例,可用于理解OpenFOAM中压力速度耦合的框架策略。本文从最基本的NS方程在笛卡尔网格下的有限体积离散开始推导,并和OpenFOAM植入的代码相对应,帮助用户理解NS方程如何在有限体积法中的离散以及瞬态PISO算法在OpenFOAM中的实现。对于流场内温度变化较大,但可忽略浮力以及重力的可压缩求解器,推荐选用rhoPimpleFoam。附加重力的可压缩浮力驱动流体求解器,推荐选用buoyantPimpleFoam

方程离散以及求解

NS方程中压力速度的耦合求解分为方程离散步、动量预测步、压力泊松方程组建、速度更新步。本节主要讨论icoFoam中的方程离散以及求解步骤。

方程离散

对于不可压缩流体,并忽略体积力源项的动量方程可表示为(Ferziger and Peric, 2012)

\begin{equation} \frac{\partial \mathbf{U}}{\partial t}+\nabla \cdot (\mathbf{U} \otimes\mathbf{U})=-\nabla \frac{p}{\rho}+\nabla \cdot(\nu \nabla \mathbf{U}), \label{mom} \end{equation}

其中的$\mathbf{U}$为速度矢量,$p$为压力,$\rho$为密度。忽略压力梯度项(在后文中即对应动量预测步骤),对方程\eqref{mom}进行欧拉全隐离散有(Jasak, 1996)

\begin{equation} \int \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}}, \label{mom1} \end{equation}
\begin{equation} \int \int {\nabla \cdot \left( {\mathbf{U}\otimes\mathbf{U}} \right)\mathrm{d}V\mathrm{d}t = \int \int {\mathbf{U}\otimes\mathbf{U}\mathrm{d}\bfS\mathrm{d}t} = } \sum {{{\left( {{\mathbf{U}^n}\otimes{\mathbf{U}^r}} \right)}_f}} \bfS_f\Delta t = \sum {\left( {F_f^n \mathbf{U}_f^r} \right)} \Delta t, \label{mom2} \end{equation}
\begin{equation} {\int \int \nabla \cdot \left( {\nu \nabla \mathbf{U}} \right)\mathrm{d}V\mathrm{d}t = \int \int {\nu \nabla \mathbf{U}\mathrm{d}\bfS\mathrm{d}t = \sum {\left( {\nu \nabla {\mathbf{U}^r}} \right)} } _f}\bfS_f\Delta t = \sum {\left( {\nu \left( {\nabla _f^ \bot {\mathbf{U}^r}} \right)} \right)} \Delta t, \label{mom3} \end{equation}

其中上标$^n$表示为当前的时间步(已知),上标$^r$表示预测时间步(待求),下标$_f$表示网格单元面上的值,$V_\rP$表示网格单元体积,$\bfS$表示网格单元的各个面的面矢量,$\Delta t$表示时间步长,$F_f$为通量,$\nu$为动力粘度。方程\eqref{mom3}中$\nabla _f^ \bot {\mathbf{U}^r}$(面法向梯度)在正交网格中进一步可以表示为:

\begin{equation} \nabla _f^ \bot {\mathbf{U}^r}{\rm{ = }}{\left( {\nabla {\mathbf{U}^r}} \right)_f}\bfS_f = \left| \bfS_f \right|\frac{{\mathbf{U}_\mathrm{N}^n - \mathbf{U}_\mathrm{P}^r}}{{\left| d \right|}}, \label{sngrad} \end{equation}

其中$d$表示网格单元体心之间的距离,下标$_\mathrm{N}$表示相邻网格单元的速度,下标$_\mathrm{P}$表示当前网格单元的速度。$\nabla _f^ \bot {\mathbf{U}^n}$在三维情况下为一个矢量。在这里需要注意的是,方程\eqref{sngrad}只对矩形网格精准,网格非正交性增加,需要特定的数值处理提高其计算精准性(Demirdžić,2015)。同时,在本篇文章中,Demirdžić对不同的非正交修正算法进行综述并对审稿人以及编辑进行批判(对怼顶级期刊JCP,学术界照样撕)。

将方程\eqref{mom1}-\eqref{sngrad}代入到方程\eqref{mom}有:

\begin{equation} \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)} . \label{momF} \end{equation}

需要注意的是,方程\eqref{mom2}中的左边第二项(对流项)是非线性的。在求解的过程中,要么选择非线性求解器,要么将对流项线性化。由于非线性求解器非常复杂,因此在OpenFOAM均采用线性化处理。具体的,在对方程\eqref{mom2}中的左边第二项对流项离散中,其中的通量$F_f$采用当前已知时间步的速度$\mathbf{U}^n$来计算,同时保留预测速度$\mathbf{U}^r_\mathrm{P}$为未知量。这种将高阶非线性项降为一阶线性项的过程即为线性化操作(Versteeg and Malalasekera, 2007; 陶文铨, 2001)。线性化的一个问题即为通量速度的信息有所滞后(Jasak, 1996)。方程\eqref{momF}中具体的如何使用高斯定理以及相关量纲问题可参考其他讨论,在此不做介绍(fvc::div的前后量纲关系)。

动量预测

方程\eqref{momF}中$\mathbf{U}_f^r$被定义为面上的预测速度,然而在计算中求得的均为网格单元体心上的速度。因此需要从体心速度进行插值来获得面速度,在此步可以引入各种插值格式。假设使用中心差分:

\begin{equation} \mathbf{U}_f^r = \frac{{\mathbf{U}_\mathrm{P}^r + \mathbf{U}_\mathrm{N}^n}}{2}. \label{linear} \end{equation}

将方程\eqref{sngrad},\eqref{linear}代入到方程\eqref{momF}有

\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| \bfS_f \right|\frac{{\mathbf{U}_\mathrm{N}^n - \mathbf{U}_\mathrm{P}^r}}{{\left| d \right|}}} \right)}, \end{equation}
\begin{equation} \left( {\frac{{{V_\mathrm{P}}}}{{\Delta t}} + \sum {\left( {\frac{{F_f^n}}{2}} \right) + \sum {\left( {\nu \frac{{\left| \bfS_f \right|}}{{\left| d \right|}}} \right)} } } \right)\mathbf{U}_\mathrm{P}^r = - \sum {\left( {\left( {\frac{{F_f^n}}{2} - \nu \frac{{\left| \bfS_f \right|}}{{\left| d \right|}}} \right)\mathbf{U}_\mathrm{N}^n} \right)} + \frac{{{V_p}}}{{\Delta t}}\mathbf{U}_\mathrm{P}^n. \end{equation}

将上式简化后有

\begin{equation} {A_\mathrm{P}}\mathbf{U}_\mathrm{P}^r{\rm{ + }}\sum {A_\mathrm{N}\mathbf{U}_\mathrm{N}^n} {\rm{ - }}E_\mathrm{P}^n=0, \label{apanmom} \end{equation}

其中

\begin{equation} A_\mathrm{P}=\left( {\frac{{{V_p}}}{{\Delta t}} + \sum {\left( {\frac{{F_f^n}}{2}} \right) + \sum {\left( {\nu \frac{{\left| \bfS_f \right|}}{{\left| d \right|}}} \right)} } } \right), \end{equation}
\begin{equation} A_\mathrm{N}=\left( {\frac{{F_f^n}}{2} - \nu \frac{{\left| \bfS_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}

可以看出在某个时间步内$A_\mathrm{P}$和$A_\mathrm{N}$保持不变。$E_\mathrm{P}$的值不变。需要注意的是,由于为不可压缩流体求解器,icoFoam中p的单位为压力除以密度($p/\rho$)的单位,即$\left[\mathrm{m}^2\cdot \mathrm{s}^{-2}\right]$。随后在下文中,$p$对应方程\eqref{mom}中的$p/\rho$。现在,在方程\eqref{apanmom}中加入上个时间步已知的未离散的压力梯度项:

\begin{equation} {A_\mathrm{P}}\mathbf{U}_\mathrm{P}^r{\rm{ + }}\sum {A_\mathrm{N}\mathbf{U}_\mathrm{N}^n} {\rm{ - }}E_\mathrm{P}= -\nabla p^n, \label{apanmomP} \end{equation}

其中$p^n$为当前时间步已知的压力。方程\eqref{apanmomP}中除了$\mathbf{U}_\mathrm{P}^r$外均为已知,求解此方程则有预测速度$\mathbf{U}_\mathrm{P}^r$。方程\eqref{apanmomP}即为OpenFOAM中的动量预测方程。

压力泊松方程

在求得预测速度$\mathbf{U}^r_\mathrm{P}$之后,再次去掉方程\eqref{apanmomP}中的压力项并进行移项。下面,定义一个速度并使用符号$\mathbf{HbyA}^r$来表示:

\begin{equation} \mathbf{HbyA}^r = \mathbf{U}_\mathrm{P}^r = \frac{1}{{{A_\mathrm{P}}}}\left( { - \sum {{A_\mathrm{N}}\mathbf{U}_\mathrm{N}^r} + E_\mathrm{P}} \right). \label{hbya} \end{equation}

由于预测速度$\mathbf{U}_\mathrm{N}^r$已经从动量预测方程\eqref{apanmomP}中求出,但目前并未涉及到连续性方程,因此此时预测的速度场$\mathbf{U}^r$(或$\mathbf{HbyA}^r$)并不符合连续性方程。如果定义$\mathbf{U}^{n+1}$为PISO循环后最终收敛的速度,那么在进行修正后的$\mathbf{U}_\mathrm{P}^{n+1}$一方面应该符合动量方程,即

\begin{equation} {A_\mathrm{P}}\mathbf{U}_\mathrm{P}^{n+1}{\rm{ + }}\sum {{A_\mathrm{N}}\mathbf{U}_\mathrm{N}^{n+1}} {\rm{ - }}E_\mathrm{P} = - \nabla {p^{n+1}}, \label{apanmomPn1} \end{equation}

其中$p^{n+1}$为PISO循环内最后求得的压力。将方程\eqref{hbya}代入到方程\eqref{apanmomPn1}有

\begin{equation} \mathbf{U}_\mathrm{P}^{n+1} = \mathbf{HbyA}^{n+1} - \frac{1}{{{A_\mathrm{P}}}}\nabla {p^{n+1}}. \end{equation}

常规方法下,面上的$\mathbf{U}_{\mathrm{P},f}^{n+1}$并不是通过体心的速度插值来获取,而是通过Rhie-Chow插值来进行获得(Rhie and Chow, 1983):

\begin{equation} \mathbf{U}_{\mathrm{P},f}^{n+1} = \mathbf{HbyA}_f^{n+1} - \frac{1}{{{A_{\mathrm{P},f}}}}\nabla_f {p^{n+1}}. \label{uf} \end{equation}

另一方面,$\mathbf{U}_\mathrm{P}^{n+1}$应该符合连续性方程:

\begin{equation} \nabla \cdot \mathbf{U}_\mathrm{P}^{n+1}=0. \end{equation}

其离散后为

\begin{equation} \sum (\mathbf{U}_{f}^{n+1} \cdot \bfS_f)=0. \label{cont} \end{equation}

将方程\eqref{uf}代入到方程\eqref{cont}有:

\begin{equation} \sum \left( \left( \mathbf{HbyA}_f^{n+1} - \frac{1}{{{A_{\mathrm{P},f}}}}\nabla_f {p^{n+1}} \right) \cdot \bfS_f \right)=0, \end{equation}
\begin{equation} \sum {\left( {\mathbf{HbyA}_f^{n+1}} \cdot \bfS_f\right) } = \sum {\left( {\frac{1}{{A_{\mathrm{P},f}}}{\nabla _f}{p^{n+1}}} \cdot \bfS_f\right) }, \label{jim} \end{equation}
\begin{equation} \nabla \cdot (\mathbf{HbyA}^{n+1}) = \nabla \cdot(\frac{1}{A_{\mathrm{P}}} \nabla p^{n+1}). \label{poss} \end{equation}

方程\eqref{poss}即为压力泊松方程。

PISO循环

可以看出,最后求解的速度$\mathbf{U}_\mathrm{P}^{n+1}$和压力$p^{n+1}$应该符合方程\eqref{apanmomPn1}和\eqref{poss},分别对应动量方程和连续性方程。然而目前,我们只有通过2.1节动量预测步骤求出来的$\mathbf{U}_\mathrm{P}^r$。$\mathbf{U}_\mathrm{P}^{n+1}$和$p^{n+1}$均为未知的。方程方程($\ref{apanmomPn1}$)和($\ref{poss}$)均为未知。1988年Issa提出了PISO算法 (Pressure Implicit with Splitting of Operators)来解决这个问题(Issa, 1986)。PISO算法在提出时主要针对不可压缩非稳态计算,其为一种非迭代的算法(Issa, 1986),在随后也被拓展用于稳态问题中(Versteeg and Malalasekera, 2007)。其通过对当前时间步内的多次修正来获得最终的$\mathbf{U}_\mathrm{P}^{n+1}$和$p^{n+1}$。

依据PISO算法,改写方程\eqref{apanmomPn1}为

\begin{equation} {A_\mathrm{P}}\mathbf{U}_\mathrm{P}^{(1)}{\rm{ + }}\sum {{A_\mathrm{N}}\mathbf{U}_\mathrm{N}^r} {\rm{ - }}E_\mathrm{P} = - \nabla {p^{(1)}}. \label{pisoapan} \end{equation}

对比原方程\eqref{apanmomPn1},方程\eqref{pisoapan}将$\mathbf{U}_\mathrm{P}^{n+1}$和$p^{n+1}$代替为$\mathbf{U}_\mathrm{N}^{(1)}$和$p^{(1)}$且保持为未知。$\mathbf{U}_\mathrm{N}^{n+1}$代替为$\mathbf{U}_\mathrm{N}^{r}$为已知。这是因为在压力循环步收敛的时候,$\mathbf{U}_\mathrm{N}^{(n)}$和$\mathbf{U}_\mathrm{N}^{n+1}$是相同的(Issa, 1986)。也正是在方程\eqref{pisoapan}中,其引入了SIMPLE算法中略去临点修正的思想,这可以通过将方程\eqref{pisoapan}和方程\eqref{apanmomP}相减来证实:

\begin{equation} A_\mathrm{P}\bfU'=- \nabla p', \label{uprime} \end{equation}

其中$\bfU'=\mathbf{U}_\mathrm{P}^{(1)}-\mathbf{U}_\mathrm{P}^{r}$,$p'=p^{(1)}-p^n$。

参考方程\eqref{uf},方程\eqref{pisoapan}可整理为

\begin{equation} \mathbf{U}_{\mathrm{P},f}^{(1)} = \mathbf{HbyA}_f^{r} - \frac{1}{{{A_{\mathrm{P},f}}}}\nabla_f {p^{(1)}}. \label{pisouf} \end{equation}

在方程\eqref{pisouf}中,求得$p^{(1)}$之后,可以更新速度$\mathbf{U}_{\mathrm{P}}^{(1)}$,但其中的$\mathbf{HbyA}_f^{r}$只能采用之前的时间步,这样带来了更新的滞后,为后续的多次求解埋下铺垫。另外,如何获取$p^{(1)}$?现在改写方程\eqref{poss}:

\begin{equation} \nabla \cdot (\mathbf{HbyA}^{r}) = \nabla \cdot(\frac{1}{A_{\mathrm{P},f}} \nabla p^{(1)}). \label{pisoposs} \end{equation}

可以看出,方程\eqref{pisoposs}即为通过连续性方程组建的压力泊松方程,求解后有$p^{(1)}$。并且此步骤进入了重要的Rhie-Chow插值原理来消除压力棋盘分布(Rhie and Chow, 1983; 陶文铨, 2001)。需要注意的是,将$p^{(1)}$回代到方程\eqref{pisouf}的时候,有$\mathbf{U}_{\mathrm{P}}^{(1)}$,但依然并不能严格满足动量方程,因为$\mathbf{HbyA}_f^{r}$的信息有所滞后;因此需要要通过压力泊松方程对$\mathbf{HbyA}_f^{r}$进行再次更新。通过求解多次的压力泊松方程,这种滞后性逐渐消除并可最终忽略(Issa, 1986)。

总而言之,PISO算法中的的迭代过程可以表示为下面几个步骤:

 1.依据初始条件求解预测速度$\mathbf{U}^r$,在这里使用$\mathbf{U}^*$表示;

 2.依据速度$\mathbf{U}^*$组建$\mathbf{HbyA}^*$;

 3.求解方程\eqref{pisoposs}求解压力$p^*$;

 4.通过方程\eqref{pisouf},依据压力$p^*$以及$\mathbf{HbyA}^*$更新速度有$\mathbf{U}^{**}$;

 5.回到第二步,重新组建$\mathbf{HbyA}^{**}$,并求解压力$p^{**}$,更新速度有$\mathbf{U}^{***}$;

代码分析

#include "fvCFD.H" //有关头文件等参考laplacianFoam, potentialFoam, scalarTransportFoam解析

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

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"
    #include "createFields.H"
    #include "initContinuityErrs.H"

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

    Info<< "\nStarting time loop\n" << endl;

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

        #include "readPISOControls.H"
        #include "CourantNo.H"

        fvVectorMatrix UEqn
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
        ); //对应方程(6),(10)

        solve(UEqn == -fvc::grad(p)); //对应方程(14)

        // --- PISO loop

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

            volVectorField HbyA("HbyA", U);
            HbyA = rAU*UEqn.H(); //方程(15)
            surfaceScalarField phiHbyA
            (
                "phiHbyA",
                (fvc::interpolate(HbyA) & mesh.Sf())
                //此处依据Rhie-chow插值原理,HbyA使用线性插值得到,
                //即需要在算例中设定interpolate(HbyA)的格式
              + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
            );
            //phiHbyA即为以HbyA来计算的通量,第一次内循环的时候,
            //此处phiHbyA并不守恒,not divergence free

            adjustPhi(phiHbyA, U, p);

            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
            {
                fvScalarMatrix pEqn
                (
                    fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
                    //方程\eqref{pisoposs}
                );

                pEqn.setReference(pRefCell, pRefValue);
                pEqn.solve(); //方程\eqref{pisoposs}

                if (nonOrth == nNonOrthCorr)
                {
                    phi = phiHbyA - pEqn.flux();
                    //使用求解的压力场修正通量场,在最后一次修正的时候通量守恒,
					Issa指出,大约需要2-3次内循环步。
                    //对应方程\eqref{pisouf},pEqn.flux()返回方程\eqref{pisouf}右边第二项,
					也为fvc::interpolate(rUA)*fvc::snGrad(p)*mag(mesh.Sf())。
                    //某些可压缩求解器其中的pEqn.flux()可能为+号,
					即为phi = phiHbyA + pEqn.flux()。
					这是因为pEqn中的laplacian项为−号
                }
            }

            #include "continuityErrs.H"

            U = HbyA - rAU*fvc::grad(p);
            //方程(17),这里并没有采用重组守恒通量phi的方法来计算速度,即U=fvc::reconstruct(phi);
            //Henry表示采用fvc::reconstruct(phi)会引起比较大的守恒误差。目前这种误差有多大是未知的。

            U.correctBoundaryConditions();
        }

        runTime.write();

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

    Info<< "End\n" << endl;

    return 0;
}

验证算例

待更新

参考文献

Komen, E., Shams, A., Camilo, L., Koren, B., 2014. Quasi-DNS capabilities of OpenFOAM for different mesh types. Computers & Fluids, 96, 87-104.

Axtmann, G., Rist, U., 2016. Scalability of OpenFOAM with Large Eddy Simulations and DNS on High-Performance Systems. In High Performance Computing in Science and Engineering, 16, 413-424. Springer.

Liu, Q., Gómez, F., Pérez, J. M., Theofilis, V., 2016. Instability and sensitivity analysis of flows using OpenFOAM. Chinese Journal of Aeronautics, 29, 316-325.

Issa, R. I., 1986. Solution of the implicitly discretised fluid flow equations by operator-splitting. Journal of Computational Physics, 62, 40-65.

Ferziger, J. H., Peric, M., 2012. Computational methods for fluid dynamics. Springer Science Business Media.

Jasak, H., 1996. Error analysis and estimation for finite volume method with applications to fluid flow. PhD Thesis. Imperial College London.

Demirdžić, I., 2015. On the discretization of the diffusion term in finite-volume continuum mechanics. Numerical Heat Transfer, Part B: Fundamentals, 68, 1-10.

对怼顶级期刊JCP!学术界照样撕!

Versteeg, H. K., Malalasekera, W., 2007. An introduction to computational fluid dynamics: the finite volume method. Pearson Education.

陶文铨., 2001. 数值传热学. 西安交通大学出版社.

fvc::div的前后量纲关系,http://cfd-china.com/topic/703

Rhie, C. M., Chow, W. L., 1983. Numerical study of the turbulent flow past an airfoil with trailing edge separation. AIAA journal, 21, 1525-1532.

更新历史
2018.07.05依据尹胜强建议修正方程26、27的对应关系 | 2018.05.21依据JimShi建议将方程22中的$\bfS$移至括号内 | 2018.05.08依据张某人建议去掉方程\eqref{poss}中的下标$_f$ | 2018.04.15依据冯强老师建议增加参考文献、更正错别字 | 2017.08.13增加方程\eqref{uprime}

东岳流体®版权所有
勘误、讨论、补充内容请前往CFD中文网

');