rhoCentralFoam解析

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


引言

气体动力学的欧拉方程主要用于描述无粘的高超音速流动。相比较与抛物线特征的Navier-Stokes方程,欧拉方程具有双曲的数学特征。简化的常系数的线性欧拉方程可将其解耦为若干对流方程求解。但通常情况下流体动力学方程都是非线性的,且积分形式的欧拉方程可以存在间断解,这种间断并带有初值的问题通常被称之为黎曼问题。精确地捕捉黎曼问题按照时间步进的间断需要特殊的数值处理。一种方法是在每个网格交界面的间断处调用黎曼求解器,但这种方法计算量较大。另一种方法是依据流动的方向采用Godunov方法进行有限体积求解,进一步的可以调用高分辨率TVD格式。还有一种方法并不需要判断流动的方向,这一类方法在高超音速领域一般被称之为中心格式。中心格式不需要判断流动的方向,也不需要黎曼求解器,因此植入简单。典型的中心格式为Lax-Friedrichs方法,其在无条件不稳定格式的基础上添加一定的耗散项,使得结果无震荡。但传统的Lax-Friedrichs方法耗散项的系数较大,虽能保证求解稳定但耗散严重。Kurganov and Tadmor(KT)格式在中心格式的基础上,考虑波的传输速度,并依据波速确定更小的耗散系数,来获得更高的精度。本文针对欧拉方程,简述密度基求解框架,重点讨论KT格式在密度基框架下的离散。

考虑下图左侧的错位Lax-Friedrichs方法,纵轴为时间$t$,横轴为$x$,错位Lax-Friedrichs方法对红框区域做积分,由于红框区域将内部的局部黎曼扇进行了包含,因此错位Lax-Friedrichs方法不需要求解黎曼问题。但由于错位Lax-Friedrichs方法做的积分区域较大,导致精度降低。右侧为KT格式在错位Lax-Friedrichs方法的基础上,考虑波的传输速度。针对网格点$i$,对$i$内的未收波动影响的红色方框,附加左侧红色方框的右侧区域以及右侧方框的左侧区域进行积分,将黎曼扇的影响区域缩小,增加格式的精度。由于黎曼扇的区域取决于波的传输速度,因此KT格式需要附加波速变量,这也是唯一附加的一个变量。可见,KT格式在有限体积法中植入非常简单。

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

一维Lax-Friedrichs方法和KT格式

考虑一维常系数对流方程

\begin{equation} \frac{\partial \rho}{\partial t}+ u\frac{\p \rho}{\p x}=0 \label{adv} \end{equation}

其中$\rho$表示传输的密度,$u$表示$x$方向传输速度。对其进行显性离散有(省略当前时间步上标$^n$)

\begin{equation} \int_x\int_y\int_z\int_t\frac{\partial \rho}{\partial t}\rd t\rd x\rd y\rd z+ \int_x\int_y\int_z\int_tu\frac{\p \rho}{\p x}\rd x\rd t\rd y\rd z=0 \end{equation}
\begin{equation} (\rho_i^{n+1}-\rho_i)\Delta x\Delta y\Delta z+\left(u_{i+\frac{1}{2}}\rho_{i+\frac{1}{2}}-u_{i-\frac{1}{2}}\rho_{i-\frac{1}{2}}\right)\Delta t\Delta y\Delta z=0 \label{adv2} \end{equation}

其中$i$表示当前网格体心,$i+\frac{1}{2}$和$i-\frac{1}{2}$表示网格单元左右面。令

\begin{equation} q\Delta x\Delta y\Delta z=m,\Delta y\Delta z=|\bfS_f| \end{equation}

其中$m$表示体积为$\Delta x\Delta y\Delta z$密度为$\rho$的液体质量,有

\begin{equation} (m_i^{n+1}-m_i)+\left(u_{i+\frac{1}{2}}\rho_{i+\frac{1}{2}}-u_{i-\frac{1}{2}}\rho_{i-\frac{1}{2}}\right)\Delta t|\bfS_f|=0 \end{equation}
\begin{equation} m_i^{n+1}=m_i-\Delta t\left(u_{i+\frac{1}{2}}|\bfS_f|\rho_{i+\frac{1}{2}}-u_{i-\frac{1}{2}}|\bfS_f|\rho_{i-\frac{1}{2}}\right) \label{adv3} \end{equation}

当前时间布的$m_i$已知,若有网格单元面上的$\rho_{i+\frac{1}{2}}$和$\rho_{i-\frac{1}{2}}$,则可求下一个时间步的结果。如果认为网格单元面上的值为网格面左右网格体心值的平均,即

\begin{equation} \rho_{i+\frac{1}{2}}=\frac{1}{2}\left(\rho_i+\rho_{i+1}\right) \label{c1} \end{equation}
\begin{equation} \rho_{i-\frac{1}{2}}=\frac{1}{2}\left(\rho_{i-1}+\rho_{i}\right) \label{c2} \end{equation}

若将方程\eqref{c1}和\eqref{c2}代入到\eqref{adv3}会发现其是无条件不稳定的。Lax-Friedrichs方法在无条件不稳定格式的基础上,将其改写为

\begin{equation} \rho_{i+\frac{1}{2}}=\frac{1}{2}\left(\rho_i+\rho_{i+1}\right)-\frac{1}{2}\left(\rho_{i+1}-\rho_i\right) \label{lax1} \end{equation}
\begin{equation} \rho_{i-\frac{1}{2}}=\frac{1}{2}\left(\rho_{i-1}+\rho_{i}\right)-\frac{1}{2}\left(\rho_{i}-\rho_{i-1}\right) \label{lax2} \end{equation}

方程\eqref{lax1}和\eqref{lax2}右边的后两项可以认为是扩散系数为$\frac{1}{2}$的数值扩散项。正是因为这俩项,才使得Lax-Friedrichs方法稳定。

Lax-Friedrichs方法添加的数值扩散项系数并不是最优的,在网格分辨率较低的情况下,扩散系数$\frac{1}{2}$较大,耗散严重。是否可以将这个扩散系数降低来降低格式的耗散性呢?答案是肯定的。Kurganov and Tadmor(KT)格式在Lax-Friedrichs方法的基础上,考虑可压缩流动中波的传输速度,并依据波速确定更小的耗散系数,来获得更高的精度。KT格式认为网格界面的通量这样计算:

\begin{equation} \rho_{i+\frac{1}{2}}=\frac{1}{2}\left(\rho_i+\rho_{i+1}\right)-\frac{1}{2}\frac{\Delta t}{\Delta x}a_{i+\frac{1}{2}}\left(\rho_{i+1}-\rho_i\right) \label{kt1} \end{equation}
\begin{equation} \rho_{i-\frac{1}{2}}=\frac{1}{2}\left(\rho_{i-1}+\rho_{i}\right)-\frac{1}{2}\frac{\Delta t}{\Delta x}a_{i-\frac{1}{2}}\left(\rho_{i}-\rho_{i-1}\right) \label{kt2} \end{equation}

其中$a_{i+\frac{1}{2}}$和$a_{i-\frac{1}{2}}$为定义在网格界面$i+\frac{1}{2}$和$i-\frac{1}{2}$上的局部的波的最大传输速度,其值为通量雅可比矩阵的最大特征值。其可以看作为一个扩散系数为$\frac{1}{2}a\frac{\Delta t}{\Delta x}$的扩散项。KT格式定义库郎数和波的传播速度有关:

\begin{equation} \mathrm{Co}=\frac{\Delta t}{\Delta x}a \label{Co} \end{equation}

因此在库郎数$\mathrm{Co}<1$的情况下,KT格式的数值耗散项系数$\frac{1}{2}a\frac{\Delta t}{\Delta x}$要小于Lax-Friedrichs方法的系数$\frac{1}{2}$。

密度基求解器

首先有连续性方程

\begin{equation} \frac{\partial \rho}{\partial t}+ \nabla\cdot(\rho\bfU)=0 \label{cont} \end{equation}

对其进行有限体积法离散有

\begin{equation} (\rho^{n+1}-\rho^n)\frac{\Delta V}{\Delta t}+ \sum\left(\bfU_f\cdot\bfS_f\rho_f\right)= (\rho^{n+1}-\rho^n)\frac{\Delta V}{\Delta t}+ \sum\left(\phi_f\rho_f\right)=0 \label{cont2} \end{equation}

其中$\phi_f$表示定义在网格单元面上的体积通量,参考方程\eqref{kt1}和\eqref{kt2}有

\begin{equation} \phi_f\rho_f=\frac{1}{2}(\phi_f^+\rho_f^++\phi_f^-\rho_f^-)-\frac{1}{2}a_\max\left(\rho_f^--\rho_f^+\right)=\left(\frac{1}{2}\phi_f^+\rho_f^+-\frac{1}{2}a_\max\rho_f^-\right)+\left(\frac{1}{2}\phi_f^-\rho_f^-+\frac{1}{2}a_\max\rho_f^+ \right) \label{cont3} \end{equation}
\begin{equation} \phi_f^+=\bfU_f^+\cdot\bfS_f,\phi_f^-=\bfU_f^-\cdot\bfS_f \label{phi} \end{equation}

其中$\phi_f^+\rho_f^+,\rho_f^+$表示面$f$对应的owner网格单元向面的插值,$\phi_f^-\rho_f^-,\rho_f^-$表示面$f$对应的neighbor网格单元向面的插值,$a_\max$表示一种人工添加的耗散系数。$a_\max$的大小决定格式的耗散性。KT格式认为对于欧拉方程,不仅需要考虑流动依靠速度传播的特性,还需要考虑波的传播,因此$a_\max$和波的传播速度有关。对于欧拉方程,有特征值为

\begin{equation} \bfU\cdot\bfS,\bfU\cdot\bfS,\bfU\cdot\bfS,\bfU\cdot\bfS+a|\bfS|,\bfU\cdot\bfS-a|\bfS| \label{eigen} \end{equation}

其中$a$表示声速:

\begin{equation} a=\sqrt{\frac{C_p}{C_v}\frac{1}{\psi}} \label{a} \end{equation}

其中$C_p/C_v$为比热容,$\psi$为可压缩性。对于定义在网格单元面$f$上方程\eqref{eigen}中的每个特征值,owner和neighbor网格单元的特征值并不相同(由于间断的存在)。$\bfU\cdot\bfS+a|\bfS|$在面$f$上左右两侧最大特征值的绝对值为

\begin{equation} a_p=\mag\left(\max \left(\bfU_f^+\cdot\bfS_f+a_f^+|\bfS|,\bfU_f^-\cdot\bfS_f+a_f^-|\bfS|,0.0\right)\right) \label{ap} \end{equation}

$\bfU\cdot\bfS-a|\bfS|$在面$f$上左右两侧最大特征值的绝对值为

\begin{equation} a_m=\mag\left(\min \left(\bfU_f^+\cdot\bfS_f-\left(a_f\cdot|\bfS|\right)^+,\bfU_f^-\cdot\bfS_f-\left(a_f\cdot|\bfS|\right)^-,0.0\right)\right) \label{am} \end{equation}

因此,面$f$上的最大的特征值绝对值为

\begin{equation} a_\max=\max(a_p,a_m) \label{amaxSf} \end{equation}

将$a_\max$代入到方程\eqref{cont2}即为使用KT格式离散的连续性方程。对其求解有密度。下面考虑动量方程:

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

动量方程的第一项可以理解为对未知量$\rho\bfU$对时间的导数,第二项可以理解为$\rho\bfU$以速度$\bfU$的传输,为方便代码理解可以改写为

\begin{equation} \frac{\partial \rho_\mathbf{U}}{\partial t}+\nabla \cdot (\rho_\mathbf{U}\mathbf{U})+\nabla p=0 \label{mom2} \end{equation}

对其进行有限体积法离散:

\begin{equation} (\rho_\bfU^{n+1}-\rho_\bfU^n)\frac{\Delta V}{\Delta t}+ \underset{\mathrm{phiUp}}{\underbrace{\sum\left(\phi_f\rho_f\bfU_f\right)+\sum\left(\bfS_fp_f\right)}}=0 \label{mom3} \end{equation}

其中第三项可以分解为

\begin{equation} \bfS_f p_f=\frac{1}{2}\bfS_f\left(p_f^++p_f^-\right) \label{gradp} \end{equation}

将其代入到方程\eqref{mom3}中有

\begin{equation} (\rho_\bfU^{n+1}-\rho_\bfU^n)\frac{\Delta V}{\Delta t}+ \sum\left(\phi_f\rho_f\bfU_f+\frac{1}{2}\left(p_f^++p_f^-\right)\bfS_f\right)=0 \label{mom4} \end{equation}

考虑能量方程

\begin{equation} \frac{\p \rho E}{\p t}+\nabla\cdot(\bfU \rho E)+\nabla\cdot(\bfU p)=0 \label{en} \end{equation}

能量方程的第一项可以理解为未知量$\rho E$对时间的导数,第二项可以理解为$\rho E$以速度$\bfU$的传输,为方便代码理解可以写为

\begin{equation} \frac{\p \rho_E}{\p t}+\underset{\mathrm{phiEp}}{\underbrace{\nabla\cdot(\bfU \rho_E+\bfU p)}}=0 \label{en2} \end{equation}

对其离散和动量方程大同小异。

密度基求解器和压力基求解器的不同在于分离式密度基求解器顺序求解密度、速度、能量方程,然后依据压力和密度的代数方程更新压力即可,不需要压力基求解器的迭代过程。

代码分析

int main(int argc, char *argv[])
{
    // --- 略去非主要代码

    while (runTime.run())
    {
        surfaceScalarField rho_pos(interpolate(rho, pos)); //owner网格的$\rho_f^+$
        surfaceScalarField rho_neg(interpolate(rho, neg)); //neighbor网格的$\rho_f^-$

        surfaceVectorField rhoU_pos(interpolate(rhoU, pos, U.name())); //owner网格的$\rho_{\bfU,f}^+$
        surfaceVectorField rhoU_neg(interpolate(rhoU, neg, U.name())); //owner网格的$\rho_{\bfU,f}^-$

        volScalarField rPsi("rPsi", 1.0/psi);
        surfaceScalarField rPsi_pos(interpolate(rPsi, pos, T.name())); //owner网格的$\frac{1}{\psi_f^+}$
        surfaceScalarField rPsi_neg(interpolate(rPsi, neg, T.name())); //neighbor网格的$\frac{1}{\psi_f^-}$

        surfaceScalarField e_pos(interpolate(e, pos, T.name())); //owner网格的$E^+$
        surfaceScalarField e_neg(interpolate(e, neg, T.name())); //neighbor网格的$E^-$

        surfaceVectorField U_pos("U_pos", rhoU_pos/rho_pos); //owner网格的$\bfU_f^+$
        surfaceVectorField U_neg("U_neg", rhoU_neg/rho_neg); //neighbor网格的$\bfU_f^-$

        surfaceScalarField p_pos("p_pos", rho_pos*rPsi_pos); //owner网格的$p_f^+$
        surfaceScalarField p_neg("p_neg", rho_neg*rPsi_neg); //neighbor网格的$p_f^-$

        surfaceScalarField phiv_pos("phiv_pos", U_pos & mesh.Sf()); //owner网格的$\phi_f^+$
        surfaceScalarField phiv_neg("phiv_neg", U_neg & mesh.Sf()); //neighbor网格的$\phi_f^-$

        volScalarField c("c", sqrt(thermo.Cp()/thermo.Cv()*rPsi)); //音速$a$
		
        //owner网格的$a_f^+|\bfS_f|$
        surfaceScalarField cSf_pos
        (
            "cSf_pos",
            interpolate(c, pos, T.name())*mesh.magSf()
        ); 
		
        //neighbor网格的$a_f^-|\bfS_f|$
        surfaceScalarField cSf_neg
        (
            "cSf_neg",
            interpolate(c, neg, T.name())*mesh.magSf()
        );

        //方程\eqref{ap}中的$a_p$
        surfaceScalarField ap
        (
            "ap",
            max(max(phiv_pos + cSf_pos, phiv_neg + cSf_neg), v_zero)
        );

        //方程\eqref{am}中的$a_m$
        surfaceScalarField am
        (
            "am",
            min(min(phiv_pos - cSf_pos, phiv_neg - cSf_neg), v_zero)
        );

        surfaceScalarField a_pos("a_pos", ap/(ap - am));

        //方程\eqref{amaxSf}中的$a_\max$
        surfaceScalarField amaxSf("amaxSf", max(mag(am), mag(ap)));

        surfaceScalarField aSf("aSf", am*a_pos);

        if (fluxScheme == "Tadmor")
        {
            aSf = -0.5*amaxSf;
            a_pos = 0.5;
        }

        surfaceScalarField a_neg("a_neg", 1.0 - a_pos);

        phiv_pos *= a_pos;
        phiv_neg *= a_neg;

        //方程\eqref{cont3}右侧括号内第一项
        surfaceScalarField aphiv_pos("aphiv_pos", phiv_pos - aSf);
	
        //方程\eqref{cont3}右侧括号内第二项
        surfaceScalarField aphiv_neg("aphiv_neg", phiv_neg + aSf);

        // Reuse amaxSf for the maximum positive and negative fluxes
        // estimated by the central scheme
        amaxSf = max(mag(aphiv_pos), mag(aphiv_neg));

        //基于音速计算库郎数
        #include "centralCourantNo.H"
        #include "readTimeControls.H"

        if (LTS)
        {
            #include "setRDeltaT.H"
        }
        else
        {
            #include "setDeltaT.H"
        }

        runTime++;

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

        //方程\eqref{cont3}中的$\phi_f\rho_f$
        phi = aphiv_pos*rho_pos + aphiv_neg*rho_neg;

        //方程\eqref{mom4}中的$\left(\phi_f\rho_f\bfU_f+\frac{1}{2}\left(p_f^++p_f^-\right)\bfS_f\right)$
        surfaceVectorField phiUp
        (
            (aphiv_pos*rhoU_pos + aphiv_neg*rhoU_neg)
          + (a_pos*p_pos + a_neg*p_neg)*mesh.Sf()
        );

        surfaceScalarField phiEp
        (
            "phiEp",
            aphiv_pos*(rho_pos*(e_pos + 0.5*magSqr(U_pos)) + p_pos)
          + aphiv_neg*(rho_neg*(e_neg + 0.5*magSqr(U_neg)) + p_neg)
          + aSf*p_pos - aSf*p_neg
        );

        volScalarField muEff("muEff", turbulence->muEff());
        volTensorField tauMC("tauMC", muEff*dev2(Foam::T(fvc::grad(U))));

        // --- 求解连续性方程获得密度
        solve(fvm::ddt(rho) + fvc::div(phi));

        // --- 求解动量方程获得$\rho_\bfU$
        solve(fvm::ddt(rhoU) + fvc::div(phiUp));

        U.ref() =
            rhoU()
           /rho();
        U.correctBoundaryConditions();
        rhoU.boundaryFieldRef() == rho.boundaryField()*U.boundaryField();

        if (!inviscid)
        {
            solve
            (
                fvm::ddt(rho, U) - fvc::ddt(rho, U)
              - fvm::laplacian(muEff, U)
              - fvc::div(tauMC)
            );
            rhoU = rho*U;
        }

        surfaceScalarField sigmaDotU
        (
            "sigmaDotU",
            (
                fvc::interpolate(muEff)*mesh.magSf()*fvc::snGrad(U)
              + fvc::dotInterpolate(mesh.Sf(), tauMC)
            )
          & (a_pos*U_pos + a_neg*U_neg)
        );

        // --- 求解能量方程获得能量
        solve
        (
            fvm::ddt(rhoE)
          + fvc::div(phiEp)
          - fvc::div(sigmaDotU)
        );

        e = rhoE/rho - 0.5*magSqr(U);
        e.correctBoundaryConditions();
        thermo.correct();
        rhoE.boundaryFieldRef() ==
            rho.boundaryField()*
            (
                e.boundaryField() + 0.5*magSqr(U.boundaryField())
            );

        if (!inviscid)
        {
            solve
            (
                fvm::ddt(rho, e) - fvc::ddt(rho, e)
              - fvm::laplacian(turbulence->alphaEff(), e)
            );
            thermo.correct();
            rhoE = rho*(e + 0.5*magSqr(U));
        }

        // --- 依据密度和可压缩性的代数方程更新压力
        p.ref() =
            rho()
           /psi();
        p.correctBoundaryConditions();
        rho.boundaryFieldRef() == psi.boundaryField()*p.boundaryField();

        turbulence->correct();

        runTime.write();

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

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

    return 0;
}

参考文献

Greenshields, C. J., Weller, H. G., Gasparini, L., Reese, J. M., 2010. Implementation of semi-discrete, non-staggered central schemes in a colocated, polyhedral, finite volume framework, for high-speed viscous flows. International journal for numerical methods in fluids, 63, 1-21.

Kurganov, A., Tadmor, E., 2000. New high-resolution central schemes for nonlinear conservation laws and convection–diffusion equations. Journal of Computational Physics, 160, 241-282.

更新历史
2018.02.23创立页面

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