• 返回主页
  • icoFoam解析

    2017.03.22:增加fvc::div的前后量纲关系的讨论链接;
    2017.03.08:统一S符号;增加Rhie-Chow原理对应;
    2016.11.23:统一HbyA符号;
    2016.05.20:有关adjustphi()的描述;规范公式化书写;强调PISO算法
    2016.05.09:有关div(phiHbyA)的描述;
    2016.04.17:有关fvc::reconstruct(phi)的描述;
    求解器类型:瞬态层流牛顿流体求解器。

    icoFoam采用的是“无重力”、“同位网格”、“瞬态PISO”算法来求解“不可压缩”、“牛顿流体”的“层流”流动。其为OpenFOAM中最简单的N-S方程求解器。本文从最基本的N-S方程在笛卡尔网格下的有限体积离散开始推导,并和OpenFOAM植入的代码相对应以供理解瞬态PISO算法如何在OpenFOAM中实现以及N-S方程如何在OpenFOAM中的离散。对于Patanka的SIMPLE算法以及Henry发明的PIMPLE算法的对比请参考其他文献$^{1,2}$。对于无重力可压缩流体求解器请参考rhoPimpleFoam解析,附加重力的可压缩流体求解器请参考buoyantPimpleFoam解析

    方程离散

    首先对于不可压缩流体并忽略体积力源项我们有动量方程: \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$为密度。忽略压力梯度项(在后文中即对应动量预测步骤),对方程($\ref{mom}$)进行全隐半离散化: \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$为动力粘度。公式($\ref{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}$其在三维情况下为一个矢量。且方程($\ref{sngrad}$)只对矩形网格精准。将公式($\ref{mom1}$)-($\ref{sngrad}$)代入公式($\ref{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} 需要注意的是:在方程($\ref{mom2}$)对流项的离散中,其中通量,$F_f$,采用当前已知时间步的速度,$\mathbf{U}^n$,来计算。并保留另一个预测速度,$\mathbf{U}^r_\mathrm{P}$,为未知量。这种把其中一个速度已知化并保留另外一个未知速度的过程即为线性化:使用当前时间步的速度对高阶非线性对流项进行离散。由于NS方程是非线性化的,要么选择非线性求解器,要么将对流项线性化$^3$。由于非线性求解器非常复杂,在OpenFOAM均采用线性化处理。线性化的一个问题即为通量速度的信息有所滞后$^3$。

    另外需要注意的是在方程($\ref{momF}$)中,并没有考虑压力的作用。有关方程$\eqref{momF}$中具体的如何使用高斯定理以及相关量纲问题可参考:fvc::div的前后量纲关系,在此不做深入介绍。

    动量预测

    在公式($\ref{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} 将公式($\ref{sngrad}$),($\ref{linear}$)代入到公式($\ref{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}$保持不变。

    我们在公式($\ref{apanmom}$)中加入上个时间步已知的未离散的压力梯度项。注意:由于为不可压缩流体,在icoFoamp(单位:m$^2$s$^{-2}$)即为$p/\rho$。在下文中的$p$对应为公式(1)中的$p/\rho$即为: \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= -\nabla p^n \label{apanmomP} \end{equation} 其中$p^n$为当前时间步已知的压力。仔细观察公式($\ref{apanmomP}$),除了$\mathbf{U}_\mathrm{P}^r$外均为已知。此方程即为OpenFOAM中大名鼎鼎的动量预测方程,求解此方程我们有预测速度$\mathbf{U}_\mathrm{P}^r$。

    压力泊松方程

    在求得预测速度$\mathbf{U}^r_\mathrm{P}$之后。再次去掉公式($\ref{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}^r} \right) \label{hbya} \end{equation} 注意,预测速度$\mathbf{U}_\mathrm{N}^r$已经从动量预测方程($\ref{apanmomP}$)求出,目前并未涉及到连续性方程。此时预测的速度场$\mathbf{U}^r$并不符合连续性方程。如果我们定义$\mathbf{U}^{n+1}$为PISO循环后最终收敛的速度,那么在进行修正后的$\mathbf{U}_\mathrm{P}^{n+1}$:


    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}^{n+1} = - \nabla {p^{n+1}} \label{apanmomPn1} \end{equation} 其中$p^{n+1}$为PISO循环内最后求得的压力。将公式($\ref{hbya}$)代入到公式($\ref{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插值建议这样进行获取面上的$\mathbf{U}_{\mathrm{P},f}^{n+1}$: \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} 2.应该符合连续性方程: \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} 将公式($\ref{uf}$)代入到公式($\ref{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}} \right) \cdot \bfS_f} = \sum {\left( {\frac{1}{{A_{\mathrm{P},f}}}{\nabla _f}{p^{n+1}}} \right) \cdot \bfS_f} \end{equation} 也即为: \begin{equation} \nabla \cdot (\mathbf{HbyA}^{n+1}) = \nabla \cdot(\frac{1}{A_{\mathrm{P},f}} \nabla p^{n+1}) \label{poss} \end{equation}

    PISO循环

    总结我们的限制条件即为$n+1$时间步的速度$\mathbf{U}_\mathrm{P}^{n+1}$和压力$p^{n+1}$应该符合方程($\ref{apanmomPn1}$)和($\ref{poss}$)。注意:在目前我们已知的只有通过动量预测求出来的$\mathbf{U}_\mathrm{P}^r$,$\mathbf{U}_\mathrm{P}^{n+1}$和$p^{n+1}$均为未知的。方程方程($\ref{apanmomPn1}$)和($\ref{poss}$)均为未知,这又回到了NS方程最初离散后的求解过程,$\mathbf{U}_\mathrm{P}^r$一无所用。1985年Issa提出了一种Pressure Implicit with Splitting of Operators (PISO) 求解技术来解决这个问题$^4$,PISO算法在当时主要针对不可压缩非稳态计算,其为一种非迭代的算法。在随后也被拓展用于稳态问题中。其通过对当前时间步内的多次修正来获得最终的$\mathbf{U}_\mathrm{P}^{n+1}$和$p^{n+1}$。

    如何更新速度?

    依据PISO循环,我们改写方程($\ref{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}^r = - \nabla {p^{(1)}} \label{pisoapan} \end{equation} 对比方程($\ref{apanmomPn1}$),方程($\ref{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}$是相同的$^3$。参考方程($\ref{uf}$),方程($\ref{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} 不难发现,在方程($\ref{pisouf}$)中,在求得$p^{(1)}$之后,可以更新速度$\mathbf{U}_{\mathrm{P}}^{(1)}$。更新的速度符合动量方程。但是,如何获取$p^{(1)}$?

    如何求解$p^{(1)}$?

    我们改写方程($\ref{poss}$)有: \begin{equation} \nabla \cdot (\mathbf{HbyA}^{r}) = \nabla \cdot(\frac{1}{A_{\mathrm{P},f}} \nabla p^{(1)}) \label{pisoposs} \end{equation} 可以发现方程($\ref{pisoposs}$)即为我们通过连续性方程组建的压力泊松方程,求解后我们有$p^{(1)}$。并且此步骤进入了重要的Rhie-Chow插值原理来消除压力波:Rhie-Chow插值原理。求解之后需要注意以下几点:

    1. 将$p^{(1)}$回代到方程($\ref{pisouf}$)我们有$\mathbf{U}_{\mathrm{P}}^{(1)}$。我们知道$\mathbf{U}_{\mathrm{P}}^{(1)}$符合动量方程,但是更新后的速度又不符合连续性方程。
    2. 仔细观察方程($\ref{pisouf}$),我们发现使用的$\mathbf{HbyA}$滞后一内迭代步(使用的是$\mathbf{HbyA}^r$而不是$\mathbf{HbyA}^{(1)}$)。这种滞后性在瞬态PISO算法中每个时间步内最终收敛的情况下是可以忽略的$^3$。

    下面我们简述PISO循环中的迭代步骤:

     1.依据初始条件求解预测速度$\mathbf{U}^r$,在进入到PISO循环中我们称之为$\mathbf{U}^{(1)}$;
     2.依据速度$\mathbf{U}^{(1)}$组建$\mathbf{HbyA}^{(1)}$;
     3.求解方程($\ref{pisoposs}$)求解压力$p^{(1)}$;
     4.通过方程($\ref{pisouf}$),依据压力$p^{(1)}$更新速度有$\mathbf{U}^{(1)}$;
     5.回到第二步,执行多次PISO循环;

    下面进入代码分析:

    #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);//此函数功能请参考:更新4
    
                for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
                {
                    fvScalarMatrix pEqn
                    (
                        fvm::laplacian(rAU, p) == fvc::div(phiHbyA)//公式(26),有关div(phiHbyA)请参考:更新3
                    );
    
                    pEqn.setReference(pRefCell, pRefValue);
                    pEqn.solve();//公式(26)
    
                    if (nonOrth == nNonOrthCorr)
                    {
                        phi = phiHbyA - pEqn.flux();//使用求解的压力场修正通量场,在最后一次修正的时候通量守恒,Issa指出,大约需要2-3次内循环步。
    对应公式(25),pEqn.flux()返回公式(25)方程右边第二项,也为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)会引起比较大的守恒误差。目前这种误差有多大是未知的。更多相关的讨论请参考:更新1更新2
    U.correctBoundaryConditions(); } runTime.write(); Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << " ClockTime = " << runTime.elapsedClockTime() << " s" << nl << endl; } Info<< "End\n" << endl; return 0; }

    参考文献
    [1]. Patankar S.V., Spalding D.B. A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows[J]. International journal of heat and mass transfer, 1972, 15(10): 1787-1806.
    [2]. 李东岳. 浅谈Piso,Simple算法[OL]. http://dyfluid.com/PISOSIMPLE.html, 2016.3.10
    [3]. Jasak H. Error analysis and estimation for the finite volume method with applications to fluid flows[D]. Imperial College London (University of London), 1996.
    [4]. Issa R.I. Solution of the implicitly discretised fluid flow equations by operator-splitting[J]. Journal of computational physics, 1986, 62(1): 40-65.


    东岳流体dyfluid.com
    勘误、讨论、补充内容请前往CFD中国