版本对应：

OpenFOAM新版的shockFluid模块，对应其他版本的rhoCentralFoam。

CFD: 可压 + 密度基

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

三维欧拉方程的最大波速

双曲系统存在不同的特征值。对于本文讨论的双曲欧拉方程，其存在5个PDE方程，同时其特征值为下面5个传输标量（前3个为相同的特征值）：

(1)\[ \bfU_f\cdot\bfn_f,\bfU_f\cdot\bfn_f,\bfU_f\cdot\bfn_f,\bfU_f\cdot\bfn_f+c_f,\bfU_f\cdot\bfn_f-c_f \label{eigen} \]

其中\(c\)表示声速：

(2)\[ c=\sqrt{\frac{C_p}{C_v}\frac{1}{\psi}} \]

其中\(C_p,C_v\)表示定压/定容比热容，\(\psi\)表示可压缩性。如果考虑网格面的大小\(|\bfS_f|\)，则欧拉方程存在下面5个传输通量（前3个为相同的通量）：

(3)\[ \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| \]

其可以理解为，针对3D欧拉方程，其传输通量除了\(\bfU_f\cdot\bfS_f\)，还存在其他两个传输通量\(\bfU_f\cdot\bfS_f+c_f|\bfS_f|\)以及\(\bfU_f\cdot\bfS_f-c_f|\bfS_f|\)。

注意：

\(c^f\)一定为一个正数。

在这些通量中，最大和最小传输通量必然在\(\bfU_f\cdot\bfS_f+c_f|\bfS_f|\)或\(\bfU_f\cdot\bfS_f-c_f|\bfS_f|\)中产生。结合网格面的两个方向性，因此传输通量的最大值与最小值必然在下面4项中产生：

(4)\[\begin{split} \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}\]

\(_\own\)表示面\(f\)的own网格单元，\(_\nei\)表示面\(f\)的nei网格单元。因为\(c^f\)一定为一个正数，那么上述方程第一行一定会产生一个最大值，下面一行会产生一个数学上的最小值。定义：

(5)\[\begin{split} ap_f=\max (\bfU_\own^f\cdot\bfS_f+c_\own^f|\bfS_f|, \bfU_\nei^f\cdot\bfS_f+c_\nei^f|\bfS_f|) \\\\ am_f=\min (\bfU_\own^f\cdot\bfS_f-c_\own^f|\bfS_f|,\bfU_\nei^f\cdot\bfS_f-c_\nei^f|\bfS_f|) \end{split}\]

进一步如果考虑\(ap_f,am_f\)的正负，那么\(\max(a_f)\)可以定义为：

(6)\[ \max(a_f) = \max (|ap_f|,|am_f|) \]

至此，三维欧拉方程的最大波速的传输速度为\(\max(a_f)\)。其为一个正值。

控制方程与算法

首先有密度方程：

(7)\[ \frac{\partial \rho}{\partial t}+ \nabla\cdot (\rho\bfU)=0 \]

其中\(\rho\)表示密度，\(\bfU\)表示速度。对方程(7)进行显性离散有：

(8)\[ \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 \]

(9)\[ (\rho_\rP^{t+\dt}-\rho_\rP)V_\rP + \sum \bfU_f\cdot\bfS_f\rho_f\Delta t=0 \]

其中下标\(_\rP\)表示当前网格点，\(_f\)表示网格单元面上的值。整理有：

(10)\[ \rho_\rP^{t+\dt}= \rho_\rP - \frac{\dt}{V_\rP}\sum \phi_f\rho_f \]

其中\(\phi_f=\bfU_f\cdot\bfS_f\)。如果当前时间步的\(\rho_\rP\)已知，同时有网格单元面上的\(\rho_f\)，则可求下一个时间步的结果。假定认为网格单元面上的值为网格面左右点的值的算数平均（同时假定网格是均匀的），即：

(11)\[ \phi_f\rho_{f}=\frac{1}{2}\left(\phi_\own^f\rho_\own^f+\phi_\nei^f\rho_\nei^f\right) \]

注意：

\(\phi_\own^f\)只有在一阶格式的情况下等于\(\phi_\own\)。如果使用高阶格式，\(\phi_\own^f\)不等于\(\phi_\own\)。\(\phi_\nei^f\)同理。

其中\(\phi_\own^f\)表示从own网格对面\(f\)进行高分辨率插值。将方程(11)代入到(10)会发现其是无条件不稳定的。Lax-Friedrichs方法在无条件不稳定格式的基础上，将其改写为

(12)\[\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}\]

其中\(\mathbf{d}\)表示网格的距离矢量。对比方程(11)，因为方程(12)右边的后两项的存在（考虑正交网格，其为一个系数为\(\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格式，网格界面的通量这样计算：

(13)\[\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}\]

参考方程(6)，其中的\(\max(a_{f})\)为网格面上的局部的波的最大传播速度。将\(\max(a_{f})\)求出后，代入到方程(13)即为使用KT格式离散的连续性方程。对其求解有密度。

注意：

从方程(13)可以看出，密度被传输的速度为\(\phi_\own^f+\max(a_{f})|\bfS_f|\)以及\(\phi_\nei^f-\max(a_{f})|\bfS_f|\)的函数，因此对于这一类求解器，计算库郎数的时候，速度要通过下式进行：

\[ \max (|\phi_\own^f+\max(a_{f})|\bfS_f||, |\phi_\nei^f-\max(a_{f})|\bfS_f||) \]

下面考虑动量方程，其可以写为：

(14)\[ \frac{\partial \rho\mathbf{U}}{\partial t}+\nabla \cdot (\rho\mathbf{U}\mathbf{U})+\nabla p=0 \]

通常情况下建议采取守恒变量的方程，我们用\(\mathcal{U}\)来表示\(\rho\mathbf{U}\)，上式可以写为：

\[ \frac{\partial \mathcal{U}}{\partial t}+\nabla \cdot (\mathcal{U}\mathbf{U})+\nabla p=0 \]

将其进行有限体积法显性离散有：

(15)\[ (\mathcal{U}^{n+1}-\mathcal{U}^n) + \frac{\Delta t}{\Delta V} \sum\left(\phi_f\mathcal{U}_f+\bfS_fp_f\right)=0 \]

对于其中的\( \bfS_f p_f\)，参考方程(11)，将其写为

(16)\[ p_f\bfS_f=\frac{1}{2}\bfS_f(p_\own^f+p_\nei^f) \]

对于其中的\(\phi_f\mathcal{U}_f\)，参考方程(11)，将其写为

(17)\[ \phi_f\mathcal{U}_f= \frac{1}{2}( \phi^f_\own\mathcal{U}^f_\own+ \phi^f_\nei\mathcal{U}^f_\nei) \]

继续参考方程(13)，其可以写为

(18)\[ \phi_f\mathcal{U}_f= \frac{1}{2}\left( \left(\phi_\own^f+\max(a_{f})|\bfS_f|\right)\mathcal{U}^f_\own + \left(\phi_\nei^f-\max(a_{f})|\bfS_f|\right)\mathcal{U}^f_\nei \right) \]

方程(16)与方程(18)的加和可以整理为：

(19)\[\begin{split} p_f\bfS_f+\phi_f\mathcal{U}_f&= \\\\ \frac{1}{2}\bfS_f(p_\own^f&+p_\nei^f) + \frac{1}{2}\left( \left(\phi_\own^f+\max(a_{f})|\bfS_f|\right)\mathcal{U}^f_\own + \left(\phi_\nei^f-\max(a_{f})|\bfS_f|\right)\mathcal{U}^f_\nei \right) \end{split}\]

将方程(19)带入到方程(1)中有\(\mathcal{U}^{t+\dt}\)。随后，有更新的速度：

(20)\[ \bfU^{t+\dt}=\frac{\mathcal{U}^{t+\dt}}{\rho^{t+\dt}} \]

下面来看能量方程，在无粘的情况下：

(21)\[ \frac{\p \rho E}{\p t}+\nabla\cdot(\bfU \rho E) = -\nabla\cdot(\bfU p) \]

类似的，将其改为为守恒变量的形式：

(22)\[ \frac{\p \mathcal{E}}{\p t} +\nabla\cdot(\bfU (\mathcal{E}+p)) =0 \]

(23)\[ \mathcal{E}=\rho E \]

其离散形式为：

(24)\[ (\mathcal{E}^{n+1}-\mathcal{E}^n) + \frac{\Delta t}{\Delta V} \sum\left(\phi_f\mathcal{E}_f+\phi_fp_f\right) =0 \]

对于其中的\( \phi_f\mathcal{E}_f\)，参考方程(11)，将其写为

(25)\[ \phi_f\mathcal{E}_f = \frac{1}{2}\left( \left(\phi_\own^f+\max(a_{f})|\bfS_f|\right) \mathcal{E}^f_\own + \left(\phi_\nei^f-\max(a_{f})|\bfS_f|\right)\mathcal{E}^f_\nei \right) \]

然而对于\(\nabla\cdot(\bfU p)\)这一项，而是参考方程(11)进行的离散：

(26)\[ \phi_f p_f=\frac{1}{2} \left( \phi_\own^fp_\own^f+\phi_\nei^fp_\nei^f \right) \]

结合方程(25)与(26)有：

(27)\[\begin{split} \phi_f\mathcal{E}_f+\phi_f p_f& =\frac{1}{2}\left( \left(\phi_\own^f+\max(a_{f})|\bfS_f|\right) \mathcal{E}^f_\own + \left(\phi_\nei^f-\max(a_{f})|\bfS_f|\right)\mathcal{E}^f_\nei \right) \\\\ &+ \frac{1}{2} \left( \phi_\own^fp_\own^f+\phi_\nei^fp_\nei^f \right) \\\\ &= \frac{1}{2}\left( \left(\phi_\own^f+\max(a_{f})|\bfS_f|\right) (\mathcal{E}^f_\own+p^f_\own)\right) \\\\ & + \frac{1}{2}\left( \left(\phi_\nei^f-\max(a_{f})|\bfS_f|\right) (\mathcal{E}^f_\nei+p^f_\nei) \right) \\\\ &-\frac{1}{2}\max(a_{f})|\bfS_f|p^f_\own +\frac{1}{2}\max(a_{f})|\bfS_f|p^f_\nei \end{split}\]

将方程(27)代入到(24)中，可以获得下一个时间步的\(\mathcal{E}\)。

附加粘性

欧拉方程附加粘性之后，变成了NS方程，其数学特征也会存在变化。为了尽可能的借用双曲欧拉方程的特征并使用密度基求解器。rhoCentralFoam采用另外一种策略来考虑欧拉方程的粘度。考虑粘度后，主要是在动量方程以及能量方程添加粘度拉普拉斯项。

1. rhoCentralFoam首先按照上文讨论的过程进行求解；

2. 随后在动量方程以及能量方程中添加粘度贡献；

具体的，首先调用中心格式求解欧拉动量方程(1)，之后求解下述方程：

(28)\[ \frac{\p\rho\bfU}{\p t} - \frac{\p\rho\bfU}{\p t} - \nabla\cdot(\mu_{eff}\nabla\bfU)=0 \]

即有考虑粘度后的速度\(\bfU\)。同时，能量方程中也要进行类似的操作，在此不进行赘述。

关键代码

surfaceScalarField ap
(
    "ap",
    max(max(phiv_pos + cSf_pos, phiv_neg + cSf_neg), v_zero)
);
surfaceScalarField am
(
    "am",
    min(min(phiv_pos - cSf_pos, phiv_neg - cSf_neg), v_zero)
);

上述代码对应公式中的\(ap,am\)。

surfaceScalarField amaxSf("amaxSf", max(mag(am), mag(ap)));
amaxSf = max(mag(aphiv_pos), mag(aphiv_neg));

上述代码对应公式中的\(\max(a_f)|\bfS_f|\)。

surfaceScalarField aphiv_pos("aphiv_pos", phiv_pos - aSf);
surfaceScalarField aphiv_neg("aphiv_neg", phiv_neg + aSf);

上述代码对应公式中的\(\phi_\own^f+\max(a_{f})|\bfS_f|\)，\(\phi_\nei^f-\max(a_{f})|\bfS_f|\)。

surfaceVectorField phiUp
(
    (aphiv_pos*rhoU_pos + aphiv_neg*rhoU_neg)
  + (a_pos*p_pos + a_neg*p_neg)*mesh.Sf()
);

上述代码对应方程(18)。

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
);

上述代码对应方程(27)。

solve
(
    fvm::ddt(rho, U) - fvc::ddt(rho, U)
  - fvm::laplacian(muEff, U)
  - fvc::div(tauMC)
);

上述代码对应方程(28)。注意其中的fvm::ddt(rho, U) - fvc::ddt(rho, U)可以约去，但同时保证了矩阵的对角占优。更详细的请参考《无痛苦NS方程笔记》。