rhoCentralFoam解析


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

1. 引言

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

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

如果打不开图像,请右键在新标签页打开图像后刷新几次
2. 控制方程与算法
2.1. Lax-Friedrichs格式

考虑一维常系数对流方程

\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_\own^{n+1}-\rho_\own)\Delta x\Delta y\Delta z+\left(u_{f,\rR}\rho_{f,\rR}-u_{f,\rL}\rho_{f,\rL}\right)\Delta t\Delta y\Delta z=0 \label{adv2} \end{equation}

其中$\own$表示当前网格体心,$f,\rR$和$f,\rL$表示网格单元左右面。令

\begin{equation} \rho\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_\own^{n+1}-m_\own)+\left(u_{f,\rR}\rho_{f,\rR}-u_{f,\rL}\rho_{f,\rL}\right)\Delta t|\bfS_f|=0 \end{equation}
\begin{equation} m_\own^{n+1}=m_\own-\Delta t\left(\phi_{f,\rR}\rho_{f,\rR}-\phi_{f,\rL}\rho_{f,\rL}\right) \label{adv3} \end{equation}
\begin{equation} \phi_{f,\rR}=u_{f,\rR}|\bfS_f|,\phi_{f,\rL}=u_{f,\rL}|\bfS_f| \label{phi1} \end{equation}

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

\begin{equation} \phi_f\rho_{f}=\frac{1}{2}\left(\phi^+\rho^++\phi^-\rho^-\right) \label{c1} \end{equation}

其中$\rho^+,\rho^-$表示从网格面$f$的正向、反向对其进行插值。将方程\eqref{c1}代入到\eqref{adv3}会发现其是无条件不稳定的。Lax-Friedrichs方法在无条件不稳定格式的基础上,将其改写为

\begin{equation} \phi_{f,\rR}\rho_{f,\rR}=\frac{1}{2}\left(\phi_{f,\rR}^+\rho_{f,\rR}^++\phi_{f,\rR}^-\rho_{f,\rR}^-\right) -\frac{1}{2}\frac{\Delta x}{\Delta t}|\bfS_f|\left(\rho_{f,\rR}^+-\rho_{f,\rR}^-\right) \label{lax1} \end{equation}
\begin{equation} \phi_{f,\rL}\rho_{f,\rL}=\frac{1}{2}\left(\phi_{f,\rL}^+\rho_{f,\rL}^++\phi_{f,\rL}^-\rho_{f,\rL}^-\right) -\frac{1}{2}\frac{\Delta x}{\Delta t}|\bfS_f|\left(\rho_{f,\rL}^+-\rho_{f,\rL}^-\right) \label{lax2} \end{equation}

注意有$\frac{\Delta x}{\Delta t}|\bfS_f|=\frac{\Delta V}{\Delta t}$。方程\eqref{lax1}和\eqref{lax2}右边的后两项可以认为是扩散系数为$\frac{1}{2} \frac{\Delta x}{\Delta t}\Delta x$的数值扩散项。正是因为这俩项,才使得Lax-Friedrichs方法稳定。且在网格单元越小的时候,数值耗散越小

2.2. 中心/中心迎风格式

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

\begin{equation} \phi_{f,\rR}\rho_{f,\rR}=\frac{1}{2}\left(\phi_{f,\rR}^+\rho_{f,\rR}^++\phi_{f,\rR}^-\rho_{f,\rR}^-\right) -\frac{1}{2}\max(a_{f,\rR})|\bfS_f|\left(\rho_{f,\rR}^+-\rho_{f,\rR}^-\right) \label{kt1} \end{equation}
\begin{equation} \phi_{f,\rL}\rho_{f,\rL}=\frac{1}{2}\left(\phi_{f,\rL}^+\rho_{f,\rL}^++\phi_{f,\rL}^-\rho_{f,\rL}^-\right) -\frac{1}{2}\max(a_{f,\rL})|\bfS_f|\left(\rho_{f,\rL}^+-\rho_{f,\rL}^-\right) \label{kt2} \end{equation}

其中$\frac{1}{2}$表示权重(KT格式)。$\max(a_{f,\rR})$和$\max(a_{f,\rL})$为定义在网格面$f,\rR$和$f,\rL$上的局部的波的最大传输速度,其值为通量雅可比矩阵的特征值绝对值的最大值。它表示一种人工添加的耗散系数。值的大小决定格式的耗散性。KT格式认为对于欧拉方程,不仅需要考虑流动依靠速度传播的特性,还需要考虑波的传播,因此$\max(a_{f})$和波的传播速度有关。方程\eqref{kt1}与\eqref{kt2}右侧的后一项可以看作为一个扩散系数为$\frac{1}{2}\max(a_{f})\Delta x$的扩散项。KT格式定义库郎数和波的传播速度有关:

\begin{equation} \mathrm{Co}=\frac{\Delta t}{\Delta x}\max(a_{f}) \label{Co} \end{equation}

因此在库郎数$\mathrm{Co}<1$的情况下,中心格式的数值耗散项系数$\max(a_{f})<\frac{\Delta x}{\Delta t}$,因此中心迎风格式的扩散系数$\frac{1}{2}\max(a_{f})\Delta x$要小于Lax-Friedrichs方法的系数$\frac{1}{2} \frac{\Delta x}{\Delta t}\Delta x$。

2.3. 密度基求解器

首先有连续性方程

\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(\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}\max(a_{f})|\bfS_f|\left(\rho_f^+-\rho_f^-\right) \\\\ =\frac{1}{2}\left(\phi_f^+-\max(a_{f})|\bfS_f|\right)\rho_f^++\frac{1}{2}\left(\phi_f^-+\max(a_{f})|\bfS_f| \right)\rho_f^- \label{cont3} \end{equation}
\begin{equation} \phi_f^+=\bfU_f^+\cdot\bfS_f,\phi_f^-=\bfU_f^-\cdot\bfS_f \label{phi} \end{equation}

欧拉方程的特征值为

\begin{equation} \bfU_f\cdot\bfS_f,\bfU_f\cdot\bfS_f,\bfU_f\cdot\bfS_f,\bfU_f\cdot\bfS_f+c_f|\bfS_f|,\bfU_f\cdot\bfS_f-c_f|\bfS_f| \label{eigen} \end{equation}
其中$c$表示声速:

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

其中$C_p/C_v$为比热容,$\psi$为可压缩性。在这些特征值中,最大和最小特征值必然在$\bfU_f\cdot\bfS_f+c_f|\bfS_f|$或$\bfU_f\cdot\bfS_f-c_f|\bfS_f|$中产生。这两项的特征值的绝对值的最大值为:

\begin{equation} \max(a_{f})=\max(|\bfU_f^+\cdot\bfS_f+c_f^+|\bfS_f||,|\bfU_f^-\cdot\bfS_f+c_f^-|\bfS_f||, \\\\ |\bfU_f^+\cdot\bfS_f-c_f^+|\bfS_f|,|\bfU_f^-\cdot\bfS_f-c_f^-|\bfS_f|) \label{amax} \end{equation}
将$\max(a_{f})$求出后,代入到方程\eqref{cont3}与\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_{\bfU}}_f+\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}
对其求解有$\rho_\bfU$。随后速度进行更新:
\begin{equation} \bfU=\frac{\rho_\bfU}{\rho} \label{rhourho} \end{equation}

2.4. 能量方程

参考CFD中的能量方程rhoCentralFoam考虑的为比能方程:

\begin{equation} \frac{\p \rho E}{\p t}+\nabla\cdot(\bfU \rho E)-\alpha_\mathrm{eff}\nabla e= -\nabla\cdot(\bfU p)+\nabla \cdot(\tau \cdot \mathbf{U}) \label{en} \end{equation}

比能方程的第一项可以理解为未知量$\rho E$对时间的导数,第二项可以理解为$\rho E$以速度$\bfU$的传输,第三项为热传导。右侧第一项表示正应力做功,第二项表示粘性力做功。为方便代码理解可以写为

\begin{equation} \frac{\p \rho_E}{\p t}+\underset{\mathrm{phiEp}}{\underbrace{\nabla\cdot(\rho_E\bfU +\bfU p)}}-\nabla \cdot(\tau \cdot \mathbf{U})-\alpha_\mathrm{eff}\nabla e=0 \label{en2} \end{equation}

在无粘的情况下,方程\eqref{en2}为

\begin{equation} \frac{\p \rho_E}{\p t}+\underset{\mathrm{phiEp}}{\underbrace{\nabla\cdot(\rho_E\bfU+\bfU p)}}=0 \label{inviscod} \end{equation}
对其进行离散有:
\begin{equation} (\rho_E^{n+1}-\rho_E^n)\frac{\Delta V}{\Delta t}+ \underset{\mathrm{phiEp}}{\underbrace{\sum\left(\phi_f\left({\rho_E}_f+p_f\right)\right)}}=0 \label{eeqn} \end{equation}
求解方程\eqref{eeqn}有$\rho_E$。无粘流的密度方程、动量方程与能量方程即欧拉方程,为一个严格的双曲系统,特征变量以不同的特征速度进行传播。积分形式的欧拉方程允许求解出现间断。若结合高分辨率格式间断更加明晰。在考虑粘度的情况下,方程\eqref{en2}并未封闭,求解器中采用一种操作符分裂的方式进行求解:

在考虑粘度的情况下,方程变为抛物特征。间断将变得光顺。在求解的时候,密度基求解器和压力基求解器的不同在于分离式密度基求解器顺序求解密度、速度、能量方程,然后依据压力和密度的代数方程更新压力即可,不需要压力基求解器的迭代过程。

3. 代码分析
while (runTime.run())
    {
        surfaceScalarField rho_pos(interpolate(rho, pos));
        surfaceScalarField rho_neg(interpolate(rho, neg));
        // 组建$\rho_f^+,\rho_f^-$

        surfaceVectorField rhoU_pos(interpolate(rhoU, pos, U.name()));
        surfaceVectorField rhoU_neg(interpolate(rhoU, neg, U.name()));
        // 调用高分辨率格式组建${\rho_\bfU}_f^+,{\rho_\bfU}_f^-$

        surfaceScalarField rhoE_pos(interpolate(rhoE, pos, T.name()));
        surfaceScalarField rhoE_neg(interpolate(rhoE, neg, T.name()));
        // 调用高分辨率格式组建${\rho_E}_f^+,{\rho_E}_f^-$

        volScalarField rPsi("rPsi", 1.0/psi);

        surfaceScalarField rPsi_pos(interpolate(rPsi, pos, T.name()));
        surfaceScalarField rPsi_neg(interpolate(rPsi, neg, T.name()));
        // 调用高分辨率格式组建$\psi_f^+,\psi_f^-$

        surfaceScalarField e_pos(interpolate(e, pos, T.name()));
        surfaceScalarField e_neg(interpolate(e, neg, T.name()));
        // 调用高分辨率格式组建$e_f^+,e_f^-$
        
        surfaceVectorField U_pos("U_pos", rhoU_pos/rho_pos);
        surfaceVectorField U_neg("U_neg", rhoU_neg/rho_neg);
        // 计算$\bfU_f^+,\bfU_f^-$,其为高分辨率格式

        surfaceScalarField p_pos("p_pos", rho_pos*rPsi_pos);
        surfaceScalarField p_neg("p_neg", rho_neg*rPsi_neg);
        // 计算$p_f^+,p_f^-$

        surfaceScalarField phi_pos("phi_pos", U_pos & mesh.Sf());
        surfaceScalarField phi_neg("phi_neg", U_neg & mesh.Sf());
        // 计算$\phi_f^+,\phi_f^-$

        volScalarField c("c", sqrt(thermo.Cp()/thermo.Cv()*rPsi));
        // 音速
        surfaceScalarField cSf_pos
        (
            "cSf_pos",
            interpolate(c, pos, T.name())*mesh.magSf()
        );
        surfaceScalarField cSf_neg
        (
            "cSf_neg",
            interpolate(c, neg, T.name())*mesh.magSf()
        );
        // 调用高分辨率格式组建$c_f^+,c_f^-$

        surfaceScalarField amaxSf
        (
            "amaxSf", 
            max
            (
                max(mag(phi_pos + cSf_pos), mag(phi_neg + cSf_neg)),
                max(mag(phi_pos - cSf_pos), mag(phi_neg - cSf_neg))
            )
        );
        // 方程\eqref{amax}

        phi_pos *= 0.5;
        phi_neg *= 0.5;

        surfaceScalarField aphi_pos("aphi_pos", phi_pos + 0.5*amaxSf);
        surfaceScalarField aphi_neg("aphi_neg", phi_neg - 0.5*amaxSf);
		
        phi = aphi_pos*rho_pos + aphi_neg*rho_neg;
        // 方程\eqref{cont3}

        surfaceVectorField phiUp
        (
            (aphi_pos*rhoU_pos + aphi_neg*rhoU_neg)
          + (0.5*p_pos + 0.5*p_neg)*mesh.Sf()
        );
        // 方程\eqref{mom3}中的phiUp

        surfaceScalarField phiEp
        (
            "phiEp",
            aphi_pos*(rho_pos*(e_pos + 0.5*magSqr(U_pos)) + p_pos)
          + aphi_neg*(rho_neg*(e_neg + 0.5*magSqr(U_neg)) + p_neg)
          - 0.5*amaxSf*p_pos + 0.5*amaxSf*p_neg
        );
        // 方程\eqref{eeqn}中的phiEp

        //- For Courant number
        amaxSf = max(mag(aphi_pos), mag(aphi_neg));

        #include "centralCourantNo.H"
        #include "readTimeControls.H"
        #include "setDeltaT.H"
        runTime++;
        Info<< "Time = " << runTime.timeName() << nl << endl;

        phi = aphi_pos*rho_pos + aphi_neg*rho_neg;

        surfaceVectorField phiUp
        (
            (aphi_pos*rhoU_pos + aphi_neg*rhoU_neg)
          + (0.5*p_pos + 0.5*p_neg)*mesh.Sf()
        );

        surfaceScalarField phiEp
        (
            "phiEp",
            aphi_pos*(rho_pos*(e_pos + 0.5*magSqr(U_pos)) + p_pos)
          + aphi_neg*(rho_neg*(e_neg + 0.5*magSqr(U_neg)) + p_neg)
          - 0.5*amaxSf*p_pos + 0.5*amaxSf*p_neg
        );
        // 保证有界性

        // --- Solve density
        solve(fvm::ddt(rho) + fvc::div(phi));

        // --- Solve momentum
        solve(fvm::ddt(rhoU) + fvc::div(phiUp));

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

        // --- Solve energy
        solve
        (
            fvm::ddt(rhoE)
          + fvc::div(phiEp)
        );

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

        p.ref() = rho()/psi();
        p.correctBoundaryConditions();
        rho.boundaryFieldRef() == psi.boundaryField()*p.boundaryField();

        runTime.write();

        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }
4. 验证算例

rhoCentralFoam自带若干算例。在OpenFOAM-7中,自带的算例主要为以下几个:

上述算例在这里可以下载。在此不做过多介绍。本文中的算例研究Shu-Osher问题,对标相关文献的研究工作 (Johnsen et al., 2010)。Shu-Osher问题主要用于检测离散格式的分辨率,研究离散格式是否能分辨小尺度的流动细节。可用于标定数值粘度的大小。研究Shu-Osher问题主要采用一维欧拉方程。问题中,马赫数为3的激波向前步进到附加周期扰动的密度场中,随着激波的步进,激波后会出现俩个频率不同的密度场扰动。高频率的密度场扰动周期为低频率周期的1/2,扰动互相连接。

4.1. 模型设置

算例网格为采用blockMesh生成的1D纯六面体结构网格。前后两侧均为固定值边界条件,上下左右为空类型(1D模拟)。算例运行1.8s后,可与文献对比密度场结果。

\begin{equation}\label{shu} (\rho,u,p)=\left\{\begin{matrix} 3.857143,2.629369,10.33333, & x<-4\\ 1+0.2\mathrm{sin}(5x),0,1, & x\geq -4 \end{matrix}\right. \end{equation}

constant文件夹的物性主要设置如下:

4.2. 算例运行

直接运行./Allrun即可。通过调用方程\eqref{Co}中不同的$a$值(需要在求解器中进行代码修改),研究数值粘度对结果的影响(如下图)。


采用不同的$a$值计算的Shu-Osher问题。黑点:文献参考值

点击下载

访问密码: x9zc9r

参考文献

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.

Johnsen, E., Johan L., Ankit V. B., William H. C., Parviz M., Britton J. O., Pradeep S. R. et al. 2010. Assessment of high-resolution methods for numerical simulations of compressible turbulence with shock waves. Journal of Computational Physics 229, 1213-1237.

更新历史
2019.8.23增加验证算例 | 2019.6.30大修能量方程部分

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