rhoCentralFoam解析


李东岳


1. 引言

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

2. 控制方程与算法

考虑速度已知的对流方程

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

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

\begin{equation} \int_{\Delta V}\int_\dt\frac{\partial \rho}{\partial t}\rd t\rd V +\int_{\Delta V}\int_\dt \nabla\cdot (\rho\bfU)\rd t\rd V=0 \end{equation}

\begin{equation} (\rho_\rP^{t+\dt}-\rho_\rP)\Delta V_\rP + \sum \bfU_f\cdot\bfS_f\rho_f\Delta t=0 \label{adv2} \end{equation}

其中下标$_\rP$表示当前网格点,$_f$表示网格单元面上的值。整理有:

\begin{equation} \rho_\rP^{t+\dt}= \rho_\rP - \frac{\dt}{\Delta V_\rP}\sum \phi_f\rho_f \label{adv3} \end{equation}

其中$\phi_f=\bfU_f\cdot\bfS_f$。当前时间步的$\rho_\rP$已知,若有网格单元面上的$\rho_f$,则可求下一个时间步的结果。如果认为网格单元面上的值为网格面左右点的值的平均,即:

\begin{equation} \phi_f\rho_{f}=\frac{1}{2}\left(\phi_\own^f\rho_\own^f+\phi_\nei^f\rho_\nei^f\right) \label{c1} \end{equation}

其中$\phi_\own^f$表示从own网格对面$f$进行高分辨率插值。将方程\eqref{c1}代入到\eqref{adv3}会发现其是无条件不稳定的。一阶Lax-Friedrichs方法在无条件不稳定格式的基础上,将其改写为

\begin{equation} \begin{split} &\phi_f\rho_{f}=\frac{1}{2}\left(\phi_\own^f\rho_\own^f+\phi_\nei^f\rho_\nei^f\right) -\frac{1}{2}\frac{|\mathbf{d}||\bfS_f|}{\Delta t}\left(\rho_\nei^f-\rho_\own^f\right) \\\\ &\phi_\own^f=\phi_\own,\phi_\nei^f=\phi_\nei,\rho_\own^f=\rho_\own,\rho_\nei^f=\rho_\nei \end{split}\label{lax1} \end{equation}

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

\begin{equation} \begin{split} \phi_f\rho_f&=\frac{1}{2}\left(\phi_\own^f\rho_\own^f+\phi_\nei^f\rho_\nei^f\right) -\frac{1}{2}\max(a_{f})|\bfS_f|\left(\rho_\nei^f-\rho_\own^f\right) \\\\ &=\frac{1}{2}\left(\phi_\own^f+\max(a_{f})|\bfS_f|\right)\rho_\own^f +\frac{1}{2}\left(\phi_\nei^f-\max(a_{f})|\bfS_f| \right)\rho_\nei^f \end{split} \label{cont3} \end{equation}

其中$\frac{1}{2}$表示权重(KT格式)。$\max(a_f)$为定义在网格面$f$上的局部的波的最大传输速度,其值为通量雅可比矩阵的特征值绝对值的最大值。同时,可以看出,在局部黎曼扇区域较小时,有:

\begin{equation} \max(a_{f})<\frac{|\bfd|}{\dt} \end{equation}

可见,KT格式的扩散系数小于Lax-Friedrichs方法。KT格式认为对于欧拉方程,不仅需要考虑流动依靠速度传播的特性,还需要考虑波的传播,因此$\max(a_{f})$和波的传播速度有关。对于欧拉方程,其特征值为

\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} \begin{split} \max(a_{f}) = \max( |\bfU_\own^f\cdot\bfS_f+c_\own^f|\bfS_f||, |\bfU_\nei^f\cdot\bfS_f+c_\nei^f|\bfS_f||, |\bfU_\own^f\cdot\bfS_f-c_\own^f|\bfS_f||, |\bfU_\nei^f\cdot\bfS_f-c_\nei^f|\bfS_f||, \end{split} \label{amax} \end{equation}

将$\max(a_{f})$求出后,代入到方程\eqref{cont3}即为使用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$的传输,这样既对守恒变量$\rho_\mathbf{U}$进行离散。为方便代码理解可以改写为

\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 t}{\Delta V}\sum\left(\phi_f{\rho_{\bfU}}_f+\bfS_fp_f\right)=0 \label{mom3} \end{equation}

其中$\sum\left(\phi_f{\rho_{\bfU}}_f+\bfS_fp_f\right)$为$\mathrm{phiUp}$,第三项可以分解为

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

对其求解有$\rho_\bfU$。随后速度进行更新:

\begin{equation} \bfU=\frac{\rho_\bfU}{\rho} \label{rhourho} \end{equation}

参考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}+\nabla\cdot(\rho_E\bfU +p\bfU)-\nabla \cdot(\tau \cdot \mathbf{U})-\nabla\cdot(\alpha_\mathrm{eff}\nabla e)=0 \label{en2} \end{equation}

其中$\nabla\cdot(\rho_E\bfU +p\bfU )$为$\mathrm{phiEp}$。在无粘的情况下,方程\eqref{en2}为

\begin{equation} \frac{\p \rho_E}{\p t}+\nabla\cdot(\rho_E\bfU+p\bfU)=0 \label{inviscod} \end{equation}

对其进行离散有:

\begin{equation} (\rho_E^{n+1}-\rho_E^n)+ \frac{\Delta t}{\Delta V}\sum\left(\phi_f\left({\rho_E}_f+p_f\right)\right)=0 \label{eeqn} \end{equation}

求解方程\eqref{eeqn}有$\rho_E$。在考虑粘度的情况下,方程\eqref{en2}并未封闭,求解器中采用一种操作符分裂的方式进行求解:首先剥离比内能项求解

\begin{equation} \frac{\p \rho_E}{\p t}+\nabla\cdot(\bfU \rho_E+\bfU p)-\nabla \cdot(\tau \cdot \mathbf{U})=0 \label{viscod1} \end{equation}

然后更新比内能:

\begin{equation} \frac{\p e}{\p t}-\frac{\p e}{\p t}-\nabla\cdot(\alpha_\mathrm{eff}\nabla e)=0 \label{viscod2} \end{equation}

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

3. 代码分析
while (runTime.run())
{
    surfaceScalarField rho_pos(interpolate(rho, pos));
    surfaceScalarField rho_neg(interpolate(rho, neg));
    // 调用高分辨率格式组建$\rho_\own^f,\rho_\nei^f$

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

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

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

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

    surfaceScalarField e_pos(interpolate(e, pos, T.name()));
    surfaceScalarField e_neg(interpolate(e, neg, T.name()));
    // 调用高分辨率格式组建$e_\own^f,e_\nei^f$

    surfaceVectorField U_pos("U_pos", rhoU_pos/rho_pos);
    surfaceVectorField U_neg("U_neg", rhoU_neg/rho_neg);
    // 计算$\bfU_\own^f={\rho_\bfU}_\own^f/\rho_\own^f,\bfU_\nei^f={\rho_\bfU}_\nei^f/\rho_\nei^f$

    surfaceScalarField p_pos("p_pos", rho_pos*rPsi_pos);
    surfaceScalarField p_neg("p_neg", rho_neg*rPsi_neg);
    // 计算$p_\own^f=\frac{1}{\psi_\own^f}\rho_\own^f,p_\nei=\frac{1}{\psi_\nei^f}\rho_\nei^f$

    surfaceScalarField phi_pos("phi_pos", U_pos & mesh.Sf());
    surfaceScalarField phi_neg("phi_neg", U_neg & mesh.Sf());
    // 计算$\phi_\own^f=\bfU_\own^f\cdot\bfS_f,\phi_\nei^f=\bfU_\nei^f\cdot\bfS_f$

    volScalarField c("c", sqrt(thermo.Cp()/thermo.Cv()*rPsi));
    // 音速
    
    surfaceScalarField cSf_pos
    (
      "cSf_pos",
      interpolate(c, pos, T.name())*mesh.magSf()
    );// 计算$c_\own^f|\bfS_f|$
    
    surfaceScalarField cSf_neg
    (
      "cSf_neg",
      interpolate(c, neg, T.name())*mesh.magSf()
    );
    // // 计算$c_\nei^f|\bfS_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;
    // 计算$\phi_\own^f=0.5\bfU_\own^f\cdot\bfS_f,\phi_\nei^f=0.5\bfU_\nei^f\cdot\bfS_f$


    surfaceScalarField aphi_pos("aphi_pos", phi_pos + 0.5*amaxSf);
    // 计算$0.5\left(\phi_\own^f + \max(a_f)|\bfS_f|\right)$
    surfaceScalarField aphi_neg("aphi_neg", phi_neg - 0.5*amaxSf);
    // 计算$0.5\left(\phi_\nei^f - \max(a_f)|\bfS_f|\right)$
    

    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;

    // --- 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中,自带的算例主要为以下几个:

上述算例已经在OpenFOAM包含。在此不做过多介绍。

4.1. Shu-Osher问题

本文中的算例研究Shu-Osher问题,其为一个经典的可压缩流算法验证问题。Shu-Osher问题主要用于检测离散格式的分辨率,研究离散格式是否能分辨小尺度的流动细节。可用于标定数值粘度的大小。研究Shu-Osher问题主要采用一维欧拉方程。问题中,马赫数为3的激波向前步进到附加周期扰动的密度场中,随着激波的步进,激波后会出现俩个频率不同的密度场扰动。高频率的密度场扰动周期为低频率周期的1/2,扰动互相连接。算例网格为采用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文件夹的物性主要设置如下:

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


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

点击下载

访问密码: x9zc9r

4.2. 方腔超音速流动

算例描述,待更新


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

东岳流体 2014 - 2020
勘误、讨论、补充内容请前往CFD中文网