0 引言

近年来随着流体力学的迅速发展，业内人士提出了许多计算流体力学的数值方法^[1-2]。迎风格式以其高分辨率和清晰的物理解而深受相关人士的青睐，其中低耗散E-CUSP是一种迎风对流差分格式^[3]，主要是将无黏通量分为守恒项和压力项，避免了矩阵运算，使得运算效率大幅度提高，但是其在间断处抹平现象比较严重，分辨率不高，为此需要耦合加权本质无振荡(WENO)格式^[4-7]来提高分辨率。WENO格式利用各个备选模板的凸组合方式进行重构，并且每个模板权重的选取依赖于该模板光滑因子的局部光滑性。Embedded-WENO格式^[8]的主要思想是在间断处WENO格式表现为低阶精度，但是仍然可以使用剩余的光滑模板重构间断模板，同时引入内模板和外模板，使得在间断处格式收敛到精度较高的内模板，在光滑处收敛到外模板，这样可以很好地消除截断误差^[9]和提高分辨率。本文通过数值实验证明了所提出的新格式E-CUSP-Embedded-WENO5具有良好的收敛效果和更高的分辨率，为二维欧拉方程的进一步研究奠定基础。

1 数值方法 1.1 控制方程

控制方程为一维欧拉方程，其守恒形式为

$ \frac{{\partial \mathit{\boldsymbol{U}}}}{{\partial t}} + \frac{{\partial \mathit{\boldsymbol{F}}}}{{\partial x}} = 0, $

(1)

式中：U=[ρ, ρu, ρe]^T为守恒型变量，ρ为密度，u为速度，e为单位质量的内能；F=[ρu, ρu²+p, (ρe+p)u]^T为数值通量，p为压强且$ p = (\gamma - 1)(\rho e - \frac{1}{2}\rho {u^2})$，γ为比热率且γ=1.4; x为空间指标；t为时间指标。

1.2 E-CUSP格式

对方程(1)在区间单元上利用有限差分法可得$\frac{{\mathit{\boldsymbol{U}}_i^{n + 1} - \mathit{\boldsymbol{U}}_i^n}}{{\Delta t}} = - \frac{1}{{\Delta x}}(\mathit{\boldsymbol{F}}_{i + \frac{1}{2}}^n - \mathit{\boldsymbol{F}}_{i - \frac{1}{2}}^n)$, 式中：Δx为空间步长；Δt为时间步长；U_iⁿ和U_i ^{n+ 1}分别表示U在t=tⁿ和t=tⁿ⁺¹时刻的空间平均值, n为时间层, i为空间节点；$\mathit{\boldsymbol{F}}_{i - \frac{1}{2}}^n$和$\mathit{\boldsymbol{F}}_{i + \frac{1}{2}}^n$分别表示穿越单元边界${x_{i - \frac{1}{2}}}$和${x_{i + \frac{1}{2}}}$的通量。E-CUSP格式的基本思想是将数值通量分解为守恒型通量和压力型通量，以避免矩阵运算，从而提高计算效率。方程(1)中雅克比矩阵可表示为

$ \mathit{\boldsymbol{A}} = \frac{{\partial \mathit{\boldsymbol{F}}}}{{\partial \mathit{\boldsymbol{U}}}} = \mathit{\boldsymbol{T \boldsymbol{\varLambda} }}{\mathit{\boldsymbol{T}}^{ - 1}}, $

式中：$\mathit{\boldsymbol{ \boldsymbol{\varLambda} }} = \left[ {\begin{array}{*{20}{c}} {u - a}&0&0\\ 0&u&0\\ 0&0&{u + a} \end{array}} \right]$；T为特征向量。进而，可得F=TΛT^-1U。因此，有

$ \mathit{\boldsymbol{F}} = \mathit{\boldsymbol{T}}\left[ {\begin{array}{*{20}{l}} u&0&0\\ 0&u&0\\ 0&0&u \end{array}} \right]{\mathit{\boldsymbol{T}}^{ - 1}}\mathit{\boldsymbol{U}} + \mathit{\boldsymbol{T}}\left[ {\begin{array}{*{20}{r}} { - a}&0&0\\ 0&0&0\\ 0&0&a \end{array}} \right]{\mathit{\boldsymbol{T}}^{ - 1}}\mathit{\boldsymbol{U}} = {\mathit{\boldsymbol{F}}^{\rm{c}}} + {\mathit{\boldsymbol{F}}^{\rm{p}}}, $

式中：F^c=u[ρ, ρu, ρe]^T表示守恒型通量；F^p=[0, p, ρu]^T表示压力型通量。这样数值通量F就被分解为守恒型通量F^c与压力型通量F^p之和。守恒型通量F^c可表示为${\mathit{\boldsymbol{F}}^{\rm{c}}} = \frac{1}{2}[{(\rho u)_{\frac{1}{2}}}(\mathit{\boldsymbol{U}}_{\rm{L}}^{\rm{c}} + \mathit{\boldsymbol{U}}_{\rm{R}}^{\rm{c}}) - |\rho u|(\mathit{\boldsymbol{U}}_{\rm{R}}^{\rm{c}} - \mathit{\boldsymbol{U}}_{\rm{L}}^{\rm{c}})]$, 式中：U_L^c表示U^c在节点${x_{i - \frac{1}{2}}}$处值，U_R^c表示U^c在节点$ {x_{i + \frac{1}{2}}}$处值，${\mathit{\boldsymbol{U}}^{\rm{c}}} = {[1, u, e]^{\rm{T}}};{(\rho u)_{\frac{1}{2}}}$表示交界面处质量通量且${(\rho u)_{\frac{1}{2}}} = ({\rho _{\rm{L}}}u_{\rm{L}}^ + + {\rho _{\rm{R}}}u_{\rm{R}}^ - )$, $u_{{\rm{L(R)}}}^ \pm = {a_{\frac{1}{2}}}\{ \frac{{{M_{{\rm{L(R)}}}} \pm |{M_{{\rm{L(R)}}}}|}}{2} + {\alpha _{{\rm{L(R)}}}}[ \pm \frac{1}{4}{({M_{{\rm{L(R)}}}} \pm 1)^2} - \frac{{{M_{{\rm{L(R)}}}} \pm |{M_{{\rm{L(R)}}}}|}}{2}]\} $, ${a_{\frac{1}{2}}} = \frac{1}{2}({a_{\rm{L}}} + {a_{\rm{R}}}), {M_{{\rm{L(R)}}}} = \frac{{{u_{{\rm{L(R)}}}}}}{{{a_{\frac{1}{2}}}}}, {\alpha _{{\rm{L(R)}}}} = \frac{{2{{({h_t}/\rho )}_{{\rm{L(R)}}}}}}{{{{({h_t}/\rho )}_{{\rm{L(R)}}}} + {{({h_t}/\rho )}_{{\rm{R(L)}}}}}}$, 其中h_t为总焓且${h_t} = e + \frac{p}{\rho }$。对于压力型通量F^p，动量方程中的压力可表示为p=(P⁺p)_L+(P^-p)_R, 其中${P^ \pm } = \frac{1}{4}{(M \pm 1)^2}(2 \pm M) \pm \alpha M{({M^2} - 1)^2}, \alpha = 3/16$；能量方程中的压力可分解为$pu = \frac{1}{2}{[p(u + {a_{\frac{1}{2}}})]_{\rm{L}}} + \frac{1}{2}{[p(u - {a_{\frac{1}{2}}})]_{\rm{R}}}$。最终E-CUSP格式中数值通量F可表示为

$ \mathit{\boldsymbol{F}} = \frac{1}{2}[{(\rho u)_{\frac{1}{2}}}(\mathit{\boldsymbol{U}}_{\rm{L}}^{\rm{c}} + \mathit{\boldsymbol{U}}_{\rm{R}}^{\rm{c}}) - |\rho u{|_{\frac{1}{2}}}(\mathit{\boldsymbol{U}}_{\rm{R}}^{\rm{c}} - \mathit{\boldsymbol{U}}_{\rm{L}}^{\rm{c}})] + [0, {P^ + }p, \frac{1}{2}p(u + {a_{\frac{1}{2}}})]_{\rm{L}}^{\rm{T}} + [0, {P^ - }p, \frac{1}{2}p(u - {a_{\frac{1}{2}}})]_{\rm{R}}^{\rm{T}}。$

1.3 嵌入WENO格式 1.3.1 WENO格式

五阶WENO重构为$\mathit{\boldsymbol{U}}_{i + \frac{1}{2}}^{{\rm{L(R)}}} = {\omega _0}{q_0} + {\omega _1}{q_1} + {\omega _2}{q_2}$, 其中q_k(k=0, 1, 2)表示不同模板的重建。非线性权重${\omega _k} = \frac{{{\alpha _k}}}{{{\alpha _0} + {\alpha _1} + {\alpha _2}}}, {\alpha _k} = \frac{{{d_k}}}{{{{({\beta _k} + \varepsilon )}^2}}}$, 其中d_k(k=0, 1, 2)为线性权重，d₀=0.1, d₁=0.6, d₂=0.3; ε是为了避免分母变为0，通常视为10^-6; β_k为光滑因子，经典光滑因子β_k ^[4]定义为

$ \left\{ \begin{array}{l} {\beta _0} = \frac{{13}}{{12}}{({\mathit{\boldsymbol{U}}_{i - 2}} - 2{\mathit{\boldsymbol{U}}_{i - 1}} + {\mathit{\boldsymbol{U}}_i})^2} + \frac{1}{4}{({\mathit{\boldsymbol{U}}_{i - 2}} - 4{\mathit{\boldsymbol{U}}_{i - 1}} + 3{\mathit{\boldsymbol{U}}_i})^2}, \\ {\beta _1} = \frac{{13}}{{12}}{({\mathit{\boldsymbol{U}}_{i - 1}} - 2{\mathit{\boldsymbol{U}}_i} + {\mathit{\boldsymbol{U}}_{i + 1}})^2} + \frac{1}{4}{({\mathit{\boldsymbol{U}}_{i - 1}} - {\mathit{\boldsymbol{U}}_{i + 1}})^2}, \\ {\beta _2} = \frac{{13}}{{12}}{({\mathit{\boldsymbol{U}}_i} - 2{\mathit{\boldsymbol{U}}_{i + 1}} + {\mathit{\boldsymbol{U}}_{i + 2}})^2} + \frac{1}{4}{(3{\mathit{\boldsymbol{U}}_i} - 4{\mathit{\boldsymbol{U}}_{i + 1}} + {\mathit{\boldsymbol{U}}_{i + 2}})^2}。\end{array} \right. $

(2)

由非线性权重ω_k可得$\frac{{{\omega _k}}}{{{\omega _l}}} = \frac{{{\alpha _k}}}{{{\alpha _l}}} = \frac{{{d_k}}}{{{d_l}}}{(\frac{{{\beta _l} + \varepsilon }}{{{\beta _k} + \varepsilon }})^2}$, 其中k, l=0, 1, 2。光滑因子β_k的Taylor展开式为

$ \left\{ \begin{array}{l} {\beta _0} = {(\mathit{\boldsymbol{U}}_i^\prime )^2}\Delta {x^2} + (\frac{{12}}{{13}}{(\mathit{\boldsymbol{U}}_i^{\prime \prime })^2} - \frac{2}{3}\mathit{\boldsymbol{U}}_i^\prime \mathit{\boldsymbol{U}}_i^{\prime \prime \prime })\Delta {x^4} - (\frac{{13}}{6}\mathit{\boldsymbol{U}}_i^{\prime \prime }\mathit{\boldsymbol{U}}_i^{\prime \prime \prime } - \frac{1}{2}\mathit{\boldsymbol{U}}_i^\prime \mathit{\boldsymbol{U}}_i^{\prime \prime \prime })\Delta {x^5} + O(\Delta {x^6}), \\ {\beta _1} = {(\mathit{\boldsymbol{U}}_i^\prime )^2}\Delta {x^2} + (\frac{{12}}{{13}}{(\mathit{\boldsymbol{U}}_i^{\prime \prime })^2} + \frac{2}{3}\mathit{\boldsymbol{U}}_i^\prime \mathit{\boldsymbol{U}}_i^{\prime \prime \prime })\Delta {x^4} + O(\Delta {x^6}), \\ {\beta _2} = {(\mathit{\boldsymbol{U}}_i^\prime )^2}\Delta {x^2} + (\frac{{12}}{{13}}{(\mathit{\boldsymbol{U}}_i^{\prime \prime })^2} - \frac{2}{3}\mathit{\boldsymbol{U}}_i^\prime \mathit{\boldsymbol{U}}_i^{\prime \prime \prime })\Delta {x^4} + (\frac{{13}}{6}\mathit{\boldsymbol{U}}_i^{\prime \prime }\mathit{\boldsymbol{U}}_i^{\prime \prime \prime } + \frac{1}{2}\mathit{\boldsymbol{U}}_i^\prime \mathit{\boldsymbol{U}}_i^{\prime \prime \prime \prime })\Delta {x^5} + O(\Delta {x^6})。\end{array} \right. $

(3)

由式(3)可得

$ \frac{{{\beta _k}}}{{{\beta _l}}} = \left\{ {\begin{array}{*{20}{l}} {1 + O(\Delta {x^3}), }&{k = 0, l = 2{\rm{ 或 }}k = 2, l = 0, }\\ {1 + O(\Delta {x^2}), }&{其他}。\end{array}} \right. $

因此，非线性权重的比值为$\frac{{{\omega _k}}}{{{\omega _l}}} = \frac{{{d_k}}}{{{d_l}}}(1 + O(\Delta {x^s}))$, 其中s≥2。由此可知，非线性权重的比值只取决于光滑因子β_k、β_l局部光滑性和线性权重的比值。

1.3.2 嵌入WENO格式

五阶WENO格式模板在x_i处可分为三个三点模板S₀={x_i-2, x_i-1, x_i}, S₁={x_i-1, x_i, x_i+1}, S₂={x_i, x_i+1, x_i+2}, S₀、S₁、S₂即为外部模板。当S₀发生间断且S₁、S₂光滑时，间断一定发生在(x_i-2, x_i-1)处，这样可以通过S₁、S₂来构造四点模板，即内部模板S_{1, 2}:=S₁∪S₂={x_i-1, x_i, x_i+1, x_i+2}; 同理，当S₂发生间断且S₀、S₁光滑时，内部模板S_{0, 1}:=S₀∪S₁={x_i-2, x_i-1, x_i, x_i+1}。假设内部格式是由嵌入模板权重α₀⁽²⁾、α₁⁽²⁾、α₁⁽⁰⁾、α₂⁽⁰⁾所定义，其中α_m⁽ⁿ⁾上标为间断模板索引，下标为光滑模板索引，数值通量${U_{i + \frac{1}{2}}}$的凸组合为$U_{i + \frac{1}{2}}^{(0, 1)} = \alpha _0^{(2)}U_{i + \frac{1}{2}}^{(0)} + \alpha _1^{(2)}U_{i + \frac{1}{2}}^{(1)}, U_{i + \frac{1}{2}}^{(1, 2)} = \alpha _1^{(0)}U_{i + \frac{1}{2}}^{(1)} + \alpha _2^{(0)}U_{i + \frac{1}{2}}^{(2)}$。若从收敛性角度出发，选择模板权重值为α₀⁽²⁾=1/4, α₁⁽²⁾=3/4, α₁⁽⁰⁾=1/2, α₂⁽⁰⁾=1/2;若从光谱角度出发，选择模板权重值为α₀⁽²⁾=1/10, α₁⁽²⁾=9/10, α₁⁽⁰⁾=7/10, α₂⁽⁰⁾=3/10。

由非线性权重的比值可得，非线性权重ω_k是线性权重d_k的重新分配，而且ω_k的选取和光滑因子的局部光滑性有关，因此引入相关比例系数c₀、c₂, 其下标表示不光滑模板索引，可定义为α₀⁽²⁾:α₁⁽²⁾=c₂d₀:d₁, α₁⁽⁰⁾:α₂⁽⁰⁾=d₁:c₀d₂, 相应的c₀、c₂值为c₀=2, c₂=2或c₀=2/3, c₂=6/7。文献[8]提出了嵌入WENO格式非线性权重的一般形式为α_k^(E)=α_k^(O)g_k, 式中：α_k^(E)为Embedded-WENO格式未归一化非线性权重；α_k^(O)为外部格式权重；g_k为修正因子且${g_k} = {a_{kk}} + \sum\limits_{l \ne k} {\frac{{{a_{kl}}{\beta _l}}}{{{\beta _k} + \varepsilon }}} $, 其中k, l=0, 1, 2，a_kl为嵌入系数且$\sum\limits_{l = 0}^{r - 1} {{a_{kl}}} = 1, k = 0, 1, \cdots , r - 1$。嵌入方程为$\frac{{{d_k}}}{{\alpha _k^{(K)}}}\sum\limits_{q \ne K} {{a_{kq}}} = \frac{{{d_l}}}{{\alpha _l^{(K)}}}\sum\limits_{t \ne K} {{a_{lt}}} $, K表示不光滑模板索引且K为0或2，因此得到$\frac{{{d_0}}}{{\alpha _0^{(2)}}}{a_{02}} = \frac{{{d_1}}}{{\alpha _1^{(2)}}}{a_{12}}, \frac{{{d_2}}}{{\alpha _2^{(0)}}}{a_{20}} = \frac{{{d_1}}}{{\alpha _1^{(0)}}}{a_{10}}$, 简化可得$\frac{{{a_{02}}}}{{{a_{12}}}} = {c_2}, \frac{{{a_{20}}}}{{{a_{10}}}} = {c_0}$。综上可得，g_k的9个系数满足: a₀₀=1-a₀₁-a₀₂, a₁₁=1-a₂₀/c₀-a₀₂/c₂, a₂₂=1-a₂₀-a₂₁, a₁₂=a₀₂/c₂, a₁₀=a₂₀/c₀, 其中a₀₁、a₀₂、a₂₀、a₂₁可自由选择。实验证明，a₀₁、a₀₂、a₂₀、a₂₁的不同选择都会影响五阶WENO格式数值效果。通过大量实验发现a₀₁=a₂₁=0, a₀₂=c₀/3, a₂₀=c₂/3时效果最佳。将五阶WENO格式的线性权重d_k作为外部格式权重，可得Embedded-WENO格式的非线性权重α_k^(E)为

$ \left\{ {\begin{array}{*{20}{l}} {\alpha _0^{{\rm{(E)}}} = \frac{1}{3}{d_0}(3 - {c_2} + {c_2}\frac{{{\beta _2}}}{{{\beta _0} + \varepsilon }}), }\\ {\alpha _1^{{\rm{(E)}}} = \frac{1}{3}{d_1}(1 + \frac{{{\beta _2}}}{{{\beta _1} + \varepsilon }} + \frac{{{\beta _0}}}{{{\beta _1} + \varepsilon }}), }\\ {\alpha _2^{{\rm{(E)}}} = \frac{1}{3}{d_2}{{(3 - {c_0} + {c_0}\frac{{{\beta _0}}}{{{\beta _2} + \varepsilon }})}_0}}。\end{array}} \right. $

(4)

假设模板均光滑，将Taylor展开式应用到(4)中，可得ω_k^(E)-d_k=1/3(ω_k^(O)-d_k)+O(Δx³), 其中ω_k^(E)为Embedded-WENO5格式归一化非线性权重。假设只有一个模板β_k光滑，则β_k=O(Δx²), 有ω_k^(E)=1+O(Δx²)。综上两种分析，模板均为光滑或者只有一个模板光滑时，Embedded-WENO5等价于WENO5格式；当两个相邻模板光滑且第三个模板存在间断时，可以获得内部模板权重，并且Embedded-WENO5格式的实验效果优于WENO5格式。

1.4 龙格-库塔方法

Embedded-WENO格式在空间离散的同时又可以保证时间的连续性，这种方法可以将偏微分方程(PDE)转化成大量耦合的常微分方程(ODE)，即$\frac{{{\rm{d}}U}}{{{\rm{d}}t}} = L(U)$，其中L为Embedded-WENO格式的计算结果。在空间离散后，将时间t通过tⁿ=nΔt离散，在时间方向上采用四阶保持强稳定龙格-库塔(SSPRK)方法^[10]，显式SSPRK方法精度高且不会因时间积分而产生伪振荡。本文采用四阶SSPRK方法，其在一个时间步长Δt上可以表示为

$ \left\{ \begin{array}{l} \begin{array}{*{20}{l}} {{u^{(0)}} = {u^n}, }\\ {{u^{(1)}} = {u^{(0)}} + \Delta tL({u^{(0)}}), } \end{array}\\ {u^{(2)}} = \frac{1}{2}{u^{(0)}} + \frac{1}{2}{u^{(1)}} + \frac{1}{4}\Delta tL({u^{(0)}}) + \frac{1}{2}\Delta tL({u^{(1)}}), \\ {u^{(3)}} = \frac{1}{9}{u^{(0)}} + \frac{2}{9}{u^{(1)}} + \frac{2}{3}{u^{(2)}} - \frac{1}{9}\Delta tL({u^{(0)}}) - \frac{1}{3}\Delta tL({u^{(1)}}) + \Delta tL({u^{(2)}}), \\ {u^{(n + 1)}} = \frac{1}{3}{u^{(1)}} + \frac{1}{3}{u^{(2)}} + \frac{1}{3}{u^{(3)}} + \frac{1}{6}\Delta tL({u^{(1)}}) + \frac{1}{6}\Delta tL({u^{(3)}})。\end{array} \right. $

(5)

2 数值实验 2.1 一维欧拉方程Lax激波管问题

在区域[0, 1]上求解初值问题

$ (\rho , u, p) = \left\{ {\begin{array}{*{20}{l}} {(0.445, 0.698, 3.528), }&{x \le 0.5, }\\ {(0.5, 0, 0.571), }&{x > 0.5}。\end{array}} \right. $

采用Neumann边界条件，取200个网格点数，计算中条件数为0.4，本算例分别利用E-CUSP格式、E-CUSP耦合五阶WENO格式(E-CUSP-WENO5)和E-CUSP-Embedded-WENO5格式，在t=0.16时刻下求解欧拉方程激波管问题的密度和速度，结果如图 1和图 2所示。可以看出，E-CUSP-Embedded-WENO5格式在间断处明显优于E-CUSP格式和E-CUSP-WENO5格式。在光滑处三者等价，但是在间断处E-CUSP-Embedded-WENO5格式捕捉效果更好，过渡带更窄，更加接近理论解，分辨率更高。

2.2 一维欧拉方程激波管问题

在区域[0, 1]上求解初值问题

$ (\rho , u, p) = \left\{ {\begin{array}{*{20}{l}} {(5.999{\kern 1pt} {\kern 1pt} {\kern 1pt} 24, 19.597{\kern 1pt} {\kern 1pt} {\kern 1pt} 5, 460.894), }&{x \le 0.4, }\\ {(5.992{\kern 1pt} {\kern 1pt} {\kern 1pt} 42, - 6.196{\kern 1pt} {\kern 1pt} {\kern 1pt} 33, 46.095), }&{x > 0.4}。\end{array}} \right. $

采用Neumann边界条件，取200个网格点数，计算中条件数为0.4, 本算例分别利用E-CUSP格式、E-CUSP-WENO5格式和E-CUSP-Embedded-WENO5格式，在t=0.035时刻下求解欧拉方程激波管问题的密度和速度，结果如图 3和图 4所示。可以看出，E-CUSP-Embedded-WENO5格式的鲁棒性更好，比E-CUSP-WENO5格式更接近理论解，分辨率更高，尤其在激波处捕捉到的过渡单元明显减少，仅需要两到三个单元。

3 结束语

本文通过求解一维欧拉方程的典型算例，分析了低耗散E-CUSP格式耦合Embedded-WENO后所得新格式E-CUSP-Embedded-WENO5的性能，将新格式与E-CUSP -WENO5格式以及E-CUSP格式的实验结果进行了分析比较，发现E-CUSP-Embedded-WENO5格式和E-CUSP-WENO5格式在光滑区域精度相同，但是由于E-CUSP-Embedded-WENO格式改变了非线性权重，使得在间断处精度有所提高，更接近理论解，分辨率更高。新格式对激波的捕捉能力更强，并且没有增加过多的计算量，这为欧拉方程的数值模拟提供了新的备选方法，也为二维欧拉方程的进一步研究奠定了基础。