采用重组模板的权重优化WENO-Z格式

柴得林,王 强,2,易 贤,2*,刘 宇

(1. 中国空气动力研究与发展中心 结冰与防除冰重点实验室, 四川 绵阳 621000;2. 中国空气动力研究与发展中心 空气动力学国家重点实验室, 四川 绵阳 621000)

随着电子计算机技术的飞速发展,计算流体动力学在流体力学领域的理论研究与工程应用中起到越来越大的作用。作为计算流体动力学的基础之一,离散格式的性能对流场数值模拟具有重要影响。特别地,加权本质无振荡(weighted essentially non-oscillatory,WENO)格式的提出极大地推进了含激波等复杂流动结构的流场的精确数值模拟。

Liu等[1]在本质无振荡(essentially non-oscillatory, ENO)[2-3]格式的基础上创造性地提出了WENO 格式,采用子模板上的低阶格式的非线性凸组合使得格式兼具高精度与本质无振荡特性,并设计了有限体积形式的3阶、4阶精度 WENO 格式。Jiang等[4]对 WENO 格式进行理论分析,将WENO格式拓展至有限差分形式,系统地设计了任意阶有限差分形式格式的光滑因子与非线性权计算方法,他们提出的 5 阶WENO 格式成为最经典的 WENO 格式之一,一般记为WENO-JS格式。

尽管WENO-JS格式具有优越的激波捕捉性能,但仍存在耗散过大,极值点处精度降阶等不足。围绕WENO-JS格式,学者们开展了大量的性能优化研究。在非线性权计算方面,Henrick等[5]指出5阶WENO-JS格式在求解双曲守恒律时并未完全满足 5 阶精度,在极值点附近发生降阶。为此,他们在非线性权的计算中引入了一个非线性权映射函数,设计了一种完全满足 5 阶精度的改进WENO格式,记为 WENO-M 格式。根据 WENO-M 格式的思想,多种新型映射函数被提出并应用于 WENO 格式的优化[6-8]。Borges等[9]从增大间断子模板的权重分配的角度开展研究,指出增大非线性加权时间断子模板的权重,可降低格式耗散,优化格式性能;在这一理论上,他们设计了一个高阶光滑因子,构造了新的非线性权,提出了耗散更低、分辨率更高的 5 阶 WENO-Z 格式。

在WENO-Z格式的基础上,Liu等[10]改进了5阶WENO-Z格式的高阶光滑因子及其在非线性权中的应用,使得格式既满足5阶精度充分条件,又具有较低耗散。Castro等[11]给出了高阶光滑因子的通用公式,将 WENO-Z 格式拓展至任意奇数阶。Acker等[12]在WENO-Z格式的非线性权公式中增加了光滑因子比值相关项,提出WENO-Z+格式,进一步提高了格式中间断子模板的权重,改善了格式对高频波的分辨率,并指出在通常的应用中,间断子模板上的权重对WENO格式实际计算性能的影响起主要因素。文献[13-16]均引入光滑因子比值优化了格式权重, 文献[17-20]对光滑因子进行了重新设计与构造。徐维铮等[21]则对3阶WENO-Z格式的光滑因子进行了多种设计与系统研究,研究表明格式在连续解非极值点处的理论精度对实际计算性能起决定性的作用,极值点处的精度影响则较小。上述研究表明改进非线性权计算方法可有效实现WENO-Z格式耗散降低,性能提升。

与上述仅改进非线性权计算方法的研究不同,模板优化是WENO格式改进的另一重要方法。Martín等[22]和HU等[23]通过在WENO构造模板中引入下迎风模板,优化权重,分别提出了WENO-SYMBO和WENO-CU6格式;Zhu等[24]则创造性地设计了由一个5点模板和两个2点模板加权得到的5 阶有限差分 WENO格式,该格式对线性权的选择更为灵活,实现更为简单。这些研究在优化模板的基础上,对非线性权计算方法,包括线性权的选取、光滑因子的计算等,进行了一定改进,最终实现了格式性能的提升。

尽管上述研究已经实现了WENO格式性能的大幅改善,模板优化方法也往往伴随着非线性权计算方法的重新设计,但鲜有研究对非线性权计算方法改进与模板优化之间的关系进行研究,鲜有研究提出可以转化为改进非线性权计算方法的模板优化方法。

本文以5阶WENO-Z格式为研究对象,借鉴文献[22-24]的优化模板思路,在WENO-Z格式构造中引入一个由3点模板重新组合形成的4点模板,取其上重构格式为对应3点模板格式的线性组合,通过这种模板重组方法实现了格式的非线性权调节与性能提升;同时借助子模板的线性组合特性将所提模板优化方法等效转化为格式的改进非线性权计算方法;采用一系列基准问题对改进格式性能提升进行数值验证。

1 WENO-Z格式

以一维双曲守恒律为研究对象

(1)

(2)

式中,h(x)为数值通量函数,其隐式定义为

(3)

(4)

5阶WENO-JS格式或WENO-Z格式的重构模板如图1所示。

图1 5阶WENO-JS/WENO-Z格式重构模板Fig.1 Reconstruction stencils of fifth-order WENO-JS/WENO-Z

S(5)={Ii-2,Ii-1,Ii,Ii+1,Ii+2}为构造大模板,可划分为3个互有重叠的3点子模板,即S0={Ii-2,Ii-1,Ii}、S1={Ii-1,Ii,Ii+1}和S2={Ii,Ii+1,Ii+2}。在各子模板上计算数值通量对应的线性格式,可得

(5)

而大模板上数值通量对应的5阶线性格式为

(6)

式中:dm为线性权重,d0=0.1,d1=0.6,d2=0.3;m=0,1,2。引入与其相对应的非线性权ωm,则可得WENO格式的一般形式

(7)

5阶WENO-JS格式的非线性权计算公式为

(8)

式中,βm为光滑因子,表征相应模板上变量的光滑程度,ε为一正数小量,防止分母为0。文献[4]给出的k阶WENO格式的光滑因子βm通用计算公式为

(9)

对于5阶WENO-JS格式各光滑因子为

(10)

而ε则取经验值10-6。

5阶WENO-Z格式的非线性权计算公式为

(11)

式中:βm为式(10)所示的WENO-JS光滑因子;τ5=|β0-β2|为高阶光滑因子;p为度量光滑因子对非线性权的影响的指数,一般取1;εZ为一正数小量,防止分母为0。参考文献[9],本文εZ取一极小量10-40使其仅起到防止分母为0的作用。较之WENO-JS格式,WENO-Z格式增大了间断模板的非线性权,降低了格式耗散,提高了格式分辨率。

2 重组模板WENO-Z格式

2.1 模板重组

图1中WENO-Z格式计算xi+1/2处的值时所用整体模板是上迎风的,为了降低格式耗散,应使得模板更接近中心模板,因而本文在WENO格式的构造中引入中心模板S3={Ii-1,Ii,Ii+1,Ii+2},如图2所示。则改进格式含有4个子模板,记改进格式为WENO-ZF,F为单纯记号,则新模板组合下,对应加权公式为

(12)

2.2 线性权与子模板S3参数计算

由图2可知子模板S3为S1、S2的组合,假设S3上线性格式为S1、S2对应线性格式的线性组合,即

图2 5阶WENO-ZF格式重构模板Fig.2 Reconstruction stencils of fifth-order WENO-ZF

(13)

(14)

同时有

(15)

由数值通量线性格式精度式(4)得

(16)

(17)

在单元Ii=[xi-1/2,xi+1/2]上采用与式(15)相同的线性组合,则有

将式(18)代入光滑因子的定义式(9)可得S3上光滑因子计算式为

2.3 非线性加权与简化

参考WENO-Z格式的加权方法,仍将WENO-ZF格式的非线性权取为

(20)

其中:m=0,1,2,3;参考文献[9],取τ5=|β0-β2|,p=1,εZF=10-40。

(21)

式(21)中,τ5、βm、εZF与WENO-ZF格式中式(20)中的定义及取值相同。对比式(20)和式(21)可知,重组模板的模板优化方案最终等效为非线性权计算方法的优化。为减少计算量,后续数值实验中采用式(21)进行WENO-ZF格式的计算。

(22)

将其代入式(18)得

(23)

将式(22)、式(23)代入光滑因子的定义式(9),可得

(24)

因而,计算非线性权时可取

(25)

这样可以减少计算量。为便于表示,下文中式(19)和式(25)计算β3的改进格式分别记为WENO-ZF1、WENO-ZF2。

2.4 近似色散关系分析

近似色散关系(approximate dispersion relation,ADR)分析方法[25]可有效分析非线性格式的谱特性,即色散耗散特性,广泛应用于非线性格式的设计与优化。采用该方法计算所得WENO-ZF1、WENO-ZF2格式的色散耗散特性如图3所示。图中Φ(φ)为波数φ对应的修正波数,其实部和虚部分别反映格式的色散与耗散。

比较图3中各格式色散与耗散特性:低频波范围内,所有格式均有理想的色散耗散;中高频波范围内,WENO-ZF1与WENO-ZF2格式的色散与耗散均小于WENO-Z格式,而WENO-JS格式的色散和耗散则明显大于这三种格式,这表明采用模板重组引入模板S3改进WENO格式的方法可以有效降低格式色散与耗散。对比WENO-ZF1与WENO-ZF2格式的色散与耗散特性,WENO-ZF2格式除对部分中频波的色散略小于WENO-ZF1格式外,其他情况下对中高频波的色散和耗散均大于WENO-ZF1格式。

(a) 色散(a) Dispersion

(b) 耗散(b) Dissipation

3 数值实验

本节采用5阶WENO-JS、WENO-Z、WENO-ZF1和WENO-ZF2格式对一系列经典的数值算例进行模拟,以验证格式的精度及其色散耗散特性、分辨率等方面的性能。需说明的是,算例中变量如无特别说明,均为无量纲量;各算例的时间离散均采用3阶TVD Runge-Kutta 格式计算。

3.1 线性对流方程

对线性对流方程

(26)

取初始条件u(x,0)=sin(πx),在网格数为20至320的均匀网格上进行计算,结束时间取t=2,计算所得L1、L∞误差与精度如表1所示。

由表1中数据可得,四种WENO格式均可达到设计的5阶精度;相同网格数下,WENO-Z、WENO-ZF1、WENO-ZF2三种格式误差较为接近,较WENO-JS格式小一个数量级。比较WENO-Z、WENO-ZF1、WENO-ZF2三种格式误差,对于L∞误差,相同网格下,WENO-ZF1格式与WENO-ZF2格式误差相当,均小于WENO-Z格式。对于L1误差,相同粗网格下,WENO-ZF1格式与WENO-ZF2格式误差相当,均小于WENO-Z格式;相同细网格下,三种格式误差基本相同。综上,该算例中WENO-ZF1和WENO-ZF2格式可达到设计收敛精度,且有较WENO-Z格式更高的分辨率。

表1 t=2时线性对流方程各WENO格式的误差与精度

3.2 一维欧拉方程

求解欧拉方程时,为了计算的稳定性,将方程投影至特征空间,对特征变量进行求解,并采用Lax-Friedrichs通量分裂方法对通量进行分裂,采用Roe平均方法计算网格单元界面上的变量。

对一维欧拉方程

(27)

进行计算,式(27)中守恒变量U=[ρ,ρu,ρE]T,通量F=[ρu,ρu2+p,(ρE+p)u]T。

3.2.1 Sod激波管

计算Sod激波管问题,其初始条件为

(28)

左右边界设置为零梯度边界条件,取200个均匀网格进行计算,模拟该问题至t=0.4,计算结果密度分布及激波与接触间断附近的密度分布放大图如图4所示。

(a) 密度分布(a) Density distributions

(b) 接触间断附近局部放大图(b) Locally enlarged plot near the contact discontinuity

(c) 激波附近局部放大图(c) Locally enlarged plot near the shock

图4表明各格式均可无振荡计算得到激波与接触间断。但比较图4(b)和图4(c)中各格式在接触间断与激波附近的结果,WENO-ZF1、WENO-ZF2格式结果接近,略优于WENO-Z格式,而WENO-JS格式抹平较大。这表明,基于模板重组的WENO-ZF1、WENO-ZF2格式耗散低于WENO-Z格式。

3.2.2 激波-密度波相互作用问题

计算激波-密度波相互作用问题,即Shu-Osher问题,该问题的初始条件为

(29)

计算中左右边界设置为零梯度边界条件,取200个均匀网格进行计算,模拟该问题至t=1.8。由于该问题无解析解,故将在网格数为2 000的均匀网格上采用WENO-JS格式计算所得结果作为基准解。图5为网格数为200时各格式计算结果的密度分布。

图5中各格式对锯齿状的声波形成的一系列小激波及x=2.4附近的激波都有较好的捕捉效果,但对高频振荡的折射熵波计算效果较差。对比激波附近放大图,就所得激波陡峭程度而言,WENO-ZF1、WENO-ZF2格式结果接近,优于WENO-Z格式,且优于WENO-JS格式。对比折射熵波附近放大图,在对幅值的计算上,WENO-ZF2结果接近甚至优于WENO-ZF1格式,两者均优于WENO-Z格式,WENO-JS格式则最差。可见,改进格式的色散耗散特性得到了优化。

(a) 密度分布(a) Density distributions

(b) 图(a)局部放大图(b) Locally enlarged plot of graph (a)

3.3 二维欧拉方程

对二维欧拉方程

(30)

做计算,式(30)中守恒变量U=[ρ,ρu,ρv,ρE]T,各方向通量F=[ρu,ρu2+p,ρuv,(ρE+p)u]T,G=[ρv,ρuv,ρv2+p,(ρE+p)v]T。特征投影、Lax-Friedrichs通量分裂等处理方法与求解一维欧拉方程相同。

3.3.1 二维黎曼问题

该问题中取计算区域为[0,1]×[0,1],初场为

(31)

计算区域划分为400×400的均匀网格,各边界均采用零梯度边界条件,计算至t=0.8。计算结果中含一道射流和四道激波等复杂结构,可以有效校验数值格式的色散耗散特性和激波捕捉特性。各格式计算结果密度分布如图6所示,图中展示0.14到 1.70均匀分布的40条密度等值线。

由图6可得,四种WENO格式计算结果均可正确反映流场分布,实现对激波的有效捕捉。但对比滑移线附近的流场,在对复杂离散涡的计算上,WENO-ZF1和WENO-ZF2格式计算得到了更多细节结构,WENO-Z格式次之,WENO-JS更次之。这表明,四种格式中,WENO-ZF1与WENO-ZF2格式耗散最小,WENO-Z格式次之,WENO-JS格式耗散最大。

(a) WENO-JS

(b) WENO-Z

(c) WENO-ZF1

(d) WENO-ZF2

3.3.2 激波-涡相互作用

激波-涡相互作用问题描述二维空间下运动的涡穿过静止的马赫数为1.1的激波时,二者相互作用的问题。计算区域为[0,2]×[0,1]。初始条件为:在x=0.5处给定一道垂直于x轴的静止的激波,激波左侧物理量为

(32)

激波右侧物理量由Rankine-Hugoniot条件计算得出。同时,激波左侧叠加一个以(0.25,0.5)为中心的涡,即在激波左侧速度、温度和熵给定值上施加涡对应的脉动量,各脉动量为

(33)

(a) WENO-ZF1

(b) WENO-ZF2

图8 t=0.6时激波-涡相互作用y=0.5压力分布Fig.8 Pressure distributions of the shock-vortex interaction problem on central line y=0.5 at t=0.6

图7与图8中计算结果表明WENO-ZF1、WENO-ZF2格式均可较好地模拟这一问题,得到基本的流场分布,得到正确的激波结构和涡核处的低压区。比较图9中激波附近和涡核附近的压力分布局部放大图,在对激波的捕捉和涡核压力的计算上,WENO-ZF1格式与WENO-ZF2格式计算结果相当,优于WENO-Z格式,WENO-JS格式计算结果的偏差最大。综上,改进的WENO-ZF1格式、WENO-ZF2格式对激波的捕捉特性和小尺度流场结构的分辨率均优于WENO-Z格式与WENO-JS格式。

(a) 激波附近压力分布(a) Pressure distributions near the shock wave

(b) 涡核附近压力分布(b) Pressure distributions near the vortex core

表2给出采用不同WENO格式计算上述算例中线性对流方程(初始分布为正弦波,网格数为320)、二维黎曼问题、激波-涡相互作用问题的相对用时,其他算例由于耗时较短,难以精确统计,未予列出。表2中分别以各算例对应WENO-Z格式所用时间作为单位时间对其他格式所用时间进行无量纲化,其中WENO-ZF格式表示采用最初的4个子模板加权公式计算的WENO-ZF格式(该格式与WENO-ZF1格式完全等价,计算结果完全相同,故在各算例结果中没有列出)。由表中数据可得,最初的WENO-ZF格式较WENO-Z格式用时长7%～19%,WENO-ZF1格式较WENO-Z格式用时长4%～16%,WENO-ZF2格式较WENO-Z格式用时长1%～4%。可见,WENO-ZF1格式的等效转化较WENO-ZF格式提高了格式效率,WENO-ZF2格式的β3计算简化则进一步提高了格式效率。

表2 计算不同算例时各WENO格式所用相对时间

4 结论

本文在5阶有限差分WENO-Z格式的基础上,提出一种模板重组技术:将原格式中的子模板组合得到一中心子模板S3,S3上逼近函数取对应原始子模板的线性组合;在包括S3的4个子模板上采用WENO-Z格式加权方法得到WENO-ZF格式。通过数学运算将WENO-ZF格式转换为原始子模板上的改进非线性权计算方法的WENO-ZF1格式。WENO-ZF和WENO-ZF1格式的构建首次实现了格式模板优化与非线性权改进二者的等价转换。

为简化计算,提高格式效率,本文又将WENO-ZF1中S3上光滑因子β3简化为原始子模板光滑因子的线性组合,得到了计算更为简洁的WENO-ZF2格式。

ADR分析表明,WENO-ZF1格式和WENO-ZF2格式色散耗散特性得到不同程度的优化。数值实验表明,WENO-ZF1格式和WENO-ZF2格式可达到理想的设计精度;较WENO-Z格式、WENO-JS格式,具有更低的耗散,对激波有更优的捕捉特性,对小尺度流场结构有更高的分辨率。格式效率上,WENO-ZF2格式效率高于WENO-ZF1格式,又高于采用4个子模板加权公式计算的最初WENO-ZF格式。