引 言

可压缩流动中存在激波、接触间断等结构. 为了捕捉这些流动结构, 激波捕捉格式应运而生. 早期的激波捕捉格式以总变差减小(total variation diminishing, TVD)格式^[1]、无振荡、无自由参数耗散(non-oscillatory and non-free parameter dissipation, NND)格式^[2]为主. TVD方法是当前工业界最常用的激波捕捉方法. 经典TVD方法在极值点附近降为一阶精度, 整体最多只能达到二阶精度, 不利于获得更细致的流场结构信息. 湍流的大涡模拟、直接数值模拟以及气动声学问题的数值模拟中通常需要捕捉这些细微结构, 这要求格式具有良好的分辨率性质. 如果问题中存在间断, 则可采用具有良好分辨率性质的高精度激波捕捉格式进行模拟.

高精度激波捕捉格式最受关注的一个类别是加权本质无振荡 (weighted essentially non-oscillatory, WENO)格式. WENO格式最初由Liu等^[3]提出, 后由Jiang等^[4]改进并达到了模板上的最优精度. 后续研究者们针对WENO格式做了大量的改进研究, 如为了获得更好流场分辨率的WENO-M^[5], WENO-Z^[6-7], WENO-NS^[8], WENO-CU^[9]等, 为了获得更好收敛性或稳定性的ZSWENO^[10], ESWENO^[11]等. 加权紧致非线性格式(weighted compact nonlinear scheme, WCNS)^[12-13]是另一类借鉴了WENO思想的激波捕捉方法, 张树海^[14]对两类格式进行了比较. 而关于非结构网格WENO格式的构造可参看文献[15-16].

WENO方法采用多个候选模板, 每个候选模板上均有一个逼近值, 对这些逼近值进行加权得到最终的逼近值. 候选模板上函数的光滑程度直接影响这些逼近值的权重. 含间断候选模板的权重几乎为0, 因此可以达到捕捉间断的效果. 在连续区域中, 这些权值需尽量接近于最优权值, 从而减小数值离散误差. 实际计算中, 变量在连续区域小尺度的波动使得WENO格式可能明显偏离线性格式, 给数值模拟带来不利影响.

Wu等^[17]提出了光滑因子单频波保常设计准则, 即光滑因子满足对单频波为常数. 这种光滑因子在一般极值点附近的变化幅度小, 从而使WENO格式更接近线性格式. 按照这一要求, Wu等^[17-18]构造了一类光滑因子, 并设计了相应的WENO-S格式. 由于光滑因子满足单频波保常准则, 因此这类格式的近似色散关系^[19]与线性基底格式一致. 同时, 该类光滑因子在小尺度波动区域变化小, 间断附近变化大, 因而WENO-S格式能更有效地区分间断与小尺度波动, 在数值模拟中能获得良好的结果.

WENO-S格式的光滑因子比经典形式的光滑因子^[4]更简洁, 所需浮点运算少, 因而具有较高的计算效率. 由于其光滑因子仅与子模板上几个函数值相关, 与子模板在大模板中的位置无关. 对于线性对流方程, 计算相邻数值通量时部分光滑因子相同, 无需重复计算. 因此, 可以通过消除冗余的光滑因子计算来提高计算效率.

本文以7阶WENO-S格式为例介绍其构造方法, 包括光滑因子的设计准则, 然后对格式特点进行分析, 接下来讨论计算效率, 针对线性对流方程给出基于消除冗余光滑因子计算的算法, 然后讨论该方法的使用条件, 并将其推广到其他问题的数值模拟中. 最后通过数值实验证明该算法的有效性.

1.   7阶WENO-S格式

1.1   光滑因子单频波保常设计准则

常见WENO格式对于小尺度波动通常表现出过多的耗散. 其原因可归结于极值点附近其权值偏离线性权值较大. 如最常用的Jiang-Shu型r点子模板光滑因子^[4]

式中, $q\left( x \right)$为该子模板上的重构函数, ${x_j}$为大模板的中点, h为空间步长. 在远离临界点的区域, 可由泰勒展开得到$\;\beta _j^{\left( {{\text{JS}}} \right)} = O\left( {{h^2}} \right),$ 临界点表示一阶导数为0的点, 极值点则是一类常见的临界点. 在临界点附近, 则有$\;\beta _j^{\left( {{\text{JS}}} \right)} = O\left( {{h^4}} \right). $显然, 该类光滑因子在极值点附近下降明显, 这使得相应的模板权重变化较大, 对数值模拟产生不利影响.

针对上述情况, 需要减小连续区域光滑因子的变化幅度. 为此, Wu等^[17]提出了光滑因子单频波保常准则. 依此准则设计的光滑因子对于单频波为常数, 相应WENO格式非线性权保持不变, 退化为线性权. 而对于其他一般连续波形, 非线性权的变化也明显减小. 该文认为, 光滑因子应该对单频波为常数的原因如下:

(1)单频波可视为一个单频振荡的函数, 每一部分都具有相同的频率和最大振幅;

(2)如图1所示, 几何意义下圆周处处同样光滑, 而单频波可作为表示圆周横坐标(或者纵坐标)关于角度的函数, 可视为函数意义下处处同样光滑.

图  1  圆周和正弦函数的光滑程度示意图

Figure  1.  Schematic diagram of the smoothness of a circle and a sine function

下载: 全尺寸图片 幻灯片

Wu等^[17-18]给出了一系列满足该准则的光滑因子. 以4点等距模板为例, 对应的光滑因子为

其中${f_j},{f_{j + 1}},{f_{j + 2}},{f_{j + 3}}$为模板上节点函数值. 对于单频波$f\left(x\right) = A\text{ sin}\left(kx + \phi \right) + c$, 其中$A,k,\phi ,c$是给定常量. 通过演算可知, 该光滑因子$\;\beta _j^{\left( {\text{S}} \right)} = 64{A^2}{\sin ^4}\left( {\dfrac{{kh}}{2}} \right){\cos ^2}\left( {\dfrac{{kh}}{2}} \right)$, 为与位置无关的常数.

1.2   7阶WENO-S格式

考虑一维双曲守恒律方程

对于等距网格${x_j} = jh$, 其半离散守恒型差分格式为

其中${\hat f_{j + \frac{1}{2}}}$和${\hat f_{j - \frac{1}{2}}}$为数值通量. ${\hat f_{j \pm \frac{1}{2}}}$可视为函数$s\left( x \right)$在${x_j} \pm \frac{h}{2}$的逼近值, 其中$s\left( x \right)$可由如下隐式给出

下文中公式均假设特征速度${{\partial f}}/{{\partial {{u}}}}$为正. 对于特征速度为负的情况, 可通过对称性获得相应的结果. 若无法确定正负, 则需要采用通量分裂方法, 如Lax-Friedrichs分裂, 时间格式可以采用高阶Runge-Kutta方法, 相关处理可参见文献[4]及其参考文献.

7点WENO格式所用模板为$S = \left\{ {{x_{j - 3}}, {x_{j - 2}}, \cdots , {x_{j + 3}}} \right\}$, 而4个候选子模板分别为${S_0} = \left\{ {{x_{j - 3}},{x_{j - 2}},{x_{j - 1}},{x_j}} \right\},$${S_1} = \left\{ {{x_{j - 2}},{x_{j - 1}},{x_j},{x_{j + 1}}} \right\},$${S_2} = \left\{ {{x_{j - 1}},{x_j},{x_{j + 1}},{x_{j + 2}}} \right\}$, ${{{S}}_3} = \left\{ {{x_j},{x_{j + 1}},{x_{j + 2}},{x_{j + 3}}} \right\}$. 4个子模板上得到的${x_{j + \frac{1}{2}}}$处的重构值为

那么7阶WENO-S格式可写为

${\omega _k}$为非线性权系数, 其公式为

其中

线性权系数${d_k}$为

当${\omega _{j,k}}$取${d_k}$时, WENO格式退化为其线性基底格式. $\varepsilon $用于防止分母为0, 这里取${10^{ - 40}}$. $\;\beta _{j,k}^{\left( {\text{S}} \right)}$为光滑因子, 其公式为

参数$\tau _j^{\left( {\text{S}} \right)}$可视为大模板上的光滑因子, 其计算公式如下

其中$ {c_k} = - {f_{j + k - 3}} + 3{f_{j + k - 2}} - 3{f_{j + k - 2}} + {f_{j + k}}, $ $ {\text{ }}k = 0,1,2,3. $该变量可在计算光滑因子$\;\beta _{j,k}^{\left( {\text{S}} \right)}$的过程中得到.

2.   7阶WENO-S格式的分析

2.1   收敛精度

7阶WENO-S格式在连续区域通常为7阶精度. 但是非线性格式在临界点附近可能发生降阶, 这里临界点指一阶导数为0的点. 临界点的阶数用${n_{{\text{cp}}}}$表示, 如${n_{{\text{cp}}}} = 1$表示$f' = 0,f'' \ne 0,{\text{ }}{n_{{\text{cp}}}} = 2$表示$f' = 0, f'' = 0,f''' \ne 0$. 根据Wu等^[17]的分析, WENO7-S在二阶及以下临界点可保持7阶精度, 而在3阶及以上临界点可能会发生降阶.

对于不小于3阶的临界点, 如果要保持7阶收敛精度, 可在7阶WENO-S格式非线性权的计算中采用与空间步长相关的$\varepsilon $值^[20].

2.2   极值点附近的模拟

极值点附近波形的数值模拟是激波捕捉方法的难点. 正是由于极值点的存在, 经典TVD格式只能达到二阶精度第1.1节中光滑因子单频波保常设计准则的出发点之一也是降低极值点附近光滑因子的变化幅度^[17].

极值点一般位于连续区域, 为获得更好的模拟结果, 通常需要WENO格式的非线性权系数尽量接近于线性权系数. 绝大多数极值点属于一阶临界点. 由Taylor展开可知, 经典光滑因子$\;{\beta ^{\left( {{\text{JS}}} \right)}}$^[4]在一阶临界点附近和一般连续区域分别为$O\left( {{h^4}} \right)$和$O\left( {{h^2}} \right)$. 本文中简略起见, 光滑因子在不涉及到具体点时省去下标. 而7阶WENO-S的光滑因子$\;{\beta ^{\left( {\text{S}} \right)}}$则保持为$O\left( {{h^4}} \right)$, 即光滑因子变化较小. 那么$\;{\beta ^{\left( {\text{S}} \right)}}$不仅对单频波为常数, 而且在一般波形的极值点附近变化也较小. 因此, 7阶WENO-S在极值点附近更接近线性格式, 从而数值模拟误差通常更小.

2.3   间断稳定性

为了消除或减小捕捉间断时的振荡现象, 要求WENO格式中含间断子模板的权系数尽量接近于0. 光滑因子$\;{\beta ^{\left( {\text{S}} \right)}}$在间断区域为$O\left( 1 \right)$, 连续区域为$O\left( {{h^4}} \right)$; 而经典光滑因子$\;{\beta ^{\left( {{\text{JS}}} \right)}}$在两类区域分别为$O\left( 1 \right)$和$O\left( {{h^2}} \right)$, 那么光滑因子$\;{\beta ^{\left( {\text{S}} \right)}}$在两类区域之间的差距更大. 若权系数公式其他部分相同, 则7阶WENO-S格式中含间断子模板的权系数更小. 这使得该格式具有良好的间断稳定性.

2.4   分辨率性质

Lele^[21]阐述了分辨率性质对差分格式的意义, 并采用Fourier方法分析了紧致差分格式的分辨率特性. 但是, Fourier方法不适用于非线性格式. 针对此问题, 文献中最常用的方法是近似色散关系(approximate dispersion relation, ADR)分析^[18]. 该方法首先使用数值方法对不同无量纲波数$\kappa $的单频波进行短时间的推进, 然后采用离散傅里叶变换结合初值计算得到有效波数$\tilde \kappa $. 有效波数的实部代表数值解的相速度, 通常也用来分析色散性能, 而虚部代表耗散率. Li等^[22]和毛枚良等^[23]中将ADR方法进行了简化, 去掉了时间推进过程.

由于WENO-S格式采用了对单频波为常数的光滑因子, 该类格式在计算单频波时退化为线性格式, 因而其近似色散关系与线性基底格式一致, 优于绝大多数同阶WENO格式. 文献[17]中图2给出了几种7点WENO格式及其基底格式的近似色散关系, 其中7阶WENO-S与线性基底格式完全重合.

Zhao等^[24]和Li等^[25-26]中在ADR方法的基础上考虑了非线性响应, 即单频波在非线性格式下产生的其他波数的杂波. Cravero等^[27]中借用热力学中的“温度”来表示非线性格式的这种性质. 由于WENO-S格式对于单频波退化为线性格式而不产生杂波, 因此在这类分析方法中, 该类格式不存在非线性响应, 具有温度0的特征.

3.   7阶WENO-S格式的计算效率

7阶WENO-S格式具有良好的计算效率. 其光滑因子$\;{\beta ^{\left( {\text{S}} \right)}}$比经典光滑因子$\;{\beta ^{\left( {{\text{JS}}} \right)}}$更简洁, 所需计算量更小. 但是该格式的计算效率还需考虑参数$ \tau \text{ } $的计算. Wu等^[17]对比了几种WENO格式单个数值通量的计算量, 并通过数值实验说明了其效率优势. 但是, 对于一些特殊问题, 还可进一步提升该格式的计算效率.

3.1   线性对流方程的效率提升方法

采用经典的WENO格式(3阶WENO除外)时, 各子模板光滑因子的计算公式形式不同, 在相邻数值通量的计算中不会出现相同的光滑因子. 而WENO-S的光滑因子在各子模板上的计算公式除下标不同外形式一致. 那么对于线性对流方程, 计算相邻数值通量时, 部分光滑因子相同. 因而可以通过消除WENO-S的冗余光滑因子计算来提高计算效率.

具体计算中, 可将所有光滑因子先行计算并存储. 另外, 考虑到大模板上的光滑因子${\tau ^{\left( {\text{S}} \right)}}$的计算, 同时需要存储的还有${c_j}' = - {f_j} + 3{f_{j + 1}} - 3{f_{j + 2}} + {f_{j + 3}}$组成的序列. 那么, 当网格点数较多时, 对于7阶WENO-S格式, 子模板光滑因子的计算量约为原来的1/4.

考虑一维线性对流方程

假设对流速度$a$ > 0.设计算网格点为$\left\{ {{x_1},{x_2}, \cdots ,{x_N}} \right\}$. 采用虚拟网格点的方法处理边界条件, 对于7阶WENO格式, 左端需设置4个虚拟点, 右端需设置3个虚拟点. 计算点和虚拟点上的通量组成的序列为$\left\{ {{f_{ - 3}},{f_{ - 2}},{f_{ - 1}}, \cdots ,{f_{N + 3}}} \right\}$. 对于$a$ < 0的情况, 可通过对称性得到.

计算数值通量$\hat f_{j + \frac{1}{2}}^{},{\text{ }}j = 0,1,2, \cdots ,N$的流程如下

(1) 计算${c_j}{{'}}$,$ \text{ }j = -3,-2,\cdots ,N $

(2) 计算子模板光滑因子$\beta _j^{\left( {\text{S}} \right)}$,$ \text{ }j = -3,-2,\cdots ,N $

(3) 采用第1.2节公式计算数值通量$\hat f_{j + \frac{1}{2}}^{}, j = 0,1,2, \cdots ,N$. 其中大模板光滑因子$\tau _j^{\left( {\text{S}} \right)}$中${c_k} = {c_{j + k - 3}}{{'}}$, 而子模板光滑因子取$\beta _{j,k}^{\left( {\text{S}} \right)} = \beta _{j + k - 3}^{\left( {\text{S}} \right)}$.

简便起见, 这里用WENO7-JS表示经典的7阶WENO格式, 而WENO7-Z在WENO7-JS的基础上采用了Z型权系数^[28], 其中幂因子取1, WENO7-S表示7阶WENO-S格式.

表1列出了计算数值通量$\hat f_{j + \frac{1}{2}}^{},{\text{ }}j = 0,1,2, \cdots ,N$时采用几种7阶WENO格式不同计算方法的浮点运算统计. 表中WENO7-JS0中光滑因子${\beta ^{\left( {{\text{JS}}} \right)}}$采用了文献[29]中的计算方法. WENO7-JS1中则采用了${\beta ^{\left( {{\text{JS}}} \right)}}$的快速计算方法^[18,30], 这种方法通过对称矩阵合同变换, 改写了计算公式, 减少了计算量, 本文WENO7-Z也采用这种计算公式. WENO7-S0表示采用标准计算方法的WENO7-S格式, 而WENO7-S1则表示本文提出的快速算法.

表  1  几种WENO格式计算N + 1个数值通量所需浮点运算量

Table  1.  The number of floating point operations required to calculate N + 1 numerical fluxes for several WENO schemes

Scheme	+, −	*, /	|·|	Total
WENO7-JS0	59(N + 1)	87(N + 1)	0	146(N + 1)
WENO7-JS1	55(N + 1)	75(N + 1)	0	130(N + 1)
WENO7-Z	62(N + 1)	75(N + 1)	N + 1	138(N + 1)
WENO7-S0	77(N + 1)	51(N + 1)	5(N + 1)	133(N + 1)
WENO7-S1	47N + 77	39N + 51	2N + 5	88N + 133

下载: 导出CSV 

| 显示表格

表1中可以看到, 采用该快速算法可以有效降低WENO-S格式的计算量, N很大时数值通量计算量降低 5/13$ \approx 38.46\% $. 实际程序运行时, 由于还有其他部分的运算, 其整体计算效率提高的程度应小于该值. 在后文的数值算例测试中, 其计算时间约节省20%.

3.2   其他适用情况

上述效率提升方法也适用于其他一些问题. 其关键在于一条网格线上用于WENO重构或插值的量可以表示为一个序列, 进行网格线上不同位置的WENO重构或插值时, 选取序列中相应连续几个位置的值, 然后进行WENO重构或插值.

除7阶WENO-S格式外, 还有一些其他WENO格式也可采用这种方法提升计算效率. 其要求是光滑因子在各子模板上的计算公式除下标不同外形式一致. 如更高阶的WENO-S格式^[18]、WENO-XS格式^[31]、WENO-LOC格式^[32]和FWENO格式^[30]等.

在稍复杂流动问题的数值模拟中, 一些额外的处理过程可能会破坏这一条件.

(1) 通量分裂

采用守恒型差分方法求解非线性方程时, 考虑到迎风性, 通常采用通量分裂的方法. 如果要采用这种效率提升方法, 则必须采用全局通量分裂. 而局部特征分裂将破坏上述条件, 从而无法采用该方法提升效率.

对于有限体积法和WCNS型^[12]的差分方法, 可采用Godunov类方法计算通量. 该类方法中进行非线性插值或重构的量不是通量, 而是原始变量或守恒变量, 避免了上述全局通量分裂的要求, 可采用该加速方法计算.

(2) 特征投影

对于双曲型偏微分方程组, 由于信号沿特征方向传播的特点, 数值求解采用特征投影处理更符合方程的特点, 在含间断问题中结果更准确. 对于非线性问题, 投影时一般采用局部特征矩阵.

不采用特征投影而直接对各个方程进行数值离散求解的方法称为组元型方法. 组元型方法一般仅用于较弱激波的捕捉, 对于强激波, 该处理可能会产生较强的波后非物理振荡, 对计算结果产生明显不利影响. 但是局部特征投影破坏了该效率提升方法的前提. 那么, 这种效率提升方法只能结合组元型求解方法使用.

某些特定问题中部分方程具有明确的特征方向, 即使剩下的方程组成的方程组采用特征投影求解, 这些方程也可以采用该方法减少计算量. 如多组分问题中的组分方程^[33]、多介质流问题中表示界面的方程^[34]等.

4.   数值算例

本节先进行极值点附近权系数偏离情况测试, 然后采用多个算例对7阶WENO-S的性能进行测试验证, 并且同时测试所提效率提升方法的有效性. 所选算例包括一维对流问题、球面波传播问题、二维旋转问题、二维小扰动传播问题及一维和二维无黏流动问题. 空间离散采用WENO7-JS、WENO7-Z和WENO7-S格式, 时间推进采用3阶TVD Runge-Kutta方法^[4]. 在进行时间效率对比时, WENO7-JS和WENO-7S在格式的后面加0或1以分别代表经典计算方法和快速计算方法.

4.1   极值点附近权系数偏离程度测试

WENO格式的权系数在极值点附近常常明显偏离线性权系数, 这可能导致数值模拟误差的增大. 本算例比较极值点附近几种WENO格式的权系数偏离程度. 下面用$\displaystyle\mathop \sum \nolimits_{k = 0}^3 \left| {{\omega _k} - {d_k}} \right|/4$表示某点通量的权系数偏离值, 统计极值点附近若干通量点的最大偏离值以及平均偏离值. 所选取的3个函数为: (1) ${f_1}\left( x \right) = {{\rm{e}}^{ - 4{x^2}}} $, (2) ${f_2}\left( x \right) = {{\rm{e}}^x} - x - 1 $, (3) ${f_3}\left( x \right) = {\sin ^4}\left( {{\text{π}} x/4} \right)$. 极值点均取$x = 0$. 该点为${f_1}\left( x \right)$和${f_2}\left( x \right)$的一阶临界点, ${f_3}\left( x \right)$的3阶临界点.

计算中取该极值点作为网格点之一, 统计该点左右各4个通量点的情况. 网格间距均取$h = 0.1$,其结果放在表2 中. 可以看到表中WENO7-S的权系数偏离值在极值点相对更小. 其中一阶临界点附近小一个量级甚至更多, 这与第2.2小节的结论相符.

表  2  权系数偏离值的平均值与最大值

Table  2.  The average values and maximum values of the deviations of weights

Function		WENO7-JS	WENO7-Z	WENO7-S
${f_1}\left( x \right)$	average	1.56 × 10−3	2.40 × 10−3	2.58 × 10−4
	minimum	3.61 × 10−2	1.00 × 10−2	7.58 × 10−4
${f_2}\left( x \right)$	average	8.49 × 10−4	9.31 × 10−6	4.38 × 10−8
	minimum	2.40 × 10−3	4.42 × 10−5	7.27 × 10−8
${f_3}\left( x \right)$	average	9.03 × 10−2	4.11 × 10−2	2.07 × 10−2
	minimum	1.87 × 10−1	1.14 × 10−1	5.43 × 10−2

下载: 导出CSV 

| 显示表格

对于${f_3}\left( x \right) = {\sin ^4}\left( {{\text{π}} x/4} \right)$的3阶临界点$x = 0$, 几种格式的权系数偏离值均较大, 这表明WENO格式在这种临界点附近与线性基底格式有明显差异, 最终导致它们可能在此问题中发生降阶. 注意当网格间距进一步减小时, WENO7-JS由于更大的ε值, 其权系数偏离值可能会小于其他两个格式, 使得其更容易在密网格下恢复7阶精度.

4.2   一维对流问题

考虑一组包括高斯波、方波、尖三角波和椭圆波在内的复杂组合包以速度1向右传播的问题^[4], 其控制方程为

初值为

其中的常数$a = 0.5,{\text{ }}z = - 0.7,{\text{ }}\delta = 0.005,$${\text{ }}\alpha = 10,{\text{ }}\beta = \lg 2/\left( {36{\delta ^2}} \right)$. 两端采用周期边界条件.

这里考虑波形的远距传输, 计算终止时间取$T = 2000$, 网格点数取400. 为减小时间离散误差的影响, CFL数取0.01. 计算结果显示在图2中. 对于这种远距传输问题, WENO7-JS耗散过大, 而WENO7-Z则在间断附近出现了明显的数值振荡, 但是在三角波峰值附近表现最好. WENO7-S格式对于其他几种波形均获得了优良的结果. WENO7-JS, WENO7-Z和WENO7-S计算结果的L₁误差依次为2.51 × 10⁻¹, 4.67 × 10⁻²和4.79 × 10⁻². 此问题中WENO7-Z的误差最小, 这可能是因为其三角波附近的优良表现和更小的间断抹平程度.

图  2  组合波的远距传输问题, 网格点数N = 400, 计算终止时间T = 2000

Figure  2.  Long-distance advection of combined waves with the number of grid points N = 400 and the end time T = 2000

下载: 全尺寸图片 幻灯片

表3给出了几种格式的计算时间, 结果表明WENO7-S0比WENO7-JS1略快, 后续算例也大多如此, 而表1中的浮点计算统计则刚好相反. WENO7-JS0和WENO7-Z之间也有类似表现. 在实际计算中, 运行时间与计算平台、浮点计算单元使用率和缓存命中率等均有较大关系^[35] . 本文计算均采用Fortran程序在CPU为AMD^® Ryzen^® R5 4500 U的计算机上运行. WENO7-S0的光滑因子具有相同的形式, 编译器可能将其处理为向量计算, 这会提升其计算速度. WENO7-Z相比WENO7-JS0在光滑因子部分的计算时间减少了, 而τ及权系数的计算时间增大了. 可能在程序运行中, 增大的这一部分时间大于减小的时间.

表  3  组合波远距传输问题的计算时间

Table  3.  Computational time of long-distance advection of combined waves

Scheme	WENO7-JS0	WENO7-JS1	WENO7-Z	WENO7-S0	WENO7-S1
time/s	636.39	580.32	645.67	546.09	416.99

下载: 导出CSV 

| 显示表格

当消除WENO7-S的冗余光滑因子计算后, 其效率得到进一步提高,该问题中WENO7-S1比WENO7-S0提高24.64%.

4.3   球面波传播问题

球面波传播问题是一个计算气动声学的标准测试算例^[36], 其控制方程为

计算区域为[5, 450], 初始条件为$u = 0$, 在$r = 5$处的边界条件为$u = \sin (\omega t)$.

取${\omega } = {{\text{π}} }/{4}$, 计算终止时间为$T = 400$, 空间步长为${\rm{d}}r = 1$. 为减小时间推进误差的影响, 时间步长取${\rm{d}}t = 0.1$. 计算结果显示在图3中. 图中几种格式计算结果的波幅有明显差距, WENO7-JS耗散最大, 而WENO7-S由于其良好的分辨率性质, 获得了最优的计算结果.

图  3  球面波传播问题, 计算终止时间为T = 400, 空间步长为dr = 1. 时间步长为dt = 0.1

Figure  3.  The spherical wave propagation problem with the end time T = 400, the grid length dr = 1, and the time step dt = 0.1

下载: 全尺寸图片 幻灯片

由于该问题计算时间过短, 计算时间差距过小. 因此, 表3给出了时间步长${\rm{d}}t$取0.01时的计算时间统计. 结果表明, WENO7-S计算效率较高, 而消除冗余计算后, 效率进一步提升21.82%.

表  4  球面波传播问题的计算时间, dt = 0.01

Table  4.  Computational time of spherical wave propagation problem with dt = 0.01

Scheme	WENO7-JS0	WENO7-JS1	WENO7-Z	WENO7-S0	WENO7-S1
time/s	1.422	1.336	1.465	1.219	0.953

下载: 导出CSV 

| 显示表格

4.4   二维旋转问题

这个问题来自于Zalesak^[37]的文章. 该问题中, 一个割开的圆柱绕着附近的轴旋转. 该运动的控制方程为

这里$u = - {{\varOmega }}\left( {y - {y_0}} \right),v = {{\varOmega }}\left( {x - {x_0}} \right)$, 其中 ${{\varOmega }}\; (=2{\text{π}} /360)$ 是旋转角速度. 这里计算区域为$\left( {x,y} \right) \in \left[ 0,10\left] \times \right[0,10 \right], ({x_0},{y_0} ) = \left( {5,5} \right)$是旋转轴的坐标. 初始时刻, 割开圆柱的中心坐标为$\left. {({x_{\text{c}}},{y_{\text{c}}}} \right) = \left( {5,7.5} \right)$, 圆柱半径为1.5, 两条割线的横坐标分别为4.75和5.25, 割开截止位置的纵坐标为8.5; 割开圆柱上的$\phi $值为3.0, 其他区域则是1.0. 根据这个控制方程, 旋转过程中割开圆柱应该保持原始的形状. 本问题将该割开圆柱旋转5个周期.

该问题中两个分量速度既可能为正又可能为负, 通常同时包含正负速度时需要进行通量分裂. 但是该问题每条网格线上特征速度可确定为正或负, 无需通量分裂处理.

采用200 × 200均匀网格进行计算, 网格点坐标为${x_i} = \left( {i - 0.5} \right) \cdot {{\Delta }}x$, ${y_j} = \left( {j - 0.5} \right) \cdot {{\Delta }}y$, 其中$\text{Δ}x和\text{Δ}y$分别为两个方向的网格间距. 图4给出了计算结果, 其中高度代表$\phi $值. 图中可以看到, WENO7-JS对圆柱边缘处的抹平最为严重, 顶部区域有塌陷现象. WENO7-Z的边缘更锐利, 但是附近有轻微的过冲现象. WENO7-S则避免了过度的抹平和边缘处过冲现象.

图  4  二维旋转问题计算结果与精确解

Figure  4.  The computational results and the exact solution of the two-dimensional rotation problem

下载: 全尺寸图片 幻灯片

为了更清晰地显示计算结果的差异, 截取直线$y = {y_{150}} = 7.475$和$x = {x_{100}} = 4.975$上的结果进行对比. 图5(a)为两条直线的位置示意图, 图5(b)和图5(c)给出了两条线上的结果. 从图5(b)中可以看到WENO7-JS在割开部分表现出明显的耗散, WENO7-Z和WENO7-S则出现了下冲, 其中WENO7-Z的下冲程度相对较大, 这也与图5(c)中区间[6, 8.5]的结果相符. 图5(b)中WENO7-Z在$x = 6.5$附近出现了微弱的上冲. 图5(c)中区间[8.5, 9]处WENO7-JS的峰值明显小于精确值3, 这也对应图4(b)中的顶部塌陷现象.

图  5  二维旋转问题两条线上的结果对比图

Figure  5.  Comparison of results on two lines of the two-dimensional rotation problem

下载: 全尺寸图片 幻灯片

采用

来表示上冲/下冲幅度. WENO7-JS, WENO7-Z和WENO7-S的$\delta $值分别为1.158 × 10⁻⁴, 0.156和5.944 × 10⁻². 这说明WENO7-Z更容易出现上冲/下冲现象. 而3种格式的L₁误差依次为2.298, 1.774和1.959. 说明尽管WENO7-Z有较明显的上冲/下冲, 但是其误差相对最小. 从图5(b)和图5(c)中也可看出, WENO7-Z的间断抹平程度最小, 而间断附近贡献了大部分误差量, 这使得其误差小于WENO7-S.

为了考察不同网格下的计算情况, 采用200 × 200, 400 × 400和800 × 800这3套均匀网格对该问题进行模拟, 其时间步长依次设为0.1, 0.05和0.025. 为减少模拟时间, 这里仅将圆柱旋转一个周期. 计算结果的L₁误差和上冲/下冲幅度由表5给出. 显然WENO7-Z的L₁误差依然最小. 而由于间断的存在, 3种格式的L₁误差收敛精度均约为1阶. 对于上冲/下冲幅度$\delta $, 200 × 200和400 × 400网格时WENO7-JS最小, 而800 × 800时WENO7-S最小. 随着网格的加密, WENO7-JS和WENO7-Z没有表现出明显的变化趋势, 而WENO7-S的$\delta $值迅速减小, 体现了其良好的间断稳定性.

表  5  几种格式不同网格下计算结果的L₁误差和上冲/下冲幅度

Table  5.  The L₁ error and up/down overshooting amplitude of the computational results with different schemes and grids

	Grid	WENO7-JS	WENO7-Z	WENO7-S
L1 error	200 × 200	1.759	1.428	1.600
	400 × 400	0.959	0.788	0.886
	800 × 800	0.524	0.430	0.484
up/down overshooting amplitude $\delta$	200 × 200	1.11 × 10−4	1.45 × 10−2	5.02 × 10−3
	400 × 400	9.86 × 10−5	9.13 × 10−3	2.72 × 10−4
	800 × 800	1.27 × 10−4	2.58 × 10−2	2.87 × 10−8

下载: 导出CSV 

| 显示表格

表6给出了200 × 200网格旋转5个周期的计算时间对比. 结果表明WENO7-S计算效率较高, 新的计算方法进一步提升了其计算效率, 该问题中提升23.48%.

表  6  二维旋转问题的计算时间

Table  6.  Computational time for the two-dimensional rotation problem

Scheme	WENO7-JS0	WENO7-JS1	WENO7-Z	WENO7-S0	WENO7-S1
time/s	114.731	107.400	117.493	99.744	76.327

下载: 导出CSV 

| 显示表格

4.5   二维小扰动传播问题

控制方程采用二维线化Euler方程

其中

${M_x} = 0.5,{\text{ }}{M_y} = 0$是$x$和$y$方向的平均流马赫数. ${\boldsymbol{U}}{{'}}$表示密度、速度及压强的扰动量. 计算区域为$\left( {x,y} \right) \in \left[ { - 100,100} \right] \times \left[ { - 100,100} \right]$. 初值为

该问题为计算气动声学的一个标准算例^[36,38].

根据迎风格式的要求, 采用了每条网格线上的全局通量分裂. x方向计算时正负通量分别为${{\boldsymbol{F}}'^ \pm } = \dfrac{1}{2}\left[ {{\boldsymbol{F}}' \pm \left( {{M_x} + 1} \right){\boldsymbol{U}}{{'}}} \right].$ y方向计算时正负通量分别为${{\boldsymbol{G}}'^ \pm } = \dfrac{1}{2}\left[ {{\boldsymbol{G}}' \pm \left( {{M_y} + 1} \right){\boldsymbol{U}}{{'}}} \right]$.

计算采用的网格为200 × 200, 终止时间为$T = 30$. 图6给出了水平中心线$y = 0$上的密度分布图. 可以看到几种WENO格式均取得了与精确解非常接近的结果. 相比而言, WENO7-S在峰值附近取得了最佳的结果. 表7给出了几种格式的计算时间. 本算例中WENO7-S计算效率最高, 且还可受益于本文计算效率提升方法, 该问题中提升17.84%.

图  6  二维小扰动传播问题计算结果, 网格200 × 200, 终止时间T = 30

Figure  6.  Computational results of two-dimensional small disturbance propagation problem, with grid 200 × 200 and end time T = 30

下载: 全尺寸图片 幻灯片

表  7  二维小扰动传播问题的计算时间

Table  7.  Computational time of two-dimensional small disturbance propagation problem

Scheme	WENO7-JS0	WENO7-JS1	WENO7-Z	WENO7-S0	WENO7-S1
time/s	8.842	7.780	8.377	7.485	6.150

下载: 导出CSV 

| 显示表格

4.6   一维无黏流动问题

控制方程为一维Euler方程

其中

式中$\rho ,u,p,E$表示密度、速度、压强和内能. 气体状态方程为

其中比热比$\gamma = 1.4$.

这里计算Shu-Osher问题, 该问题中马赫数3的激波与正弦波相互干扰, 诱发高频振动, 其初值条件为

计算终止时间$T = 1.8$, 这里取网格点为200.

分别采用组元型和特征型方法结合几种几种WENO格式进行计算. 采用5阶WENO-JS在4000个均匀网格点的计算结果作为参考解, 图7(a)为参考解的密度分布图. 图7(b)对比了几种WENO格式的计算结果. 图例中几种WENO格式省略了名称前缀WENO7, 并用后缀Comp表示组元型方法, Char表示特征型方法.

图  7  Shu-Osher 问题计算结果的密度对比

Figure  7.  The density distributions of computational results of Shu-Osher problem

下载: 全尺寸图片 幻灯片

从图7(b)中可以看出, 在高频震荡区域, 特征型的结果更好, 这在WENO7-JS格式的结果中体现得更为明显. 在图中高频震荡左侧, WENO7-Z和WENO7-S的组元型结果相对更接近参考解. 几种WENO格式中, WENO7-S的结果最好, 其中特征型的结果比组元型的结果略好.

该问题计算时间较短, 为此将网格数设为2000, 此时特征型WENO7-JS的计算时间为5.317 s. 对于组元型计算方法, 表8给出了几种方法的计算时间. WENO7-JS特征型的计算时间约为组元型的1.6倍. 组元型方法中WENO7-S格式计算时间最短, 其中本文加速方法提升了其25.26%的计算速度.

表  8  Shu-Osher问题计算时间, N = 2000

Table  8.  Computational time of Shu-Osher problem, N = 2000

Scheme	WENO7-JS0	WENO7-JS1	WENO7-Z	WENO7-S0	WENO7-S1
time/s	3.330	3.037	3.317	3.015	2.255

下载: 导出CSV 

| 显示表格

该问题中激波马赫数为3, 并非弱激波. 此时组元型和特征型的结果之间的差距可能相对较明显, 如WENO7-JS格式. 但是用WENO7-S格式计算时两种方法的差距较小, 可考虑采用组元型方法并结合本文加速方法.

4.7   二维无黏流动问题

控制方程为守恒形式的二维Euler方程

其中

式中$\rho ,u,v,p,E$表示密度、$x$向速度、$y$向速度、压强和内能. 气体状态方程为

其中比热比$\gamma = 1.4$.

这里采用二维激波旋涡相互作用问题作为测试算例^[4] , 该问题描述了一个静止的激波和旋涡的相互作用. 计算区域为$\left[ {0,2} \right] \times \left[ {0,1} \right]$. 一个马赫数为1.1的静止激波放置在$x$ = 0.5处且与$x$轴垂直. 左状态为$\left( {\rho ,u,v,p} \right) = \left( {1,{\text{ }}1.1\sqrt \gamma ,{\text{ }}0,{\text{ }}1} \right)$, 右状态由Rankine-Hugoniot关系给出. 一个中心在$\left( {{x_c},{y_c}} \right) = \left( {0.25,0.5} \right)$的小旋涡放在激波左侧. 这个旋涡可被视为速度$\left( {u,v} \right)$、熵$S = \log \dfrac{p}{{{\rho ^\gamma }}}$, 温度$T = p/\rho $的小扰动, 即

其中 $\tau = r/R$, $r = \sqrt {{{\left( {x - {x_c}} \right)}^2} + {{\left( {y - {y_c}} \right)}^2}} $. 这里$\varepsilon $表示旋涡的强度, $\alpha $控制旋涡的衰减速度, $R$是旋涡最大强度处的半径. 该问题中取$\varepsilon = 0.3,{\text{ }}R = 0.05,{\text{ }}\alpha = 0.204.$ 计算终止时间取$T = 0.6,$采用均匀网格$200 \times 100.$

由于特征投影方法采用了流场当地的投影矩阵, WENO7-S1的效率提升方法无法使用. 这里采用组元型方法, 该方法计算量比特征投影方法小, 但是不宜用于强激波问题. 该问题来流马赫数为1.1, 说明激波较弱, 计算结果对比也表明两种方法差距很小. 另外, 特征型WENO格式中特征投影及逆投影所需计算量很大. 本算例采用WENO7-JS1结合特征型方法时计算时间为41.893 s, 约为该格式结合组元型方法计算时间13.195 s的3倍. 这高于一维问题中的1.6倍. 原因是特征投影矩阵在二维问题中是四维, 其乘法计算量显著高于一维问题中的三维矩阵. 可以预计, 在三维问题中, 特征型算法的计算时间相比组元型多3倍以上.

采用组元型方法计算时, $x$方向正负通量分别为${F^ \pm } = \dfrac{1}{2}\left( {F \pm {\lambda _x}U} \right),$ 其中${\lambda _x}$为该$x$向网格线上所有矩阵$\partial F/\partial U$特征值绝对值的最大值. $y$方向计算时正负通量分别为${G^ \pm } = \dfrac{1}{2}\left( {G \pm {\lambda _y}U} \right)$, 其中${\lambda _y}$为该$y$向网格线上所有矩阵$\partial G/\partial U$特征值绝对值的最大值.

图8中给出了计算结果的压强云图, 其等值线范围为1.19 ~ 1.37, 共90条. 其中图8(a)为特征型WENO7-JS的结果, 可以看到在较弱的激波下, 组元型和特征型结果相差较小. 组元型方法的计算结果中, WENO7-JS和WENO7-S较为相近, 而WENO7-Z在激波附近聚集了更多等值线. 这可能是由于WENO7-Z在间断附近耗散更小, 激波波后非物理振荡更为明显, 而这种振荡会导致相应值的等值线多次出现.

图  8  激波旋涡相互作用问题, 压强等值线范围1.19-1.37, 共90条, 网格200 × 100, 终止时间T = 0.6

Figure  8.  Shock vortex interaction problem. 90 pressure isolines ranging from 1.19 to 1.37. Component-wise seventh-order WENO schemes with grid 200 × 100 and end time T = 0.6

下载: 全尺寸图片 幻灯片

表9给出了几种格式的计算时间. 结果表明WENO7-S 格式效率相对较高, WENO7-S1的效率提升方法进一步减少了21.66%的计算时间.

表  9  激波旋涡相互作用问题的计算时间

Table  9.  Computational time of shock vortex interaction problem

Scheme	WENO7-JS0	WENO7-JS1	WENO7-Z	WENO7-S0	WENO7-S1
time/s	14.350	13.195	14.520	13.007	10.159

下载: 导出CSV 

| 显示表格

5.   结论

WENO-S格式具有许多优良的性质. 特别是其光滑因子对单频波为常数, 使得其是一类满足近似色散关系与线性基底格式一致的激波捕捉格式. 另外, 该格式能有效区分间断与小尺度波动, 在连续区域误差小, 间断附近也能保证良好的稳定性, 在数值模拟中能获得良好的结果.

由于WENO-S格式的光滑因子在各子模板上的计算公式除下标不同外形式一致, 在计算线性对流方程的相邻数值通量时, 部分光滑因子相同. 本文基于7阶WENO-S格式, 介绍了如何通过避免冗余光滑因子计算提高格式计算效率. 除WENO-S格式外, 还有一些其他WENO格式也可采用这种方法提升计算效率.

这种方法可推广至多种问题, 其应用条件是一条网格线上用于WENO重构或插值的量可以表示为一个序列. 因此在非线性方程中使用通量分裂时只能采用全局通量分裂. 非线性方程的另一类处理方法是对原始变量或守恒变量直接进行WENO重构或插值, 然后采用Godunov类方法求得数值通量. 这时也可以采用避免冗余光滑因子计算的算法.

对于可压缩流动问题, 局部特征投影求解会破坏这种消除冗余计算的方法. 因此, 在可压缩流动模拟中, 该方法只能结合组元型方法求解. 某些流动问题中部分方程的求解不受该限制, 仍可采用该方法减小计算量, 如多组分问题中的组分方程、多介质流问题中表示界面的方程等.

数值实验结果表明7阶WENO-S格式同时具有良好的小尺度结构分辨率和激波捕捉能力. 在计算效率方面, 本文所提方法能有效减少7阶WENO-S格式约20%的计算时间.