高精度WCNS格式加权策略改进与数值验证

杨 强, 李伟鹏

( 上海交通大学 航空航天学院,上海 200240)

可压缩流动广泛存在于超声速/高超声速飞行器中,激波间断和湍流共存且相互耦合作用给高可靠性数值模拟带来了巨大挑战.近40年,人们发展了众多高精度激波捕获方法,例如总变差不增(TVD)方法[1]、ASUM[2]、基本无振荡/加权基本无振荡(ENO/WENO)方法[3-5]、加权非线性紧致格式(WCNS)方法[6]、基于目标本质无振荡(TENO)方法[7]和AFTENO[8]等.这些格式在激波间断附近增加局部耗散来抑制非物理数值振荡,而在光滑湍流区域减少耗散以保持对多尺度湍流结构的分辨率.经典激波捕捉格式,如TVD方法和ENO方法,虽然可以解决激波间断无数值振荡,但对湍流光滑区域产生了过多的数值耗散.WENO方法在ENO方法的基础上,引入子模板光滑因子,通过非线性加权方法来动态调整局部耗散[5],使得高精度数值格式得到进一步发展.Borges等[9]在非线性权值中引入更高阶的光滑指标,提出计算成本更低、精度更高的WENO-Z格式.虽然WENO格式的精度可以设计为任意高阶[10-13],但WENO格式对光谱区域的分辨率仍然不能令人满意[14],它的加权机制和光滑度指标设定依旧对小尺度湍流结构中引入过多的数值耗散[15].1999年,Deng等[16]提出非线性紧致格式(CNS),后续引入加权技术,构造了 WCNS[6].近年来,Fu等[7]提出了另一类非线性加权框架,即TENO 方法,通过类ENO方法的模板选择来控制非线性耗散,在无间断区域(包括光滑临界点)保证最优的高阶空间重构,在间断区域通过移除跨越间断的模板来避免数值振荡,这种方法有利于候选模板的使用和高阶格式的改进.当子模板中存在多个相邻的不连续间断,该方法可退化为三阶重构.即使对于高阶版本,数值鲁棒性也得到了明显改善,且避免了阶跃退化问题.Peng等[17]在此基础上进一步改进了TENO格式的CT自适应函数策略,发展了TENO5-LAD格式,实验结果证明新方法能够更好地抑制不连续点附近的数值振荡,并能以较低的额外计算成本保持TENO的低耗散特性.而Ye等[8]结合AFWENO格式和TENO格式的优点构建了可以实现任意高阶精度的AFTENO格式,在激波捕获和波数分辨率方面,相较于AFWENO-JS和AFWENO-Z有着明显的优势.之后,在加权紧致格式方面,Zhang等[18]将WCNS与TENO格式组合发展了TCNS,实验结果显示TCNS能够在模拟中捕获更丰富的波结构,特别是在对瑞利-泰勒不稳定性问题的模拟中,TCNS在较粗网格上取得了与WCNS-JS[5]相似甚至更好的结果.

将TENO格式的模板选择和加权策略引入到WCNS的构造过程,发展了WCNS-T格式,并在此基础上测验7阶精度下的数值表现情况.新格式在扩展模板宽度的同时改变了模板选择过程,当全局模板存在多个间断时,WCNS7-T格式将不会使用所有候选子模板,通过引入截断函数,将光滑和非光滑候选模板完全分离,并将非光滑候选模板的权重直接设为0,彻底剔除非光滑模板的影响.该模板策略能使格式在高阶精度时相较原来能保持更高的数值稳健性.同时,光滑度因子的新评测策略较之前能更准确判定子模板的光滑程度,降低了TENO格式中复杂光滑指标所带来的不必要的计算量,提高了格式捕捉激波的能力.一维和二维数值算例结果表明,WCNS7-T格式相比于传统的WCNS7-Z格式具有更高的精度和分辨率.

1 数值方法

1.1 控制方程

针对黎曼问题,差分格式的推导仅考虑一维守恒型Euler方程:

∂Q/∂t+∂E/∂x=0

(1)

守恒变量Q和E为

(2)

式中:ρ为密度;u为来流速度;p为压力;γ为绝热指数.

在每个网格节点xi处,对式(1)进行有限差分格式离散得:

(3)

(4)

1.2 WCNS7-Z格式构造

WCNS是在CNS基础上构造的,在解决CNS三对角线反演效率低于显式格式的问题同时,使用类似WENO格式加权的思想,充分利用计算模板的信息,对光滑模板和包含间断的模板分配对应的权值系数,从而使格式更高效.WCNS的构造方法为3个步骤:① 特征变量节点到半节点的非线性加权插值;② 半节点的通量求解;③ 节点-半节点到节点的混合差分求导.重点介绍第1步骤中的非线性加权插值过程.

与5阶格式的模板不同,在7阶WCNS[19]的构造中使用4个4点子模板:(xi-3,xi-2,xi-1,xi),(xi-2,xi-1,xi,xi+1),(xi,xi+1,xi+2,xi+3),(xi,xi+1,xi+2,xi+3),如图1所示.图中:L为全局模板;K1、K2、K3、K4分别是4个子模板.每个子模板上的4阶紧致格式设定为特定偏置形式:

图1 WCNS7-Z格式的模板构造Fig.1 Candidate stencils of WCNS7-Z scheme

(5)

设定子模板的线性权值为

(6)

组合可得到7阶紧致插值格式:

(7)

经过非线性权w0、w1、w2、w3再组合后可得:

(8)

对于非线性权wk的使用策略,Jiang等[5]设计的经典权函数方案对去除跨越间断模板带来的振荡现象有很好的效果,但也引入了较大的非线性误差,且在极值点附近格式的精度不能得到保证,Borges等[9]引入更高阶的光滑指标,并因此提出了计算成本更低、精度更高的Z型权函数,显著地降低了由非线性权带来的格式耗散.权重方案采用形式简洁的Z型权函数:

(9)

式中:k=0,1,2,3;ε=10-12为防止分母出现零除引入的小量;τ7为全局模板光滑度指标,

τ7=|β0+3β1-3β2-β3|=O(h7)

(10)

式中:O(h7)是数学上表示“与h7同阶的项”的一种记法,并且可明显地表示出截断误差的量级.

非线性权重中局部平滑度指标βk的大小决定了每一个子模板的光滑度,通式为

(11)

对于k=0的子模板:

(12)

对于k=1的子模板:

(13)

对于k=2的子模板:

(14)

对于k=3的子模板:

(15)

1.3 WCNS7-T格式构造

不同于传统的权函数方案策略,TENO格式的模板加权策略在优化WENO格式模板选择的同时,更为准确地评估了模板光滑度,起到了间断检测器的作用,并在重构中以完全去除振荡模板来取代WENO格式通过非线性加权策略减小振荡模板影响的方式,使格式具有更加健壮的间断捕捉能力.在高阶的TENO方案中,如Hiejima[23]所采取的策略一样,使用了递增宽度的模板序列.图2展示了WCNS7-T格式的模板构造,为了获得7阶空间精度,使用了3个3点模板和2个4点模板,从而使模板的选择更加紧致.

图2 WCNS7-T格式的模板构造Fig.2 Candidate stencils of WCNS7-T scheme

传统WCNS中插值过程计算量大且复杂,对光滑区域来说这种插值没有必要.利用评估策略发现间断点,对于格式计算效率和精度都有非常重要的意义.将TENO格式的模板加权和宽度递增思想引进WCNS构造,在构造间断检测器的同时,更快地实现线性格式和非线性格式之间的重建切换.因此,在WCNS7-T格式中将表征插值模板光滑度的权重分配方案更替为TENO的光滑度评测策略,有利于增强格式的激波捕捉能力.优化后WCNS7-T格式的权重方案为

(16)

式中:k=0, 1, …,K-3;γk、χk、δk、q皆为光滑因子指示参数,经截断函数CT过滤得到最终非线性权值wk.

光滑因子计算式与WCNS的不同,具体形式为

对于非光滑的模板,根据类ENO模板的选择,δk=0,这样振荡模板将会被剔除;而对于光滑模板,δk=1,格式将恢复为线性方案.显然,截断参数CT的取值至关重要,而控制耗散特性的CT取值依赖于不同的问题,且在实际应用中,参数调整通常非常耗时.Fu等[7]给出了CT的自适应函数取值公式.对于本文中的算例,控制耗散特性的CT值在几组实验下对比测算后,选定为10-7,其中K阶精度τK[7]对应各个阶数的计算量过大,尤其在高阶时计算成本过高,本文将其值依旧限定为1,ε=10-12,其余光滑度指标参数与TENO给定的值保持一致.

2 数值测试算例

为了验证WCNS7-T格式计算表现和性能,开展了一维和二维问题的算例验证,计算结果与传统的WCNS7-Z格式进行对比.求解过程中采用了HLLC通量分裂方法[24],时间推进采用3阶TVD Runge-Kutta格式[25],实验的参考精确解由WCNS5-JS通过网格加密获得.

2.1 数值精度测试

为验证WCNS7-T格式的计算性能,首先以文献[8]中的数值精度实验为算例,给定初始条件为

(17)

计算域[x,y]同样设定为[0, π]×[0, π],u、v分别为x、y方向上的速度分量.各个方向上使用周期边界条件,精确的密度解为ρ(x,y,t)=1+0.2sin(x+y-2t),计算时间直到t=2 s,表1中给出了WCNS7-Z和WCNS7-T格式L1和L∞的误差(E1和E2)和精度(P1和P2)比较,网格单元划分Nx从102到3202.可知,相较于WCNS7-Z格式,WCNS7-T格式体现出更好的网格收敛性,且格式的设计精度也得到了较好的数值验证.

表1 二维欧拉方程的WCNS7-Z 和 WCNS7-T 格式的误差和精度Tab.1 Errors and accuracy with two dimensional Euler equations for WCNS7-Z and WCNS7-T schemes

2.2 Sod问题

一维Sod激波管问题是检验差分格式性质的经典算例之一.该问题包含一道激波、一个接触间断和一个膨胀波,问题的具体描述参考文献[26].给定初始条件为

计算域为[-0.5, 0.5],网格均匀划分为101个网格点,其中Δx=0.01.在边界处设置零梯度边界条件,以收敛条件判断数CFL为0.4和定时间步长Δt=0.01 s,计算到t=0.2 s.计算结果如图3所示,WCNS7-Z格式和WCNS7-T格式都能很好地捕捉间断,不会产生虚假的数值振荡,表明在间断处两种格式都引入了足够且稳定的耗散.然而,WCNS7-T格式要比WCNS7-Z格式在激波和间断处斜率更陡峭,对激波间断位置的分辨率更准确.计算结果表明,WCNS7-T格式引入的数值耗散相对较小,对激波捕获的精度更高.

图3 Sod激波管问题的的密度分布结果Fig.3 Density distribution results for the Sod problem

2.3 Lax问题

Lax问题是左稀疏波、右激波类型的黎曼问题,同样能够用于检验数值算法捕捉间断能力.设定初始条件为

(ρ,u,p)=

计算域为[0, 1],均匀划分为101个网格点,其中Δx= 0.01,CFL为0.4,Δt=0.01 s,计算到t=0.14 s,密度分布计算结果如图4所示.由放大视图可见,数值解未观察到过冲或振荡问题,与精确解较为吻合.WCNS7-T通过加权组合多个低阶紧致格式达到了高阶紧致格式,在间断附近采用TENO格式完全剔除振荡模板的方法,获得了比WCNS7-Z更加陡峭的间断捕获,表明WCNS7-T格式具有更弱的局部耗散特性使其能够更高效地捕捉间断.

图4 Lax问题的密度分布结果Fig.4 Density distribution results for the Lax problem

2.4 Shu-Osher问题

Shu-Osher问题涉及向右移动的超音速激波与正弦分布密度场相互作用, 构造一个既有平滑结构又有不连续结构的流场,用于测试格式捕获流场间断和解析平滑多尺度结构的能力.设定初始条件为

(ρ,u,p)=

接触间断由两个强激波的碰撞产生,以 4 001 网格点计算结果为参考,终止时间为t=1.8 s,两种格式计算网格取为400,CFL数设为0.5,测试结果如图5所示.从WCNS7-Z和WCNS7-T格式的计算对比结果可见,WCNS7-T格式的分辨率更高,对高频波动的求解有着更好的波形和幅度,对可分辨频域的波峰波谷区域抹平也更少,表明WCNS7-T格式比经典的WCNS7-Z格式具有更优的间断捕捉能力和脉动分辨率.

图5 Shu-Osher问题的密度分布结果Fig.5 Density distribution results for the Shu-Osher problem

2.5 Titarev-Toro问题

Titarev-Toro问题是对Shu-Osher问题的强扩展,初始条件给定为

(ρ,u,p)=

计算域为[-5, 5],均分为1 001个网格点,其中Δx=0.01,边界设置为零梯度边界条件.给定时间步长Δt=0.5 s,CFL数设为0.5,终止时间为t=5 s.由于没有理论上的精确解,所以参考结果由WCNS5-JS格式在 5 001 个网格点计算获得.图6给出了两种格式计算的密度分布, 与Shu-Osher问题测试算例相似,WCNS7-T格式在激波后解析出了更好的密度波形,而WCNS7-Z格式结果中小尺度波结构较为模糊,数值耗散也较强,分辨率也较低,表明WCNS7-T格式比WCNS7-Z格式能更优、更好地解析局部大振幅密度波动问题.

图6 Titarev-Toro问题的密度分布结果Fig.6 Density distribution results for the Titarev-Toro problem

2.6 二维黎曼问题

Lax等[27]使用的19个二维黎曼问题中的第3个问题和第6个问题所产生的流场都具有丰富的小尺度特征,非常适合用于测验数值格式求解精细流场结构的能力,因此选择其中的这两个问题来评估目前WCNS7-T格式的性能.

2.6.1算例3初始条件给定为

(ρ,u,v,p)=

设定计算域为[0, 1]×[0, 1],网格尺寸划分为 1 001×1 001,终止时间为t=0.3 s,边界条件都设定为无反射边界条件.Kelvin-Helmholtz不稳定性在滑移线上产生的小尺度复杂结构情况可用来评估数值格式的耗散表现.图7所示为黎曼问题3的密度轮廓结果.从图7可以看出,对比于WCNS7-Z格式,WCNS7-T格式在滑移线上捕获的涡量更多,能够更好地捕捉到激波间断导致的流动不稳定性.此外,在这种数值情况下,WCNS7-T的结果与文献[28-30]中一样在引用TENO非线性加权后出现了对称性不足的现象,本文推测可能是控制格式耗散表现的CT值的取值问题所致.

图7 黎曼问题3密度轮廓结果Fig.7 Density contour results for Riemann problem 3

2.6.2算例6第2个黎曼的计算域、网格划分和边界条件不变,终止时间设定为t=0.3 s.初始条件为

图8所示为黎曼问题6的密度轮廓结果.由图可见,两种格式都解析出了接触线上丰富的小尺度流动结构,而WCNS7-T格式的结果有更丰富的小尺度结构,流场分辨率更高,产生的耗散也小于传统的WCNS7-Z格式.表明WCNS7-T格式经过TENO非线性加权后数值耗散相比之前有所减小,且模板选择策略的改变也进一步降低了整体的数值耗散,使格式在流场计算时能更准确模拟流场细节结构.

图8 黎曼问题6密度轮廓结果Fig.8 Density contour results for Riemann problem 6

3 结论

基于TENO格式的模板选择和加权策略,将其间断检测和光滑度因子引入WCNS构造中,发展了一种7阶精度的WCNS7-T格式,在抑制数值耗散的同时增强了激波捕捉能力,通过一维和二维问题开展算例测试,验证了WCNS7-T格式的良好性能,重要结论总结如下:

(1) WCNS7-T格式体现出了流场光滑区域耗散小、激波间断区域分辨率高、计算结果更准确的特点,表明经权重优化后格式不仅保持了原WCNS的稳定,也强化了小扰动分辨率和激波捕获能力.

(2) WCNS7-T格式以完全剔除振荡模板来取代传统的加权分配,而在光滑区域恢复最优非线性插值,不仅降低了格式在间断处需要更复杂的平滑指标所带来的较高计算量,保证了格式的数值精度,也提高了格式的间断检测能力,增强了加权方法抑制数值振荡的能力.