Pirozzoli[8]使用简单的开关函数将守恒型紧致格式与WENO格式相结合，得到了第1个紧致WENO混合格式。然而这一格式存在诸多问题：首先，开关函数虽能有效实现2种子格式之间的切换，但突然的切换会造成明显的数值振荡；其次，在求解双曲守恒率方程组时，Pirozzoli的格式使用了通量型(Component-wise)方法，通量型方法直接使用通量项的各个分量来重构对应维度上的数值通量，这种处理方法虽然足够简单、高效，但并不符合双曲守恒率方程组的数学特性。Ren等[9]从3个方面改进了紧致WENO混合格式：首先，用精心设计的权重算子代替开关函数，实现了2种格式之间的光滑过渡。其次，采用特征型(Characteristic-wise)方法求解双曲守恒率方程组。特征型方法依据双曲守恒率方程组的数学特性，先将通量项的数据投影到特征空间中，并在特征空间中重构数值通量，再将所得到的数值通量投影回物理空间。特征型方法的计算量较大，但可以使格式获得更高的分辨率和更好的数值稳定性。最后，Ren等[9]还使用了经过熵修正的Roe型格式，耗散性比通量分裂方法更低。本文将Ren等[9]发展的这一混合格式称为HCW-R。
此后，研究者们又对紧致型混合格式做出了一系列改进。王来和吴颂平[10]构造了无自由参数的权重算子。李彦苏等[11]采用耗散性更低的紧致重构加权基本无振荡(CRWENO)格式来代替WENO子格式，以进一步提高特征型紧致WENO混合格式[9]的分辨率。Peng和Shen[12]使用新的权重算子结合五阶CRWENO格式和七阶紧致格式，构造了一个具有七阶精度的紧致型混合格式。
数值试验表明，特征型紧致WENO混合格式[9]具有十分优异的分辨率特性，但由于需要处理块状三对角方程组，其计算代价非常高昂。针对这一问题，武从海等[13]改进了这一混合格式[9]：只在计算WENO通量时使用特征型方法，对于紧致格式和混合权重的计算则都使用通量型方法，从而格式只需要处理三对角方程组，计算效率大大提高。然而这一方法并不符合双曲守恒率方程组的数学特性，会降低格式的分辨率并产生数值振荡。于是文献[13]又通过修正混合权重、增大WENO子格式的比重来抑制数值振荡，但这又进一步折损了格式的分辨率。
本文用傅德薰和马延文[14]的五阶迎风紧致格式CS5-F代替文献[8-9]所使用的紧致格式CS5-P来构造新的混合格式；结合通量分裂方法以增强格式的稳定性；并仍然采用特征型方法。新格式被命名为HCW-E。相较于HCW-R，HCW-E的分辨率稍差；但由于可推进求解而无需处理三对角或块状三对角方程组，其计算效率显著提高。综合而言，当花费相同的计算时间，HCW-E可以获得显著优于HCW-R的数值结果。
1 控制方程及其离散 1.1 Euler方程组 一维Euler方程组形式如下：
 | (1) |
式中：
 | (2) |
 | (3) |
其中：U为守恒变量；F为通量项；ρ、u和p分别为流体的密度、速度和压强；E为单位体积流体的总能量；γ为气体绝热指数，通常取γ=1.4。
Jacobian矩阵A=?F/?U有3个不同的特征值λ(i)(i=1, 2, 3)和相应的3个线性无关的特征向量r(i)(i=1, 2, 3)。3个特征向量组成右矩阵R=(r(1), r(2), r(3))，左矩阵L=R－1的行向量用符号l(i)(i=1, 2, 3)表示。将Jacobian矩阵对角化LAR=Λ，Λ=diag (λ(1), λ(2), λ(3))，方程组可转化为在方向r(i)(i=1, 2, 3)上的3个独立的波动方程：
 | (4) |
1.2 数值离散 考虑间距为Δx的均匀网格G={j|j=1, 2, …, n}。xj(j=1, 2, …, n)是网格节点，xj+1/2=
步骤1 ?在网格界面xj+1/2上获得Jacobian矩阵的某种平均，Aj+1/2=A(Uj, Uj+1)。本文采用Roe平均，也可使用简单的算术平均。
步骤2 ?对模板Sj+1/2={k|k=j－l, j－l+1, …, j+m}上的守恒变量和通量项进行局部特征投影，投影后的守恒变量和通量项分别为Qk=Lj+1/2Uk和Wk=Lj+1/2Fk。Qk的分量表示为qk(i)=lj+1/2(i)Uk，Wk的分量表示为wk(i)=lj+1/2(i)Fk, i=1, 2, 3。
步骤3 ?使用通量分裂技术，将投影后的通量项分为正负两部分，wk(i)±=wk(i)±λj+1/2(i)*qk(i), i=1, 2, 3。对于当地Lax-Friedrich (LLF)分裂，
步骤4 ?在3个特征方向上分别对正负通量进行迎风重构，得到网格界面xj+1/2上的数值通量
步骤5 ?将所得的数值通量投影回物理空间，得到
步骤6 ?使用适宜的时间推进方法求解半离散化方程组：
 | (5) |
本文采用3阶总变差不增(TVD)的Runge-Kutta方法[2]。
2 数值格式 2.1 紧致格式 Ren等[9]采用CS5-P[8](式(6))来构造混合格式：
 | (6) |
式中：fj为网格节点xj上的通量值;
CS5-P具有十分优异的分辨率特性，但耦合3个网格界面上的数值通量使得该格式的计算代价非常高昂：当求解m维双曲守恒律系统时，通量型方法需求解m个三对角方程组，特征型方法则需求解块状三对角方程组。
本文采用CS5-F[14] (式(7))来构造混合格式：
 | (7) |
CS5-F可写为迎风推进的形式：
 | (8) |
在实际应用中，先使用显式格式计算边界一侧的数值通量，再按照公式向流场内部推进求解。因此，CS5-F虽然作为紧致格式具有显著优于同阶精度显式格式的分辨率，其计算代价却与显式格式无异。
以上2种紧致格式的分辨率特性参见文献[8, 14]。CS5-P的分辨率要略优于CS5-F。但由于CS5-F无需求解三对角或块状三对角方程组，其计算效率要比CS5-P高得多。
2.2 WENO格式 五阶WENO格式的数值通量是3个低阶数值通量的加权平均：
 | (9) |
式中：
 | (10) |
权重ωk, k=1, 2, 3的计算方式决定了WENO格式的数值性能。特征型紧致WENO混合格式[7]使用WENO-JS格式[2]作为子格式。本文则采用耗散更低的WENO-Z格式[4]，其权重形式如下：
 | (11) |
式中:(d0, d1, d2)=(1/10, 6/10, 3/10)为理想权重；ε=10－40用来避免分母为0；a≥1可调控格式的数值耗散，本文取a=1以减小数值耗散；βk为光滑度量因子，其具体形式如下：
 | (12) |
2.3 混合格式 文献[9]给出如下形式的权重算子：
 | (13) |
式中：σj+1/2为混合格式中紧致格式所占的比重；0 < rc≤1为一自由参数，选取较小的rc可减少WENO通量的使用率，从而能提高格式的分辨率和计算效率，但也有可能影响格式的稳定性，本文取rc=0.5，以平衡格式的分辨率与稳定性。
本文使用文献[9]的权重算子结合CS5-F和WENO-Z格式，得到如式(14)形式的混合格式：
 | (14) |
在实际应用中，先使用WENO-Z格式计算边界一侧的数值通量，再按照式(14)向流场内部推进求解。将采用CS5-F和CS5-P混合格式分别标识为HCW-E和HCW-R。2种混合格式均采用同样的WENO-Z格式作为激波捕捉子格式。注意到，在计算权重σj+1/2时，HCW-E所使用的模板与HCW-R的一样，因此在这2种混合格式中紧致格式所占的比重也是相当的。可以预见，当网格数目相同时，HCW-R的分辨率必然高于HCW-E。但由于HCW-E格式求解方式更为高效，其有可能以相同计算代价获得优于HCW-R的数值结果。
3 数值试验和结果分析 本节将用几个数值试验来验证新格式的数值性能和计算效率。
3.1 Lax问题 在区间[－5, 5]上求解一维Euler方程组，初始条件为[15]
 | (15) |
式中：ρ、u、p分别为无量纲的密度、x方向速度、压强。网格数目为N=100，CFL=0.2，计算终止时间(无量纲)为t=1.3。图 1给出计算结果的密度分布曲线，并与精确解(标记为“exact”)作对比。表 1给出不同格式的计算耗时。
 |
| 图 1 Lax问题的密度分布曲线 Fig. 1 Density distribution curves of Lax problem |
| 图选项 |
表 1 Lax问题的计算耗时 Table 1 Computational time for Lax problem
| 格式 | WENO-Z(N=100) | HCW-R(N=100) | HCW-E(N=100) | HCW-E(N=120) |
| 计算耗时/s | 2.03 | 3.20 | 1.16 | 1.61 |
表选项
图 1(a)比较了HCW-E格式和WENO-Z格式的计算结果。HCW-E的计算结果未见任何数值振荡。这说明WENO子格式在间断处占主导，赋予混合格式ENO特性。同时，HCW-E的结果显著优于WENO-Z格式，这说明尽管在间断处的比重十分微小，低耗散紧致格式的加入仍能显著改善格式对间断(尤其是接触间断)分辨率。图 1(b)比较了HCW-E和HCW-R这2种混合格式的计算结果。当网格数目同为N=100时，两者的计算结果几乎没有差异，但后者的计算耗时要显著高于前者，几乎是前者的3倍(见表 1)。另一方面，网格数目为N=120时，HCW-E的计算结果要明显优于N=100时HCW-R的计算结果，计算耗时却显著小于后者(仅为后者的1/2)。注意到，当网格数目相同时，HCW-E的计算耗时甚至比WENO-Z格式还要小得多。这是因为当取rc < 1时，在流场数据足够光滑(rj+1/2/rc≥1)的地方，混合格式无需求解WENO通量，从而大大节省了计算量。
3.2 Osher-Shu问题 在区间[－5, 5]上求解一维Euler方程组，初始条件为[16]
 | (16) |
式中：x为坐标值。CFL=0.2，计算终止时间为t=1.8(无量纲)。图 2给出计算结果的密度分布曲线。精确解(exact)是网格数目N=2 000时WENO-Z格式的计算结果。表 2给出不同格式的计算耗时。
 |
| 图 2 Osher-Shu问题的密度分布曲线 Fig. 2 Density distribution curves of Osher-Shu problem |
| 图选项 |
表 2 Osher-Shu问题的计算耗时 Table 2 Computational time for Osher-Shu problem
| 格式 | WENO-Z(N=200) | HCW-R(N=200) | HCW-E(N=200) | HCW-E(N=250) |
| 计算耗时/s | 7.81 | 17.57 | 6.36 | 9.83 |
表选项
图 2(a)比较了在网格数目同为N=200时WENO-Z格式与HCW-E格式的数值结果。两者的差异十分显著。这说明，在光滑区域紧致子格式占主导，赋予混合格式高分辨率特性。图 2(b)比较了HCW-E和HCW-R 2种混合格式的计算结果。当网格数目同为N=200时，由于HCW-R所使用的紧致格式分辨率更高，其数值结果要略优于HCW-E，但差异并不明显。而2种混合格式在计算耗时上的差异却十分明显(见表 2)。网格数目为N=200时HCW-R的计算耗时甚至比网格数目为N=250时HCW-E的计算耗时还要高得多，但后者的数值结果却要明显优于前者(见图 2(b))。综合而言，HCW-E具有更高的计算效率。
3.3 二维Riemann问题 在正方形区域[0, 1]×[0, 1]内求解二维Euler方程组，初始条件为[17]
 | (17) |
式中：v为无量纲的y方向的速度。CFL=0.2，计算终止时间为t=0.8(无量纲)。图 3给出计算结果的密度等值线，等值线分布是从0.14~1.7之间等间距的40条。表 3给出不同格式的计算耗时。
 |
| 图 3 二维Riemann问题的密度等值线 Fig. 3 Density contours of 2D Riemann problem |
| 图选项 |
表 3 二维Riemann问题的计算耗时 Table 3 Computational time for 2D Riemann problem
| 格式 | WENO-Z(N=400×400) | HCW-R(N=400×400) | HCW-E(N=400×400) | HCW-E(N=600×600) |
| 计算耗时/s | 5188 | 23439 | 4473 | 17347 |
表选项
相比于一维Riemann问题，二维问题的流场结构要复杂得多。WENO-Z格式无力捕捉精细的流场结构，紧致格式的加入则极大地提升了格式对细微结构的分辨率。当网格数目同为N=400×400时，HCW-R和HCW-E的计算结果差别不大。而当N=600×600时，HCW-E的计算结果显然优于N=400×400时HCW-R的计算结果，但计算耗时却比后者小得多(见表 3)。
3.4 双马赫反射问题 在计算区域[0, 4]×[0, 1]内求解二维Euler方程组，初始条件为[18]
 | (18) |
计算域的上边界为激波传播的精确解，左边界及底部x≤1/6部分为入流边界条件，右边界为出流边界条件，底部x>1/6的部分为壁面边界条件。计算终止时间(无量纲)为t=0.2，CFL=0.2。表 4给出不同格式的计算耗时。图 4给出双马赫杆附近区域的密度等值线，等值线为2~22之间的50等份。
表 4 双马赫反射问题的计算耗时 Table 4 Computational time for double Mach reflection problem
| 格式 | WENO-Z(N=960×240) | HCW-R(N=960×240) | HCW-E(N=960×240) | HCW-E(N=1 600×400) |
| 计算耗时/s | 7371 | 31460 | 6804 | 28118 |
表选项
 |
| 图 4 双马赫杆区域的密度等值线 Fig. 4 Density contours at double Mach stems region |
| 图选项 |
在双马赫反射问题中，WENO-Z格式的计算结果与混合格式相比差异巨大。当网格数目同为N=960×240时，HCW-R的数值结果确实优于HCW-E，但计算耗时却比N=1 600×400时的HCW-E还高(见表 4)。而后者的计算结果显然优于前者。在本问题中，HCW-E的计算耗时仍然小于WENO-Z，但差异不如在二维Riemann中明显。这是因为在双马赫反射问题中间断占主导的区域更大，从而WENO-Z格式的比重也更大。注意到，2种混合格式在计算效率上的差异随着方程组维数的增大而增大。在一维问题中，HCW-R的计算耗时大约是HCW-E的2~3倍；而在二维问题中，前者的耗时约为后者的5倍。可以预期，在三维问题中两者的差异将变得更大。
4 结论 1) 本文运用文献[9]的权重算子将文献[14]的紧致格式(CS5-F)与WENO-Z格式相结合，构造了一个高效率的特征型紧致WENO混合格式(HCW-E)。
2) 相比于采用CS5-P的混合格式(HCW-R)[9]，HCW-E的分辨率稍差。但由于HCW-E可推进求解而无需处理块状三对角方程组，其计算效率要远远高于HCW-R；从而花费相同的计算代价，HCW-E可以获得比HCW-R好得多的数值结果。新格式的效率优势在求解高维问题时更为明显。
参考文献
| [1] | LIU X D, OSHER S, CHAN T. Weighted essentially non-oscillatory schemes[J]. Journal of Computational Physics, 1994, 115(1): 200-212. DOI:10.1006/jcph.1994.1187 |
| [2] | JIANG G S, SHU C W. Efficient implementation of weighted ENO schemes[J]. Journal of Computational Physics, 1996, 126(1): 202-228. |
| [3] | HENRICK A K, ASLAM T D, POWERS J M. Mapped weighted essentially non-oscillatory schemes:Achieving optimal order near critical points[J]. Journal of Computational Physics, 2005, 207(2): 542-567. DOI:10.1016/j.jcp.2005.01.023 |
| [4] | BORGES R, CARMONA M, COSTA B, et al. An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws[J]. Journal of Computational Physics, 2008, 227(6): 3191-3211. DOI:10.1016/j.jcp.2007.11.038 |
| [5] | ACKER F, BORGES R, COSTA B, et al. An improved WENO-Z scheme[J]. Journal of Computational Physics, 2016, 313(1): 726-753. |
| [6] | 骆信, 吴颂平. 改进的五阶WENO-Z+格式[J]. 力学学报, 2019, 51(6): 1927-1939. LUO X, WU S P. An improved fifth-order WENO-Z+ scheme[J]. Chinese Journal of Theoretical and Applied Mechanics, 2019, 51(6): 1927-1939. (in Chinese) |
| [7] | LELE S K. Compact finite difference schemes with spectral-like resolution[J]. Journal of Computational Physics, 1992, 103(1): 16-42. |
| [8] | PIROZZOLI S. Conservative hybrid compact-WENO schemes for shock-turbulence interaction[J]. Journal of Computational Physics, 2002, 178(1): 81-117. |
| [9] | REN Y X, LIU M E, ZHANG H X. A characteristic-wise hybrid compact-WENO scheme for solving hyperbolic conservation laws[J]. Journal of Computational Physics, 2003, 192(2): 365-386. |
| [10] | 王来, 吴颂平. 无自由参数型混合格式[J]. 北京航空航天大学学报, 2015, 41(2): 318-322. WANG L, WU S P. Hybrid finite difference schemes without free parameters[J]. Journal of Beijing University of Aeronautics and Astronautics, 2015, 41(2): 318-322. (in Chinese) |
| [11] | LI Y S, YAN C, YU J, et al. A new high-accuracy scheme for compressible turbulent flows[J]. International Journal of Computational Fluid Dynamics, 2017, 31(9): 362-378. DOI:10.1080/10618562.2017.1365844 |
| [12] | PENG J, SHEN Y Q. A novel weighting switch function for uniformly high-order hybrid shock-capturing schemes[J]. International Journal for Numerical Methods in Fluids, 2017, 83(9): 681-703. DOI:10.1002/fld.4285 |
| [13] | 武从海, 赵宁, 田琳琳. 一种改进的紧致WENO混合格式[J]. 空气动力学学报, 2013, 31(4): 477-481. WU C H, ZHAO N, TIAN L L. An improved hybrid compact-WENO scheme[J]. Acta Aerodynamica Sinica, 2013, 31(4): 477-481. (in Chinese) |
| [14] | FU D X, MA Y W. A high order accurate difference scheme for complex flow fields[J]. Journal of Computational Physics, 1997, 134(1): 1-15. |
| [15] | LAX P D. Weak solutions of nonlinear hyperbolic equations and their numerical computation[J]. Communications on Pure and Applied Mathematics, 1954, 7(1): 159-193. DOI:10.1002/cpa.3160070112 |
| [16] | SHU C W, OSHER S. Efficient implementation of essentially non-oscillatory shock-capturing schemes, Ⅱ[J]. Journal of Computational Physics, 1989, 83(1): 32-78. |
| [17] | LAX P D, LIU X D. Solution of two-dimensional Riemann problems of gas dynamics by positive schemes[J]. SIAM Journal on Scientific Computing, 1998, 19(2): 319-340. |
| [18] | WOODWARD P, COLELLA P. The numerical simulation of two-dimensional fluid flow with strong shocks[J]. Journal of Computational Physics, 1984, 54(1): 115-173. |
