2. 中国空气动力研究与发展中心 计算空气动力研究所, 四川 绵阳 621000
2. Computational Aerodynamics Institute of China Aerodynamics Research and Development Center, Mianyang 621000, China
随着所研究的流动问题越来越复杂，研究和应用对差分格式的计算精度也提出了更高的要求。相比于传统的二阶精度格式(如TVD格式)，高阶精度格式具有更小的数值耗散和色散误差，可以更好地刻画流场的精细结构。
高阶有限差分格式是高阶精度格式的重要组成部分。国外自从Harten提出ENO[1]和Liu、Osher和Chan等提出WENO[2]格式以来，众多学者从不同的思路出发，提出各种形式的高阶格式。其中，Jiang & Shu[3]提出了一种有效地实现WENO的方式，并将其用于有限差分，其五阶形式(WENO5) 成为最流行的高阶格式之一，其非线性权求解过程及光滑度量因子(Indicator of Smoothness, IS) 形式也逐渐成为后续研究的参考标准。尽管拥有很多的优势，但WENO5还存在一定的缺陷。Henrick等[4]首先指出，WENO5在极值点不能保持五阶精度，他们分析了保持精度的充分必要条件，并提出了一种精心设计的映射函数对格式的非线性权进行修正，相应的格式称为WENO-M。WENO-M可以保持精度一致，且分辨率有所提高。对非线性权进行改进并提出新格式的还有Borges[5]等提出的WENO-Z格式，Ha[6]等提出的WENO-NS格式，Ping Fan[7]等提出的WENO-η格式，Hui Feng[8]等提出的WENO-PM格式等。最近，Li等[9]从另一种构造映射函数的思路出发，独立地提出了一种分段多项式映射函数，对WENO5格式的非线性权进行修正，提出了WENO-PPM格式。新格式在保持计算稳定性的同时，进一步提高了格式的分辨率。国内张涵信[10]较早地提出了构造差分格式应该满足的四个原则，提出了一种无振荡、无自由参数耗散算法，即二阶精度的NND格式。由于NND形式简单、物理概念清楚、易于编程和计算量小，在国内工程领域得到了广泛的应用。随后，张涵信及其团队进一步发展了一系列高精度计算格式，包括ENN[11]、WNND[12]、WENN[13]、中心型三阶格式[14]等。
相比于基础和机理研究所采用的高阶格式，工程应用的需求通常具有以下特点：1) 实际问题中网格结构复杂，希望计算格式使用较小的网格模板；2) 由于实际问题流动的复杂性，要求格式具有类TVD格式的强鲁棒性；3) 格式具有较好的耗散特性，能够提升对旋涡等流动结构的描述能力。已有计算显示(如直升机旋翼的翼梢涡、涡流发生器产生流向涡的数值模拟等)，当采用TVD格式对这些结构进行数值模拟时，结构往往很快被耗散掉[14]。因此有必要进一步发展符合工程需求的三阶格式，在保证计算稳定性的前提下，改善计算的精度和分辨能力，且能够方便地推广到有限体积程序。
本文采用映射函数的思想，将五次分段多项式映射函数应用到中心型三阶格式上，构造了新型三阶格式SWENO3-PPM5；并与Henrick等提出的映射函数作对比(SWENO3-M)。通过对典型算例的计算，验证比较了新格式的性能。
1 数值计算方法 1.1 中心型三阶格式李沁等[14]通过分析格式的耗散关系，提出了一类耗散可调的中心型三阶格式，提升了三阶格式的性能。本文基于此格式，做进一步的改进。首先简单介绍中心型三阶格式。
以一维双曲守恒律方程为例：
|
(1) |
考虑守恒型差分格式，式(1) 的数值求解可以写成如下半离散化形式：
|
(2) |
其中
|
(3) |
计算时，采用通量矢量分裂，f(u) 被分为两部分：
|
(4) |
数值通量表示为：
|
(5) |
这里仅给出
|
(6) |
本文计算中取p=2，ε=10-16。其中，ISk代表第k个模板的光滑度量因子。参考WENO系列格式非线性加权的基本框架ISk具体形式为：
|
(7) |
由于中心型三阶格式基于NND格式增加了一个下风模板，格式的线性权系数中就引入了一个自由参数，记为n3rd。线性权系数见式(8)。在约化波数为π处，中心型三阶线性格式的耗散为WENO3线性格式的1/n3rd。通过改变n3rd的值就可以调整中心型三阶格式的耗散。具体分析可参考文献[14]。本文的计算中取n3rd=10。
|
(8) |
为了增强格式的稳定性，在使用(7) 式计算得到各个模板的光滑因子IS后，要对最下风模板进行一次修正：
|
(9) |
Henrick等[4]最先指出，在极值0点(f′ j=0) 时，WENO5格式的五阶精度不能保持。他们提出了一个精度保持的充分条件(式(10))，并提出了一种设计的映射函数(式(11))。
|
(10) |
|
(11) |
通过使用映射函数来修正原来的非线性权ωk，可以使WENO5格式的精度保持五阶。这个函数定义域为[0, 1]，且满足：
|
(12) |
李沁等[9]在对映射函数的性质进行了分析后，认为映射函数应满足以下属性：经过映射以后的非线性权ω在定义域边界处应尽快地收敛到0或者1；映射函数需满足精度要求和单调性；在ω=Ck处，映射以后的非线性权分布要足够平坦和光滑。并进一步探寻拥有更好性能的可行的映射函数，提出了分段多项式映射函数。本文采用其中的五次多项式作为映射函数。具体方法介绍如下。首先将定义域分为[0, Ck]和[Ck, 1]两部分，在左、右区间分别构造多项式作为映射函数，分别记为gkL和gkR。具体构造条件如下:
对于gkL，
|
(13a) |
对于gkR，
|
(13b) |
这样在左、右两个区间就各得到一个五次多项式。假设a=ω/Ck且b=1/(Ck-1)，映射函数的形式为:
|
(14a) |
|
(14b) |
此函数将使用式(6) 计算得到的非线性权ωk进行一次映射，再将得到的值进行归一化作为最终的非线性权。图 1给出了Henrick型映射权函数与分段多项式型映射函数的分布。可以看出，相比于Henrick型映射权函数，采用五次分段多项式映射函数方法使非线性权在ω=Ck处的分布更加平坦，即更加接近于格式的理想线性权值，这将有利于提高格式分辨精细结构的能力。
|
| 图 1 映射函数分布 Fig. 1 The distribution of mapping functions |
本文将分段多项式映射函数(式(14)) 运用于中心型三阶格式，相应的格式记为SWENO3-PPM5。并与采用Henrick等提出的映射函数的三阶格式作比较，记为SWENO3-M格式。
1.3 不同格式的近似色散关系(ADR) 对比图 2给出了三种非线性格式的近似色散关系[15](Approximate Dispersion Relations, ADR) 对比，可以看到：采用映射函数方法的中心型三阶格式的色散和耗散关系得到明显的提高，且SWENO3-PPM5的表现略优于SWENO3-M。
|
| 图 2 不同格式的ADR对比 Fig. 2 Comparison of ADR |
2 数值计算结果与分析
本节计算了典型的一维和二维算例来检测新型三阶格式的性能，对比了SWENO3-PPM5、SWENO3-M、WENO3、NND四种格式。所有的算例均使用Steger-Warming分裂方式，并在计算中使用特征投影。时间推进方面，对于非定常问题使用具有TVD性质的三阶Runge-Kutta格式，对于定常问题则使用LU-SGS格式。
2.1 Shu-Osher问题该算例描述了一维情形下马赫数为3.0的右行激波和密度波发生干扰作用，其目的在于分辨主激波波后的系列小尺度结构，通常用于验证格式的分辨率。初始条件为：
|
计算采用801个网格点的均匀网格，结果如图 3和图 4所示。图中把使用WENO5在1601个网格点上的计算结果作为准确解。从计算结果来看，由于中心型三阶格式的耗散低于WENO3格式，且使用映射函数可以使格式的权值更加接近于线性权，因此SWENO3-PPM5和SWENO3-M的分辨率要高于WENO3格式和NND格式；且由于PPM5映射函数相比于Henrick提出的函数，可以使非线性权在ω=Ck处的分布更加平坦，即更加接近于线性权，因此SWENO3-PPM5的分辨率也要略高于SWENO3-M。
|
| 图 3 t=1.8时，不同格式的全局密度分布 Fig. 3 Global distribution of density at t=1.8 |
|
| 图 4 t=1.8时，不同格式的局部密度分布 Fig. 4 Local distribution of density at t=1.8 |
2.2 激波/层流边界层干扰问题
激波/层流边界层干扰流动包含丰富的流动结构，常作为检验格式性能的算例之一。来流状态为M∞=2.0，T∞=293K，Re∞=2.96×105，Pr=0.72，β=32.585°，β为激波入射角；壁面为黏性绝热条件；进出口皆为超声速边界。计算域为：0.0≤x≤2.02，0.0≤y≤1.30，网格为nx×ny=102×121，x方向均匀分布，y方向在壁面附近进行网格加密。
为了研究这几种格式分辨膨胀波系的能力，本文取y=0.3045处的压力分布(即图 5中红线位置)，并与WENO5的计算结果进行比较，如图 6所示。可以看出，在x=1.4附近，NND格式完全不能分辨该区域；WENO3格式亦不能有效分辨该膨胀波系；虽然与模板较宽的WENO5格式的结果有一定差距，但是SWENO3-PPM5和SWENO3-M格式也能够较好分辨到此区域，分辨率明显优于NND格式和WENO3格式。
|
| 图 5 SWENO3-PPM5，压力等值线分布 Fig. 5 SWENO3-PPM5, distribution of pressure contour |
|
| 图 6 沿y=0.3045处的压力分布 Fig. 6 Distribution of pressure along y=0.3045 |
2.3 激波/旋涡干扰问题
本算例数值模拟一个Ms=3.0的强激波与一个强旋涡之间的相互作用，与Borges[5]和Inoue[16]使用的算例类似。计算使用Euler方程。为了简便，将坐标系固定在激波上，并把一个等熵涡叠加在激波前的流场中，初始条件为：
|
其中uθ代表旋涡诱导的周向速度。r0=(4, 0) 代表初始时旋涡中心所在位置。γ=1.4。p(r) 和ρ(r) 分别代表激波前区域初始的压力和密度分布，初始激波面位于x=0.0处。旋涡强度Γ=0.4。以上三个量分别采用u∞、ρ∞u∞2、ρ∞进行无量纲化。激波后的初始流场采用Rankine-Hugoniot关系计算得到。计算域为[-20, 10]×[-15, 15]。并使用1501×1501的均匀网格。计算时间步长Δt=0.002，计算至无量纲时间t=14.0。
图 7给出了四种格式计算得到的密度等值线。可以看到，SWENO3-PPM5和SWENO3-M格式分辨出的流场结构明显比WENO3和NND更加丰富。由于初始流场含有激波间断，为了适应Euler方程的解算器，在x=-3.6处会形成垂直方向的二次干扰波。NND和WENO3格式由于耗散较大，并没有较好地捕捉到此干扰波；而经过映射后的中心型三阶格式则较为清晰地分辨出此处的干扰波。同样还有图中方框区域内的流场结构，可以看出，SWENO3-PPM5分辨出的结构较SWENO3-M更加细致。
|
| 图 7 不同格式的密度等值线分布 Fig. 7 Distribution of density contour |
2.4 绕25°/55°尖双锥流动问题
绕尖双锥流动是一类典型的高超声速绕流。本文对绕25°/55°尖双锥流动进行模拟，来流条件为M∞=9.59，T∞=185.6K，Tw=293.3K，Re=139436m-1。网格为nstream×nnormal=262×148。计算结果与Gnoffo[17]给出的实验和数值模拟结果进行对比。Gnoffo数值模拟时使用的网格是512×256。
图 8给出了几种格式计算得到的密度梯度等值线和流线分布。从密度梯度等值线可以看出，几种激波捕捉格式基本上都捕捉到了流场中的激波、接触间断面、超声速射流区、分离区等复杂的流动现象，但是SWENO3-PPM5和SWENO3-M分辨出的结构要更加清晰一些，尤其是在分离区及超声速射流区处；且在分离区，SWENO3-PPM5分辨出的旋涡数最多，NND格式最少。
|
| 图 8 不同格式的密度梯度等值线和流线分布 Fig. 8 Distributions of density gradient contour and streamline |
图 9和图 10比较了几种格式计算的壁面压力和热流分布。可以看出，SWENO3-PPM5的计算出的分离涡较大，且壁面压力和热流分布更加接近于Gnoffo的计算和实验结果，显示了SWENO3-PPM5较好的热流预测能力。
|
| 图 9 不同格式壁面压力分布 Fig. 9 Distribution of pressure at the surface |
|
| 图 10 不同格式壁面热流分布 Fig. 10 Distribution of heat flux at the surface |
3 结论
本文采用映射函数的思想，将一种五次分段多项式映射函数运用到中心型三阶格式上，构造了SWENO3-PPM5格式，同时与Henrick等提出的映射函数作对比(SWENO3-M格式)。Shu-Osher算例、激波/层流边界层干扰、激波/旋涡干扰以及绕25°/55°尖双锥流动的计算表明，使用映射函数的新三阶格式具有比NND、WENO3更高的分辨能力，可以有效地捕捉到流场中的激波、接触间断面等复杂的流动现象。且由于较Henrick的映射函数，分段多项式映射函数可以使最终的非线性权更加接近于格式的线性权，因此SWENO3-PPM5的分辨能力也要稍微高于SWENO3-M。
| [1] | Harten A, Enquist B, Osher S, et al. Uniformly high order essentially non-oscillatory schemes[J]. Journal of Computational Physics, 1987, 71:231–303. DOI:10.1016/0021-9991(87)90031-3 |
| [2] | Liu X, Osher S, Chan T. Weighted essentially non-oscillatory schemes[J]. Journal of Computational Physics, 1994, 115:200–212. DOI:10.1006/jcph.1994.1187 |
| [3] | Jiang G, Shu Qiwang. Efficient implementation of weighted ENO schemes[J]. Journal of Computational Physics, 1996, 126:202–228. DOI:10.1006/jcph.1996.0130 |
| [4] | 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:542–567. DOI:10.1016/j.jcp.2005.01.023 |
| [5] | Borges R, Carmona M, Costa B, et al. An improved weighted essentially non-oscillatory scheme for hyperbolilc conservation laws[J]. Journal of Computational Physics, 2008, 227:3191–3211. DOI:10.1016/j.jcp.2007.11.038 |
| [6] | Ha Y, Kim C H, Lee Y J, et al. An improved weighted essentially non-oscillatory scheme with a new smoothness indicator[J]. Journal of Computational Physics, 2013, 232:68–86. DOI:10.1016/j.jcp.2012.06.016 |
| [7] | Ping Fan, Yiqing Shen, Naolin Tian, et al. A new smoothness indicator for improving the weighted essentially non-oscillatory scheme[J]. Journal of Computational Physics, 2014, 269:329–354. DOI:10.1016/j.jcp.2014.03.032 |
| [8] | Hui Feng, Fuxing Hu, Rong Wang. A new mapped weighted essentially non-oscillatory scheme[J]. Journal of Scientific Computing, 2012, 51:449–473. DOI:10.1007/s10915-011-9518-y |
| [9] | Li Q, Liu P X, Zhang H X. Piecewise polynomial mapping method and corresponding WENO scheme with improved resolution[J]. Communication in Computational Physics, 2015, 18(5):1417–1444. DOI:10.4208/cicp.150215.250515a |
| [10] |
Zhang H X. Non-oscillatory and non-free-parameter dissipation difference scheme[J].
Acta Aerodynamica Sinica, 1988, 6(2):143–165.
(in Chinese) 张涵信. 无波动、无自由参数的耗散差分格式[J]. 空气动力学学报, 1988, 6(2) : 143–165. |
| [11] |
Zhang H X, He G H, Zhang L. Some important problems for high order accurate difference scheme solving gas dynamic equations[J].
Acta Aerodynamica Sinica, 1993, 11(4):341–356.
(in Chinese) 张涵信, 贺国宏, 张雷, 等. 高精度求解气动方程的几个问题[J]. 空气动力学学报, 1993, 11(4) : 341–356. |
| [12] |
Liu W. Nonlinear dynamics analysis for mechanism of slender wing rock and study of numerical simulation method[D]. National University of Defense Technology, 2004.(in Chinese) 刘伟.细长机翼摇滚机理的非线性动力学分析及数值模拟方法研究[D].湖南:国防科学技术大学, 2004. |
| [13] |
Li Q, Zhang H X, Gao S C. Numerical simulations on supersonic shear layer flow[J].
Acta Aerodynamica Sinica, 2000, 18(增刊):67–77.
(in Chinese) 李沁, 张涵信, 高树椿. 关于超声速剪切流动的数值模拟[J]. 空气动力学学报, 2000, 18(增刊) : 67–77. |
| [14] |
Li Q, Sun D, Zheng Y K, et al. On a class of center-type third order difference scheme orienting to engineering utilizations[J].
Acta Aerodynamica Sinica, 2013, 31(4):466–472.
(in Chinese) 李沁, 孙东, 郑永康, 等. 一类中心型三阶格式及应用[J]. 空气动力学学报, 2013, 31(4) : 466–472. |
| [15] | Sergio Pirozzoli. On the spectral properties of shock-capturing schemes[J]. Journal of Computational Physics, 2006, 219:489–497. DOI:10.1016/j.jcp.2006.07.009 |
| [16] | Inoue O, Hattori Y. Sound generation by shock vortex interactions[J]. Journal of Fluid Mechanics, 1990, 380:81–116. |
| [17] | Gnoffo P. CFD validation studies for hypersonic flow prediction[R]. AIAA 2001-1025. |