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

柴得林; 王强; 易贤; 刘宇; CHAI Delin; WANG Qiang; YI Xian; LIU Yu

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

doi: 10.11887/j.cn.202401020

柴得林¹ ，王强^1,2 ，易贤^1,2 ，刘宇¹

1. 中国空气动力研究与发展中心结冰与防除冰重点实验室，四川绵阳 621000

2. 中国空气动力研究与发展中心空气动力学国家重点实验室，四川绵阳 621000

基金项目: 国家科技重大专项资助项目（J2019-Ⅲ-0010-0054）；国家自然科学基金资助项目（12172372）

详细信息

作者简介

柴得林（1991—），男，山西临汾人，工程师，博士，E-mail: chaidelin@cardc.cn

通讯作者

易贤（1977—），男，四川金堂人，研究员，博士，博士生导师，E-mail: yixian@cardc.cn

中图分类号: V211；O241

文献标识码: A

文章编号: 1001-2486(2024)01-187-11

Weight-optimized WENO-Z schemes with reorganized stencil

CHAI Delin¹ ， WANG Qiang^1,2 ， YI Xian^1,2 ， LIU Yu¹

1. Key Laboratory of Icing and Anti/De-Icing， China Aerodynamics Research and Development Center, Mianyang 621000 , China

2. State Key Laboratory of Aerodynamics， China Aerodynamics Research and Development Center, Mianyang 621000 , China

摘要

针对精确模拟含激波等复杂流动结构的流场对高精度格式的低耗散低色散要求，基于5阶有限差分WENO-Z格式，提出一种模板重组技术。在计算WENO非线性权时，引入一个由3点模板重新组合的4点模板，优化原格式中各模板的权重分配，进而提出了两种改进WENO-Z格式。采用近似色散关系分析方法对改进前后格式色散与耗散特性进行了对比与分析。分析表明：两种改进格式耗散有不同程度的降低。数值实验表明：改进格式具有更优越的激波捕捉性能，对小尺度流场结构具有更高的分辨率。

关键词

模板重组 / 权重 / WENO格式

Abstract

For developing high-order scheme with low dissipation and low dispersion to accurately simulate the flow field with complex structures such as shock wave, a method of stencil reorganization was proposed based on the fifth-order finite difference WENO-Z scheme. A four-point central stencil recombined by three-point stencils is introduced in the calculation of WENO nonlinear weight in order to optimize the weight allocation of each template in the original format, and two improved WENO-Z schemes were proposed. The dispersion and dissipation properties of the schemes were compared and analyzed via the approximate dispersion relation analysis, which shows that the dissipation of the two improved schemes decreases at different extent. The numerical experiments show that the improved schemes have better shock capture property and higher resolution for small-scale flow structures.

Keywords

stencil reorganization / weight / WENO scheme

1 WENO-Z格式 2 重组模板WENO-Z格式 2.1 模板重组 2.2 线性权d^ZF_m与子模板S₃参数计算 2.3 非线性加权与简化 2.4 近似色散关系分析 3 数值实验 3.1 线性对流方程 3.2 一维欧拉方程 3.2.1 Sod激波管 3.2.2 激波-密度波相互作用问题 3.3 二维欧拉方程 3.3.1 二维黎曼问题 3.3.2 激波-涡相互作用 4 结论

随着电子计算机技术的飞速发展，计算流体动力学在流体力学领域的理论研究与工程应用中起到越来越大的作用。作为计算流体动力学的基础之一，离散格式的性能对流场数值模拟具有重要影响。特别地，加权本质无振荡（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格式

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

\frac{\partial u}{\partial t} + \frac{\partial f}{\partial x} = 0 x \in [a, b]

(1)

式中，u为守恒变量，f=f（u）为通量函数。将计算区域[a，b]划分为N个等间距网格，各网格单元长度为Δx=x_i_+1/2-x_i_-1/2，且a=x_1/2，b=x_N_+1/2，网格中心

x_{i} = \frac{1}{2} (x_{i + 1 / 2} + x_{i - 1 / 2})

，网格单元I_i=[x_i_-1/2，x_i_+1/2]。在上述网格上进行空间离散，可得式（1）的半离散形式

\frac{d u_{i}}{d t} = - \frac{h (x_{i + 1 / 2}) - h (x_{i - 1 / 2})}{Δ x}

(2)

式中，h（x）为数值通量函数，其隐式定义为

f (x) = \frac{1}{Δ x} \int_{x - Δ x / 2}^{x + Δ x / 2} h (ξ) d ξ

(3)

根据式（3），则可在给定模板S_i={I_i-r，I_i-r₊₁，···，I_i，···，I_i₊_s}上得到h（x）的唯一的r+s次多项式逼近函数

\hat{f} （ x ）

，该函数满足

\hat{f} (x) = h (x) + O (Δ x^{r + s + 1})

(4)

记

\hat{f} （ x ）

在x_i_+1/2处的值为

{\hat{f}}_{i + 1 / 2}

，则

{\hat{f}}_{i + 1 / 2}

即为数值通量对应的线性格式。

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

图15阶WENO-JS/WENO-Z格式重构模板

Fig.1Reconstruction stencils of fifth-order WENO-JS/WENO-Z

S^（5）={I_i_-2，I_i_-1，I_i，I_i₊₁，I_i₊₂}为构造大模板，可划分为3个互有重叠的3点子模板，即S₀={I_i_-2，I_i_-1，I_i}、S₁={I_i_-1，I_i，I_i₊₁}和S₂={I_i，I_i₊₁，I_i₊₂}。在各子模板上计算数值通量对应的线性格式，可得

\{\begin{matrix} {\hat{f}}_{i + 1 / 2}^{0} = \frac{1}{6} (2 f_{i - 2} - 7 f_{i - 1} + 11 f_{i}) \\ {\hat{f}}_{i + 1 / 2}^{1} = \frac{1}{6} (- f_{i - 1} + 5 f_{i} + 2 f_{i + 1}) \\ {\hat{f}}_{i + 1 / 2}^{2} = \frac{1}{6} (2 f_{i} + 5 f_{i + 1} - f_{i + 2}) \end{matrix}

(5)

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

\begin{matrix} {\hat{f}}_{i + 1 / 2}^{linear} = \frac{1}{60} (2 f_{i - 2} - 13 f_{i - 1} + 47 f_{i} + 27 f_{i + 1} - 3 f_{i + 2}) \\ = \sum_{m = 0}^{2} d_{m} {\hat{f}}_{i + 1 / 2}^{m} \end{matrix}

(6)

式中：d_m为线性权重，d₀=0.1，d₁=0.6，d₂=0.3；m=0，1，2。引入与其相对应的非线性权ω_m，则可得WENO格式的一般形式

{\hat{f}}_{i + 1 / 2} = \sum_{m = 0}^{2} ω_{m} {\hat{f}}_{i + 1 / 2}^{m} m = 0,1, 2

(7)

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

\{\begin{matrix} ω_{m} = \frac{α_{m}}{\sum_{l = 0}^{2} α_{l}} \\ α_{m} = \frac{d_{m}}{{(β_{m} + ε)}^{2}} \end{matrix}

(8)

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

β_{m} = \sum_{l = 1}^{k - 1} Δ x_{i}^{2 l - 1} \int_{x_{i - 1 / 2}}^{x_{i + 1 / 2}} {(\frac{d^{l}}{d x^{l}} {\hat{f}}_{i}^{m} (x))}^{2} d x

(9)

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

\{\begin{matrix} β_{0} = \frac{13}{12} {(f_{i - 2} - 2 f_{i - 1} + f_{i})}^{2} + \frac{1}{4} {(f_{i - 2} - 4 f_{i - 1} + 3 f_{i})}^{2} \\ β_{1} = \frac{13}{12} {(f_{i - 1} - 2 f_{i} + f_{i + 1})}^{2} + \frac{1}{4} {(f_{i - 1} - f_{i + 1})}^{2} \\ β_{2} = \frac{13}{12} {(f_{i} - 2 f_{i + 1} + f_{i + 2})}^{2} + \frac{1}{4} {(3 f_{i} - 4 f_{i + 1} + f_{i + 2})}^{2} \end{matrix}

(10)

而ε则取经验值10^-6。

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

\{\begin{matrix} ω_{m}^{Z} = \frac{α_{m}^{Z}}{\sum_{l = 0}^{2} α_{l}^{Z}} \\ α_{m}^{Z} = d_{m} [1 + {(\frac{τ_{5}}{β_{m} + ε^{Z}})}^{p}] \end{matrix}

(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格式计算x_i_+1/2处的值时所用整体模板是上迎风的，为了降低格式耗散，应使得模板更接近中心模板，因而本文在WENO格式的构造中引入中心模板S₃={I_i_-1，I_i，I_i₊₁，I_i₊₂}，如图2所示。则改进格式含有4个子模板，记改进格式为WENO-ZF，F为单纯记号，则新模板组合下，对应加权公式为

(12)

图25阶WENO-ZF格式重构模板

Fig.2Reconstruction stencils of fifth-order WENO-ZF

2.2 线性权d^ZF_m与子模板S₃参数计算

由图2可知子模板S₃为S₁、S₂的组合，假设S₃上线性格式为S₁、S₂对应线性格式的线性组合，即

(13)

由WENO-Z格式的非线性权计算式（11）易知，WENO-Z格式中S₀、S₁、S₂三个子模板上非线性权的差异受各模板线性权和光滑因子两个因素影响。假设在WENO-ZF格式中子模板S₀、S₁、S₂上线性权完全相同，即

d_{0}^{Z F} = d_{1}^{Z F} = d_{2}^{Z F}

，这使得其非线性权完全表示原格式中模板光滑因子的差异对非线性权的影响；而4点子模板S₃的非线性权则主要表示原格式中子模板上线性权的不同对非线性权的影响。此外，为了计算稳定，改进格式应为各子模板对应格式的凸组合，即各线性权应为非负数。综上根据改进前后子模板的线性权对应关系，易得

\{\begin{matrix} d_{0}^{Z F} = d_{1}^{Z F} = d_{2}^{Z F} = 0.1 \\ d_{3}^{Z F} = 0.7 \end{matrix}

(14)

同时有

\begin{matrix} {\hat{f}}_{i + 1 / 2}^{3} = \frac{5}{7} {\hat{f}}_{i + 1 / 2}^{1} + \frac{2}{7} {\hat{f}}_{i + 1 / 2}^{2} \\ = \frac{1}{42} (- 5 f_{i} + 29 f_{i + 1} + 20 f_{i + 1} - 2 f_{i + 2}) \end{matrix}

(15)

由数值通量线性格式精度式（4）得

(16)

将式（16）代入

{\hat{f}}_{i + 1 / 2}^{3}

的定义式，即式（13）或式（15），可得

{\hat{f}}_{i + 1 / 2}^{3} = h_{i + 1 / 2} + O (Δ x^{3})

(17)

即

{\hat{f}}_{i + 1 / 2}^{3}

为一个3阶格式。需说明的是，调整模板S₃上通量格式的线性组合式（13）中的系数α，可使得通量格式最高达到4阶，形成4阶中心格式。为了一定程度上分离子模板上线性权和光滑因子对非线性权的影响，基于WENO-ZF格式中子模板S₀、S₁、S₂上线性权

d_{0}^{Z F} = d_{1}^{Z F} = d_{2}^{Z F}

的假设，本文取α=5/7得到3阶精度通量格式，如式（15）所示。

在单元I_i=[x_i_-1/2，x_i_+1/2]上采用与式（15）相同的线性组合，则有

\begin{matrix} {\hat{f}}^{3} (x) = \frac{5}{7} {\hat{f}}^{1} (x) + \frac{2}{7} {\hat{f}}^{2} (x) \\ = \frac{- 5 f_{i - 1} + 176 f_{i} - f_{i + 1} - 2 f_{i + 2}}{168} + \\ \frac{- 5 f_{i - 1} - 6 f_{i} + 13 f_{i + 1} - 2 f_{i + 2}}{14 Δ x} x + \\ \frac{5 f_{i - 1} - 8 f_{i} + f_{i + 1} + 2 f_{i + 2}}{14 Δ x^{2}} x^{2} \end{matrix}

(18)

将式（18）代入光滑因子的定义式（9）可得S₃上光滑因子计算式为

\begin{matrix} β_{3} = \frac{13}{3} {(\frac{5 f_{i - 1} - 8 f_{i} + f_{i + 1} + 2 f_{i + 2}}{14})}^{2} + \\ {(\frac{- 5 f_{i - 1} - 6 f_{i} + 13 f_{i + 1} - 2 f_{i + 2}}{14})}^{2} \end{matrix}

(19)

2.3 非线性加权与简化

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

\{\begin{matrix} ω_{m}^{Z F} = \frac{α_{m}^{Z F}}{\sum_{l = 0}^{3} α_{l}^{Z F}} \\ α_{m}^{Z F} = d_{m}^{Z F} [1 + {(\frac{τ_{5}}{β_{m} + ε^{Z F}})}^{p}] \end{matrix}

(20)

其中：m=0，1，2，3；参考文献^[9]，取

τ_{5} = |β_{0} - β_{2}| ， p = 1 ， ε^{Z F} = 10^{- 40}

。

子模板线性格式计算式（5）和式（15）、WENO加权计算式（12）、子模板光滑因子计算式（10）和式（19）和非线性权计算式（20）组成5阶WENO-ZF格式的整个计算过程。显然，当式（12）中非线性权ω^ZF_m取为线性权d^ZF_m时，WENO-ZF格式变为5阶线性格式，具有5阶精度。

对比WENO-ZF格式与WENO-Z格式计算过程，前者引入更多计算环节，尤其是权重β₃的计算。为了保证格式效率，下面对WENO-ZF格式做一定的简化。因

{\hat{f}}_{i + 1 / 2}^{3}

为为

{\hat{f}}_{i + 1 / 2}^{1}

为、

{\hat{f}}_{i + 1 / 2}^{2}

为的线性组合，则WENO-ZF经过化简仍可表示为3个子模板S₀、S₁、S₂上加权所得函数，化简过程省略，WENO-ZF格式等效计算公式为

\{\begin{matrix} α_{0}^{Z F 1} = d_{0} (1 + \frac{τ_{5}}{β_{0} + ε^{Z F}}) \\ α_{1}^{Z F 1} = d_{1} (1 + \frac{1}{6} \frac{τ_{5}}{β_{1} + ε^{Z F}} + \frac{5}{6} \frac{τ_{5}}{β_{3} + ε^{Z F}}) \\ α_{2}^{Z F 1} = d_{2} (1 + \frac{1}{3} \frac{τ_{5}}{β_{2} + ε^{Z F}} + \frac{2}{3} \frac{τ_{5}}{β_{3} + ε^{Z F}}) \\ ω_{m}^{Z F 1} = \frac{α_{m}^{Z F 1}}{\sum_{l = 0}^{2} α_{l}^{Z F 1}} \\ {\hat{f}}_{i + 1 / 2}^{Z F 1} = \sum_{m = 0}^{2} ω_{m}^{Z F 1} {\hat{f}}_{i + 1 / 2}^{m} m = 0,1, 2 \end{matrix}

(21)

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

接下来尝试对β₃进行简化计算，S₁、S₂上逼近多项式函数

{\hat{f}}^{1} （ x ）

、

{\hat{f}}^{2} （ x ）

可表示为

\{\begin{matrix} {\hat{f}}^{1} (x) = b_{10} + b_{11} x + b_{12} x^{2} \\ {\hat{f}}^{2} (x) = b_{20} + b_{21} x + b_{22} x^{2} \end{matrix}

(22)

将其代入式（18）得

\begin{matrix} {\hat{f}}^{3} (x) = (\frac{5}{7} b_{10} + \frac{2}{7} b_{20}) + (\frac{5}{7} b_{11} + \frac{2}{7} b_{21}) x + \\ (\frac{5}{7} b_{12} + \frac{2}{7} b_{22}) x^{2} \end{matrix}

(23)

将式（22）、式（23）代入光滑因子的定义式（9），可得

\begin{matrix} β_{3} = {(\frac{5}{7} b_{11} + \frac{2}{7} b_{21})}^{2} Δ x^{2} + \frac{13}{3} {(\frac{5}{7} b_{12} + \frac{2}{7} b_{22})}^{2} Δ x^{4} \\ = \frac{5}{7} (b_{11}^{2} Δ x^{2} + \frac{13}{3} b_{12}^{2} Δ x^{4}) - \frac{10}{49} {(b_{11} - b_{21})}^{2} Δ x^{2} + \\ \frac{2}{7} (b_{21}^{2} Δ x^{2} + \frac{13}{3} b_{22}^{2} Δ x^{4}) - \frac{130}{147} {(b_{12} - b_{22})}^{2} Δ x^{4} \\ = \frac{5}{7} β_{1} + \frac{2}{7} β_{2} + k_{1} f^{''' 2} Δ x^{6} + k_{2} f^{''' 2} Δ x^{6} + O (Δ x^{7}) \\ = \frac{5}{7} β_{1} + \frac{2}{7} β_{2} + O (Δ x^{6}) \end{matrix}

(24)

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

β_{3} = \frac{5}{7} β_{1} + \frac{2}{7} β_{2}

(25)

这样可以减少计算量。为便于表示，下文中式（19）和式（25）计算β₃的改进格式分别记为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格式的色散和耗散则明显大于这三种格式，这表明采用模板重组引入模板S₃改进WENO格式的方法可以有效降低格式色散与耗散。对比WENO-ZF1与WENO-ZF2格式的色散与耗散特性，WENO-ZF2格式除对部分中频波的色散略小于WENO-ZF1格式外，其他情况下对中高频波的色散和耗散均大于WENO-ZF1格式。

图3各格式ADR分析色散耗散特性

Fig.3ADR dispersion and dissipation of the schemes

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，计算所得L₁、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格式。对于L₁误差，相同粗网格下，WENO-ZF1格式与WENO-ZF2格式误差相当，均小于WENO-Z格式；相同细网格下，三种格式误差基本相同。综上，该算例中WENO-ZF1和WENO-ZF2格式可达到设计收敛精度，且有较WENO-Z格式更高的分辨率。

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

Tab.1 Difference WENO scheme error and order of the linear transport equation at t=2

3.2 一维欧拉方程

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

对一维欧拉方程

\frac{\partial U}{\partial t} + \frac{\partial F}{\partial x} = 0

(27)

进行计算，式（27）中守恒变量U=[ρ，ρu，ρE]^T，通量F=[ρu，ρu²+p，（ρE+p）u]^T。

3.2.1 Sod激波管

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

(28)

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

图4t=0.4时Sod激波管密度分布

Fig.4Density distributions of the Sod shock tube at t=0.4

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

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

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

(ρ, u, p) = \{\begin{matrix} (3.857143,2.629369,10.333333) \\ - 5 ⩽ x < - 4 \\ (1 + 0.2 s i n (5 x), 0,1.0) \\ - 4 ⩽ x ⩽ 5 \end{matrix}

(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格式则最差。可见，改进格式的色散耗散特性得到了优化。

图5t=1.8时激波-密度波相互作用密度分布

Fig.5Density distributions of the shock-density wave interaction problem at t=1.8

3.3 二维欧拉方程

对二维欧拉方程

\frac{\partial U}{\partial t} + \frac{\partial F}{\partial x} + \frac{\partial G}{\partial y} = 0

(30)

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

3.3.1 二维黎曼问题

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

(ρ, u, v, p) = \{\begin{matrix} (1.5,0, 0,1.5) \\ 0.8 ⩽ x ⩽ 1,0.8 ⩽ y ⩽ 1 \\ (0.5323,1.206,0, 0.3) \\ 0 ⩽ x < 0.8,0.8 ⩽ y ⩽ 1 \\ (0.138,1.206,1.206,0.029) \\ 0 ⩽ x < 0.8,0 ⩽ y < 0.8 \\ (0.5323,0, 1.206,0.3) \\ 0.8 ⩽ x ⩽ 1,0 ⩽ y < 0.8 \end{matrix}

(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格式耗散最大。

图6t=0.8时二维黎曼问题密度分布

Fig.6Density distributions of 2D Riemann problem at t=0.8

3.3.2 激波-涡相互作用

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

(ρ, u, v, p) = (1.0,1.1 \sqrt{γ}, 0,1.0)

(32)

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

\{\begin{matrix} \tilde{u} = ε_{c} τ e^{α (1 - τ^{2})} s i n θ \\ \tilde{v} = - ε_{c} τ e^{α (1 - τ^{2})} c o s θ \\ \tilde{T} = - \frac{(γ - 1) ε_{c}^{2} e^{2 α (1 - τ^{2})}}{4 α γ} \\ \tilde{S} = 0 \end{matrix}

(33)

其中，τ=r/r_c，

r = \sqrt{{(x - x_{c})}^{2} + {(y - y_{c})}^{2}}

，ε_c=0.3，α=0.204，r_c=0.05，x_c=0.25，y_c=0.5。上下边界取反射边界条件，左右边界取零梯度边界条件。取250×100的均匀网格计算该问题，结束时间t=0.6。同时取1 000×400均匀网格上WENO-JS的计算结果为基准解。因各格式计算结果流场分布接近，图7仅给出WENO-ZF1、WENO-ZF2对应结果压力分布（图中显示从1.19 到 1.37均匀分布的90条压力等值线）。图8给出通过涡核的直线y=0.5上的压力分布。图9为图8压力分布局部放大图。

图7t=0.6时激波-涡相互作用压力分布

Fig.7Pressure distributions of the shock-vortex interaction problem at t=0.6

图8t=0.6时激波-涡相互作用y=0.5压力分布

Fig.8Pressure 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格式。

图9图8局部放大图

Fig.9Locally enlarged plot of Fig.8

表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格式所用相对时间

Tab.2 Relative time of different WENO schemes spent on different cases

4 结论

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

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

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

图15阶WENO-JS/WENO-Z格式重构模板

Fig.1Reconstruction stencils of fifth-order WENO-JS/WENO-Z

下载: 全尺寸图片

图25阶WENO-ZF格式重构模板

Fig.2Reconstruction stencils of fifth-order WENO-ZF

下载: 全尺寸图片

图3各格式ADR分析色散耗散特性

Fig.3ADR dispersion and dissipation of the schemes

下载: 全尺寸图片

图4t=0.4时Sod激波管密度分布

Fig.4Density distributions of the Sod shock tube at t=0.4

下载: 全尺寸图片

图5t=1.8时激波-密度波相互作用密度分布

Fig.5Density distributions of the shock-density wave interaction problem at t=1.8

下载: 全尺寸图片

图6t=0.8时二维黎曼问题密度分布

Fig.6Density distributions of 2D Riemann problem at t=0.8

下载: 全尺寸图片

图7t=0.6时激波-涡相互作用压力分布

Fig.7Pressure distributions of the shock-vortex interaction problem at t=0.6

下载: 全尺寸图片

图8t=0.6时激波-涡相互作用y=0.5压力分布

Fig.8Pressure distributions of the shock-vortex interaction problem on central line y=0.5 at t=0.6

下载: 全尺寸图片

图9图8局部放大图

Fig.9Locally enlarged plot of Fig.8

下载: 全尺寸图片

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

下载: 全尺寸图片

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

下载: 全尺寸图片

图15阶WENO-JS/WENO-Z格式重构模板

Fig.1Reconstruction stencils of fifth-order WENO-JS/WENO-Z

图25阶WENO-ZF格式重构模板

Fig.2Reconstruction stencils of fifth-order WENO-ZF

图3各格式ADR分析色散耗散特性

Fig.3ADR dispersion and dissipation of the schemes

图4t=0.4时Sod激波管密度分布

Fig.4Density distributions of the Sod shock tube at t=0.4

图5t=1.8时激波-密度波相互作用密度分布

Fig.5Density distributions of the shock-density wave interaction problem at t=1.8

图6t=0.8时二维黎曼问题密度分布

Fig.6Density distributions of 2D Riemann problem at t=0.8

图7t=0.6时激波-涡相互作用压力分布

Fig.7Pressure distributions of the shock-vortex interaction problem at t=0.6

图8t=0.6时激波-涡相互作用y=0.5压力分布

Fig.8Pressure distributions of the shock-vortex interaction problem on central line y=0.5 at t=0.6

图9图8局部放大图

Fig.9Locally enlarged plot of Fig.8

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

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

图15阶WENO-JS/WENO-Z格式重构模板

Fig.1Reconstruction stencils of fifth-order WENO-JS/WENO-Z

图25阶WENO-ZF格式重构模板

Fig.2Reconstruction stencils of fifth-order WENO-ZF

图3各格式ADR分析色散耗散特性

Fig.3ADR dispersion and dissipation of the schemes

图4t=0.4时Sod激波管密度分布

Fig.4Density distributions of the Sod shock tube at t=0.4

图5t=1.8时激波-密度波相互作用密度分布

Fig.5Density distributions of the shock-density wave interaction problem at t=1.8

图6t=0.8时二维黎曼问题密度分布

Fig.6Density distributions of 2D Riemann problem at t=0.8

图7t=0.6时激波-涡相互作用压力分布

Fig.7Pressure distributions of the shock-vortex interaction problem at t=0.6

图8t=0.6时激波-涡相互作用y=0.5压力分布

Fig.8Pressure distributions of the shock-vortex interaction problem on central line y=0.5 at t=0.6

图9图8局部放大图

Fig.9Locally enlarged plot of Fig.8

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

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

LIU X D, OSHER S, CHAN T. Weighted essentially non-oscillatory schemes[J]. Journal of Computational Physics,1994,115(1):200-212.

SHU C W, OSHER S. Efficient implementation of essentially non-oscillatory shock-capturing schemes[J]. Journal of Computational Physics,1988,77(2):439-471.

SHU C W, OSHER S. Efficient implementation of essentially non-oscillatory shock-capturing schemes,Ⅱ[J]. Journal of Computational Physics,1989,83(1):32-78.

JIANG G S, SHU C W. Efficient implementation of weighted ENO schemes[J]. Journal of Computational Physics,1996,126(1):202-228.

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.

FENG H, HU F X, WANG R. A new mapped weighted essentially non-oscillatory scheme[J]. Journal of Scientific Computing,2012,51(2):449-473.

FENG H, HUANG C, WANG R. An improved mapped weighted essentially non-oscillatory scheme[J]. Applied Mathematics and Computation,2014,232:453-468.

WANG R, FENG H, HUANG C. A new mapped weighted essentially non-oscillatory method using rational mapping function[J]. Journal of Scientific Computing,2016,67(2):540-580.

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.

LIU S P, SHEN Y Q, ZENG F J,et al. A new weighting method for improving the WENO-Z scheme[J]. International Journal for Numerical Methods in Fluids,2018,87(6):271-291.

CASTRO M, COSTA B, DON W S. High order weighted essentially non-oscillatory WENO-Z schemes for hyperbolic conservation laws[J]. Journal of Computational Physics,2011,230(5):1766-1792.

ACKER F, DE R BORGES R B, COSTA B. An improved WENO-Z scheme[J]. Journal of Computational Physics,2016,313:726-753.

骆信, 吴颂平. 改进的五阶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)

SHEN Y Q, ZHA G C. Improvement of the WENO scheme smoothness estimator[J]. International Journal for Numerical Methods in Fluids,2010,64(6):653-675.

侯中喜, 易仕和, 李桦. 高精度高分辨率WENO格式分析与改进[J]. 国防科技大学学报,2003,25(1):17-20. HOU Z X, YI S H, LI H. Analysis and improvement of high precision,high resolution WENO schemes[J]. Journal of National University of Defense Technology,2003,25(1):17-20.(in Chinese)

柴得林, 李雨润, 孙中国, 等. WENO格式权重分析与改进[J]. 国防科技大学学报,2016,38(4):39-45. CHAI D L, LI Y R, SUN Z G,et al. Analysis and improvement of weights in WENO schemes[J]. Journal of National University of Defense Technology,2016,38(4):39-45.(in Chinese)

BHISE A A, GANDE N R, SAMALA R. An efficient hybrid WENO scheme with a problem independent discontinuity locator[J]. International Journal for Numerical Methods in Fluids,2019,91(1):1-28.

RATHAN S, GANDE N R, BHISE A A. Simple smoothness indicator WENO-Z scheme for hyperbolic conservation laws[J]. Applied Numerical Mathematics,2020,157:255-275.

WANG Y H, DU Y L, ZHAO K L,et al. A low-dissipation third-order weighted essentially nonoscillatory scheme with a new reference smoothness indicator[J]. International Journal for Numerical Methods in Fluids,2020,92(9):1212-1234.

FARDIPOUR K, MANSOUR K. A modified seventh-order WENO scheme with new nonlinear weights for hyperbolic conservation laws[J]. Computers & Mathematics with Applications,2019,78(12):3748-3769.

徐维铮, 孔祥韶, 郑成, 等. 一种三阶WENO-Z格式改进方法[J]. 北京航空航天大学学报,2017,43(12):2400-2405. XU W Z, KONG X S, ZHENG C,et al. An improved method for third-order WENO-Z scheme[J]. Journal of Beijing University of Aeronautics and Astronautics,2017,43(12):2400-2405.(in Chinese)

MARTÍN M P, TAYLOR E M, WU M,et al. A bandwidth-optimized WENO scheme for the effective direct numerical simulation of compressible turbulence[J]. Journal of Computational Physics,2006,220(1):270-289.

HU X Y, WANG Q, ADAMS N A. An adaptive central-upwind weighted essentially non-oscillatory scheme[J]. Journal of Computational Physics,2010,229(23):8952-8965.

ZHU J, QIU J X. A new fifth order finite difference WENO scheme for solving hyperbolic conservation laws[J]. Journal of Computational Physics,2016,318:110-121.

PIROZZOLI S. On the spectral properties of shock-capturing schemes[J]. Journal of Computational Physics,2006,219(2):489-497.

首页

期刊介绍

投稿指南

编委会

出版声明

开放获取声明

联系我们

期刊订阅

Rss

AI检索

English

1 WENO-Z格式

2 重组模板WENO-Z格式

3 数值实验

4 结论