一种提高极值点处精度的三阶WENO-Z改进格式及应用

徐维铮; 吴卫国

doi:10.11858/gywlxb.20170696

一种提高极值点处精度的三阶WENO-Z改进格式及应用

doi: 10.11858/gywlxb.20170696

徐维铮^{1, 2,},
吴卫国^{1, 2,*, ,}

1.
武汉理工大学高性能舰船技术教育部重点实验室, 湖北武汉 430063
2.
武汉理工大学交通学院船舶、海洋与结构工程系, 湖北武汉 430063

基金项目:

装备预研教育部联合基金（青年人才） 6141A020331

国家自然科学基金 51409202

中央高校基本科研业务费 2016-YB-016

详细信息

作者简介:
徐维铮(1991-), 男, 博士, 主要从事爆炸波数值计算方法及程序开发研究.E-mail:xuweizheng@whut.edu.cn

通讯作者:
吴卫国(1960-), 男, 教授, 博士生导师, 主要从事结构动力学及舰船抗爆抗冲击研究.E-mail:mailjt@163.com

中图分类号: O357.1
计量
- 文章访问数: 7680
- HTML全文浏览量: 3698
- PDF下载量: 214
出版历程
- 收稿日期: 2017-12-29
- 修回日期: 2018-02-11

An Improved Third-Order WENO-Z Scheme for Achieving Optimal Order near Critical Points and Its Application

XU Weizheng^{1, 2
,},
WU Weiguo^{1, 2,*
, ,}

1.
Key Laboratory of High Performance Ship Technology of Ministry of Education, Wuhan University of Technology, Wuhan 430063, China
2.
Departments of Naval Architecture, Ocean and Structural Engineering, School of Transportation, Wuhan University of Technology, Wuhan 430063, China

摘要

摘要: 高精度、高分辨率激波捕捉格式对含激波流场的数值模拟具有重要意义。为了提高三阶WENO-Z格式在极值点处的计算精度，首先通过理论推导给出三阶WENO格式满足收敛精度的充分条件。采用泰勒级数展开的方式，推导确定所构造格式的参数。通过精度测试证明改进格式在光滑流场区域能收敛到三阶精度。选用Sod激波管、Rayleigh-Taylor不稳定性等经典算例证实了提出的改进格式WENO-NN3相较其他格式（WENO-SJ3、WENO-Z3和WENO-N3）具有精度高、耗散低、对流场结构分辨率高的特性。
- 三阶WENO格式 /
- 光滑因子 /
- 高精度 /
- 高分辨率 /
- 双曲守恒律
Abstract: A high-precision and resolution shock capturing scheme is of great significance for numerical simulation of the complex flow field containing shock waves.In this study, to improve the convergence accuracy of the conventional third-order WENO-Z scheme at the critical points, we firstly derived the sufficient conditions for satisfying the convergence precision of the third-order WENO scheme from the theoretical derivation, then determined the parameters of the constructed scheme using the Taylor series expansion for satisfying the sufficient conditions, and proved using the accuracy test that the proposed scheme converges to the third order precision in smooth flow field including the critical points.Furthermore, we selected the Sod shock tube, the Rayleigh-Taylor instability and some other classic examples, verifying that the improved scheme WENO-NN3 was capable of giving more precision and high resolution results of the complex flow field structures compared with other WENO schemes such as the WENO-JS3, WENO-Z3, and WENO-N3.
- third-order WENO scheme /
- smoothness indicators /
- high precision /
- high resolution /
- hyperbolic conservation law

HTML全文

冲击波流场包含激波间断和光滑流场等复杂结构，对冲击波流场的模拟需要高精度、低耗散的激波捕捉格式。Liu等^[1]于1994年，在ENO格式^[2]构造思想的基础上首次提出WENO格式(Weighted Essentially Non-Oscillatory Scheme)的构造方法。Jiang等^[3]在改进原始WENO格式光滑因子的基础上，提出了三阶、五阶WENO格式的构造框架。目前，WENO格式作为一种典型的高精度激波捕捉格式，对流场内的激波间断具有较高的分辨率，适于求解包含激波、膨胀波以及接触间断等复杂结构的流场，并被大量研究者采用^[4-7]。然而WENO格式在极值点处会降阶，Henrick等^[8]首先指出五阶WENO格式在连续解极值点处精度会降低(五阶精度WENO格式在一阶极值点处会降为三阶)的缺陷，并提出基于映射函数的五阶WENO-M格式，之后，Borges等^[9]采用线性组合低阶模板光滑因子构造高阶全局光滑因子的方法，提出五阶WENO-Z格式。在文献[8-9]改进思想的基础上，出现了大量的改进高阶WENO格式(主要集中在五阶、七阶和九阶)在极值点处精度的文献^[10-14]。针对三阶WENO格式的改进，Yamaleev等^[15]在改进三阶WENO-Z格式中全局光滑因子的基础上，通过理论推导给出三阶能量稳定WENO格式(ESWENO)。Wu等^[16-17]首先根据理论推导指出传统的三阶WENO-Z格式^[18]在极值点处会降阶，其根据Hu等^[19]提出的构造思想，通过线性组合局部模板和全局模板光滑因子的方式，推导给出改进的三阶WENO格式(WENO-N3、WENO-NP3)。

在上述研究的基础上，提出一种改进的三阶WENO-Z格式，提高格式在极值点处的计算精度。选取Sod激波管、激波与熵波相互作用、Rayleigh-Taylor等经典算例，考察了改进格式WENO-NN3的计算性能。研究表明，改进格式不仅提高了极值点处的计算精度，还降低了格式的耗散，提高了对流场结构的分辨率。

1. 控制方程

含激波流场采用可压缩欧拉方程进行描述，其具体形式为

${\mathit{\boldsymbol{Q}}_t} + {\mathit{\boldsymbol{E}}_x} + {\mathit{\boldsymbol{F}}_y} = 0$

(1)

$\mathit{\boldsymbol{Q}} = \left( {\begin{array}{*{20}{c}} \rho \\ {\rho u}\\ {\rho v}\\ E \end{array}} \right),\mathit{\boldsymbol{E}} = \left( {\begin{array}{*{20}{c}} {\rho u}\\ {\rho {u^2} + p}\\ {\rho uv}\\ {u\left( {E + p} \right)} \end{array}} \right),\mathit{\boldsymbol{F}} = \left( {\begin{array}{*{20}{c}} {\rho v}\\ {\rho vu}\\ {\rho {v^2} + p}\\ {v\left( {E + p} \right)} \end{array}} \right)$

(2)

$E = \rho e + \frac{1}{2}\rho {u^2} + \frac{1}{2}\rho {v^2}$

(3)

$p = \left( {\gamma - 1} \right)\rho e$

(4)

式中:ρ是密度; u、v是x、y方向上的速度分量; p为流体压力; E是单位体积流体的总能量; e是比内能; γ表示气体的绝热指数，本文中取为1.4。

在方程(1)的每个方向上均可以看成是一个双曲守恒律方程

$\frac{{\partial u}}{{\partial t}} + \frac{{\partial f\left( u \right)}}{{\partial x}} = 0$

(5)

例如，针对x方向，方程(5)的数值离散形式为

$\frac{{{\rm{d}}{u_i}}}{{{\rm{d}}t}} = - \frac{{\partial f}}{{\partial x}}\left| {_{x = {x_i}}} \right. = - \frac{{{h_{i + 1/2}} - {h_{i - 1/2}}}}{{\Delta x}} \approx = - \frac{1}{{\Delta x}}\left( {{{\bar f}_{i + 1/2}} - {{\bar f}_{i - 1/2}}} \right)$

(6)

式中:f_i－1/2、f_i+1/2分别为单元(x_i-1/2, x_i+1/2)的左、右对流项数值通量，Δx为x方向的均匀网格间距，在本数值计算中，x、y两个方向的网格间距相同。控制方程的空间离散采用下述的数值方法，时间项离散采用三阶TVD龙格库塔格式^[20]。

2. 数值方法

2.1 三阶WENO-JS格式

2.1.1 数值重构过程

三阶WENO格式(WENO-JS3)的数值离散和推导过程如下^[3](为了简洁，仅仅给出右界面通量f_i+1/2的重构过程，左界面的数值通量重构与此相同，只是向左更换了一组模板)。对于三阶WENO格式，f_i+1/2的2种重构方式分别为

$f_{i + 1/2}^0 = - \frac{1}{2}{f_{i - 1}} + \frac{3}{2}{f_i},\;\;\;f_{i + 1/2}^1 = \frac{1}{2}{f_i} + \frac{1}{2}{f_{i + 1}}$

(7)

利用上述两种二阶通量的凸组合计算最终具有三阶精度的数值通量f_i+1/2，即

${{\bar f}_{i + 1/2}} = \sum\limits_{k = 0}^1 {{\omega _k}f_{i + 1/2}^k} = {\omega _0}f_{i + 1/2}^0 + {\omega _1}f_{i + 1/2}^1$

(8)

对于光滑情形，有 ${\omega _0} = {d_0} = \frac{1}{3}$ ， ${\omega _1} = {d_1} = \frac{2}{3}$ 。(8)式所给出的形式既适合光滑流场也适合含间断流场，对于含激波间断流场，式中的非线性权函数ω_k需要根据下式求得

${\omega _k} = \frac{{{\alpha _k}}}{{\sum\limits_{s = 0}^1 {{\alpha _s}} }},\;\;\;\;{\alpha _k} = \frac{{{d_k}}}{{{{\left( {\varepsilon + {\beta _k}} \right)}^2}}},\;\;\;k = 0,1$

(9)

式中:参数ε取值为10^-6。光滑因子β_k(k=0, 1)的表达式为

${\beta _0} = {\left( {{f_{i - 1}} - {f_i}} \right)^2},\;\;\;\;{\beta _1} = {\left( {{f_i} - {f_{i + 1}}} \right)^2},$

(10)

2.1.2 收敛精度的充分条件

在文献[15, 21]研究的基础上，对三阶WENO格式的收敛精度充分条件进行详细的推导证明。对f_i±1/2^k进行泰勒级数展开可得

$f_{i \pm 1/2}^k = {h_{i \pm 1/2}} + {A_k}\Delta {x^2} + o\left( {\Delta {x^3}} \right),\;\;\;\;k = 0,1$

(11)

根据(8)式、(11)式可得

${{\bar f}_{i \pm 1/2}} = \sum\limits_{k = 0}^1 {{\omega _k}f_{i + 1/2}^k} = \sum\limits_{k = 0}^1 {{d_k}f_{i \pm 1/2}^k} + \sum\limits_{k = 0}^1 {\left( {{\omega _k} - {d_k}} \right)f_{i \pm 1/2}^k}$

(12)

紧接着根据(12)式，可得单元(x_i-1/2, x_i+1/2)右界面通量f_i+1/2

$\begin{array}{l} {{\bar f}_{i + 1/2}} = \sum\limits_{k = 0}^1 {{d_k}f_{i + 1/2}^k} + \sum\limits_{k = 0}^1 {\left( {\omega _k^ + - {d_k}} \right)f_{i + 1/2}^k} \\ \;\;\;\;\;\;\;\; = \sum\limits_{k = 0}^1 {{d_k}f_{i + 1/2}^k} + \sum\limits_{k = 0}^1 {\left( {\omega _k^ + - {d_k}} \right)\left( {{h_{i + 1/2}} + {A_k}\Delta {x^2} + o\left( {\Delta {x^3}} \right)} \right)} \\ \;\;\;\;\;\;\;\; = {h_{i + 1/2}} + {B^ + }\Delta {x^3} + o\left( {\Delta {x^4}} \right) + {h_{i + 1/2}}\sum\limits_{k = 0}^1 {\left( {\omega _k^ + - {d_k}} \right)} + \\ \;\;\;\;\;\;\;\;\;\;\;\Delta {x^2}\sum\limits_{k = 0}^1 {{A_k}\left( {\omega _k^ + - {d_k}} \right)} + \sum\limits_{k = 0}^{r - 1} {\left( {\omega _k^ + - {d_k}} \right)o\left( {\Delta {x^3}} \right)} \end{array}$

(13)

同理，左界面通量f_i－1/2为

$\begin{array}{l} {{\bar f}_{i - 1/2}} = \sum\limits_{k = 0}^1 {{d_k}f_{i - 1/2}^k} + \sum\limits_{k = 0}^1 {\left( {\omega _k^ - - {d_k}} \right)f_{i - 1/2}^k} \\ \;\;\;\;\;\;\;\; = \sum\limits_{k = 0}^1 {{d_k}f_{i - 1/2}^k} + \sum\limits_{k = 0}^1 {\left( {\omega _k^ - - {d_k}} \right)\left( {{h_{i - 1/2}} + {A_k}\Delta {x^2} + o\left( {\Delta {x^3}} \right)} \right)} \\ \;\;\;\;\;\;\;\; = {h_{i - 1/2}} + {B^ - }\Delta {x^3} + o\left( {\Delta {x^4}} \right) + {h_{i - 1/2}}\sum\limits_{k = 0}^1 {\left( {\omega _k^ - - {d_k}} \right)} + \\ \;\;\;\;\;\;\;\;\;\;\;\Delta {x^2}\sum\limits_{k = 0}^1 {{A_k}\left( {\omega _k^ - - {d_k}} \right)} + \sum\limits_{k = 0}^{r - 1} {\left( {\omega _k^ - - {d_k}} \right)o\left( {\Delta {x^3}} \right)} \end{array}$

(14)

式中:B⁺=B^－，ω_k^－、ω_k⁺分别对应左界面通量和右界面通量。进一步根据(6)式可得

$\begin{array}{*{35}{l}} \frac{{{{\bar{f}}}_{i+1/2}}-{{{\bar{f}}}_{i-1/2}}}{\Delta x}= \\ \frac{{{h}_{i+1/2}}-{{h}_{i-1/2}}}{\Delta x}+\frac{1}{\Delta x}\left[ {{h}_{i+1/2}}\sum\limits_{k=0}^{1}{\left( \omega _{k}^{+}-{{d}_{k}} \right)}-{{h}_{i-1/2}}\sum\limits_{k=0}^{1}{\left( \omega _{k}^{-}-{{d}_{k}} \right)} \right]+ \\ \Delta x\sum\limits_{k=0}^{1}{{{A}_{k}}\left( \omega _{k}^{+}-\omega _{k}^{-} \right)}+\sum\limits_{k=0}^{r-1}{\left( \omega _{k}^{+}-{{d}_{k}} \right)o\left( \Delta {{x}^{2}} \right)}-\sum\limits_{k=0}^{r-1}{\left( \omega _{k}^{-}-{{d}_{k}} \right)o\left( \Delta {{x}^{2}} \right)} \\ ={{{{f}'}}_{i}}+o\left( \Delta {{x}^{3}} \right)+\frac{1}{\Delta x}\left[ {{h}_{i+1/2}}\sum\limits_{k=0}^{1}{\left( \omega _{k}^{+}-{{d}_{k}} \right)}-{{h}_{i-1/2}}\sum\limits_{k=0}^{1}{\left( \omega _{k}^{-}-{{d}_{k}} \right)} \right]+ \\ \Delta x\sum\limits_{k=0}^{1}{{{A}_{k}}\left( \omega _{k}^{+}-\omega _{k}^{-} \right)}+\sum\limits_{k=0}^{r-1}{\left( \omega _{k}^{+}-{{d}_{k}} \right)o\left( \Delta {{x}^{2}} \right)}-\sum\limits_{k=0}^{r-1}{\left( \omega _{k}^{-}-{{d}_{k}} \right)o\left( \Delta {{x}^{2}} \right)} \\ \end{array}$

(15)

根据(15)式可知，三阶WENO格式满足三阶收敛精度的充分必要条件

$\left\{ \begin{array}{l} \sum\limits_{k = 0}^1 {\left( {\omega _k^ \pm - {d_k}} \right) = o\left( {\Delta {x^4}} \right)} \\ \sum\limits_{k = 0}^1 {{A_k}\left( {\omega _k^ + - \omega _k^ - } \right) = o\left( {\Delta {x^2}} \right)} \\ \left( {\omega _k^ \pm - {d_k}} \right) = o\left( {\Delta x} \right) \end{array} \right.$

(16)

然而很难根据(16)式进行格式非线性权重的设计，因此需要寻求更简洁的限制条件。根据(16)式，可以给出如下简洁的充分条件

$\omega _k^ \pm - {d_k} = o\left( {\Delta {x^2}} \right)$

(17)

2.2 三阶WENO-NP3格式

首先给出三阶WENO-Z格式^[18](WENO-Z3)的具体形式

$\alpha = {d_k}\left( {1 + \frac{\tau }{{{\beta _k} + \varepsilon }}} \right)$

(18)

$\tau = \left| {{\beta _0} - {\beta _1}} \right|$

(19)

Wu等^[16]通过研究发现WENO-Z3格式在一阶极值点处精度降为一阶，并提出可改善分辨率的WENO-N3格式。为了提高三阶WENO-Z格式在极值点处精度，其在文献[17]中紧接着提出WENO-NP3格式，具体形式为

${\alpha _k} = {d_k}\left( {1 + \frac{{{\tau _{{\rm{NP}}}}}}{{{\beta _k} + \varepsilon }}} \right),\;\;\;\;k = 0,1$

(20)

${\tau _{{\rm{NP}}}} = {\left| {\frac{{{\beta _0} + {\beta _1}}}{2} - {\beta _3}} \right|^p},\;\;\;\;\varepsilon = {10^{ - 40}},\;\;\;p = 3/2$

(21)

式中:β₃表示三阶WENO格式全局模板(x_i-1, x_i, x_i+1)的光滑因子

${\beta _3} = \frac{{13}}{{12}}{\left( {{f_{i - 1}} - 2{f_i} + {f_{i + 1}}} \right)^2} + \frac{1}{4}{\left( {{f_{i - 1}} - {f_{i + 1}}} \right)^2}$

(22)

在x_i处对β₃进行泰勒级数展开，可得

$\begin{array}{*{20}{c}} {{\beta _3} = \\f{'}_i^2{h^2} + \left( {\frac{{13}}{{12}}f{''}_i^2 + \frac{1}{3}{{f'}_i}{{f'''}_i}} \right){h^4} + \frac{1}{{60}}{{f'}_i}f_i^{\left( 5 \right)}{h^6} + \frac{{13}}{{72}}{{f''}_i}f_i^{\left( 4 \right)}{h^6} + \frac{1}{{36}}f{'''}_i^2{h^6} + }\\ {\frac{1}{{2\;520}}{{f'}_i}f_i^{\left( 7 \right)}{h^8} + \frac{{13}}{{2\;160}}{{f''}_i}f_i^{\left( 6 \right)}{h^8} + \frac{1}{{360}}{{f'''}_i}f_i^{\left( 5 \right)}{h^8} + \frac{{13}}{{1\;728}}f_i^{\left( 4 \right)2}{h^8} + o\left( {{h^{10}}} \right)} \end{array}$

(23)

2.3 改进的三阶WENO-Z格式

在三阶WENO-NP3构造思想启发下，本研究提出如下的构造方法(通过对权函数中的光滑因子进行非线性处理)，并采用泰勒级数展开的方式进行理论推导，给出最终的构造格式

${\alpha _k} = {d_k}\left[ {1 + \frac{{{\tau _{\rm{N}}}}}{{{{\left( {{\beta _k} + \varepsilon } \right)}^p}}}} \right],\;\;\;\;k = 0,1$

(24)

${\tau _{\rm{N}}} = \left| {\frac{{{\beta _0} + {\beta _1}}}{2} - {\beta _3}} \right|,\;\;\;\;\varepsilon = {10^{ - 40}}$

(25)

在这里详细给出右界面通量f_i+1/2所对应权函数ω_k⁺的计算过程，将(10)式中的光滑因子在x_i处进行泰勒级数展开，可得

$\left\{ \begin{array}{l} {\beta _0} = f{'}_i^2{h^2} - {{f'}_i}{{f''}_i}{h^3} + \left( {\frac{1}{4}f{''}_i^2 + \frac{1}{3}{{f'}_i}{{f'''}_i}} \right){h^4} - \frac{1}{{12}}{{f'}_i}f_i^{\left( 4 \right)}{h^5} -\\ \frac{1}{6}{{f''}_i}{{f'''}_i}{h^5} + \frac{1}{{60}}{{f'}_i}f_i^{\left( 5 \right)}{h^6} + \frac{1}{{24}}{{f''}_i}f_i^{\left( 4 \right)}{h^6} + \frac{1}{{36}}f{'''}_i^2{h^6} - \frac{1}{{360}}{{f'}_i}f_i^{\left( 6 \right)}{h^7} - \\\frac{1}{{120}}{{f''}_i}f_i^{\left( 5 \right)}{h^7} - \frac{1}{{72}}{{f'''}_i}f_i^{\left( 4 \right)}{h^7} + \frac{1}{{2\;520}}{{f'}_i}f_i^{\left( 7 \right)}{h^8} + \frac{1}{{720}}{{f''}_i}f_i^{\left( 6 \right)}{h^8} +\\ \frac{1}{{360}}{{f'''}_i}f_i^{\left( 5 \right)}{h^8} + \frac{1}{{576}}f_i^{\left( 4 \right)2}{h^8} + o\left( {{h^9}} \right)\\ {\beta _1} = f{'}_i^2{h^2} - {{f'}_i}{{f''}_i}{h^3} + \left( {\frac{1}{4}f{''}_i^2 + \frac{1}{3}{{f'}_i}{{f'''}_i}} \right){h^4} + \frac{1}{{12}}{{f'}_i}f_i^{\left( 4 \right)}{h^5} +\\ \frac{1}{6}{{f''}_i}{{f'''}_i}{h^5} + \frac{1}{{60}}{{f'}_i}f_i^{\left( 5 \right)}{h^6} + \frac{1}{{24}}{{f''}_i}f_i^{\left( 4 \right)}{h^6} + \frac{1}{{36}}f{'''}_i^2{h^6} + \frac{1}{{360}}{{f'}_i}f_i^{\left( 6 \right)}{h^7} + \\\frac{1}{{120}}{{f''}_i}f_i^{\left( 5 \right)}{h^7} + \frac{1}{{72}}{{f''}_i}f_i^{\left( 4 \right)}{h^7} + \frac{1}{{2\;520}}{{f'}_i}f_i^{\left( 7 \right)}{h^8} +\frac{1}{{720}}{{f''}_i}f_i^{\left( 7 \right)}{h^8} + \\\frac{1}{{360}}{{f'''}_i}f_i^{\left( 6 \right)}{h^8} + \frac{1}{{576}}f_i^{\left( 4 \right)2}{h^8} + o\left( {{h^9}} \right) \end{array} \right.$

(26)

根据(23)式、(26)式可知，在一阶极值点处( ${f'_i} = 0, {f''_i} \ne 0, {f'''_i} \ne 0$ )，可得

$\left\{ \begin{array}{l} {\beta _0} = \frac{1}{4}f{''}_i^2{h^4} - \frac{{{{f''}_i}{{f'''}_i}}}{6}{h^5} + o\left( {{h^6}} \right)\\ {\beta _1} = \frac{1}{4}f{''}_i^2{h^4} + \frac{{{{f''}_i}{{f'''}_i}}}{6}{h^5} + o\left( {{h^6}} \right) \end{array} \right.$

(27)

$\left| {\frac{{{\beta _0} + {\beta _1}}}{2} - {\beta _3}} \right| = \frac{5}{6}f{''}_i^2{h^4} + \frac{5}{{36}}{{f''}_i}f_i^{\left( 4 \right)}{h^6} + o\left( {{h^8}} \right)$

(28)

将(27)式、(28)式代入(24)式中，可得

$\alpha _0^ + = {d_0}\left( {1 + \frac{{{\tau _{\rm{N}}}}}{{\beta _0^p}}} \right) = \\{d_0}\left[ {1 + \frac{5}{6}{4^p}{{\left( {{{f''}_i}} \right)}^{2 - 2p}}{h^{4 - 4p}} + \frac{{5p}}{9}{4^p}{{\left( {{{f''}_i}} \right)}^{1 - 2p}}{{f'''}_i}{h^{5 - 4p}} + o\left( {{h^{6 - 4p}}} \right)} \right]$

(29)

$\alpha _1^ + = {d_1}\left( {1 + \frac{{{\tau _{\rm{N}}}}}{{\beta _1^p}}} \right) = \\{d_1}\left[ {1 + \frac{5}{6}{4^p}{{\left( {{{f''}_i}} \right)}^{2 - 2p}}{h^{4 - 4p}} - \frac{{5p}}{9}{4^p}{{\left( {{{f''}_i}} \right)}^{1 - 2p}}{{f'''}_i}{h^{5 - 4p}} + o\left( {{h^{6 - 4p}}} \right)} \right]$

(30)

再根据(9)式中的加权法则，可得

$\omega _0^ + = \frac{{\alpha _0^ + }}{{\alpha _0^ + + \alpha _1^ + }} = \frac{1}{3} + \frac{{20p}}{{81}}{4^p}{\left( {{{f''}_i}} \right)^{1 - 2p}}{{f'''}_i}{h^{5 - 4p}} + o\left( {{h^{6 - 4p}}} \right)$

(31)

$\omega _1^ + = \frac{{\alpha _1^ + }}{{\alpha _0^ + + \alpha _1^ + }} = \frac{2}{3} - \frac{{20p}}{{81}}{4^p}{\left( {{{f''}_i}} \right)^{1 - 2p}}{{f'''}_i}{h^{5 - 4p}} + o\left( {{h^{6 - 4p}}} \right)$

(32)

同理可得，左界面通量f_i－1/2所对应权函数ω_k^－

$\omega _0^ - = \frac{1}{3} + \frac{{40p}}{{243}}{\left( {\frac{4}{9}} \right)^p}{\left( {{{f''}_i}} \right)^{1 - 2p}}{{f'''}_i}{h^{5 - 4p}} + o\left( {{h^{6 - 4p}}} \right)$

(33)

$\omega _0^ - = \frac{2}{3} - \frac{{40p}}{{243}}{\left( {\frac{4}{9}} \right)^p}{\left( {{{f''}_i}} \right)^{1 - 2p}}{{f'''}_i}{h^{5 - 4p}} + o\left( {{h^{6 - 4p}}} \right)$

(34)

根据(31)式~(34)式，可知

$\omega _k^ \pm - {d_k} = o\left( {{h^{5 - 4p}}} \right)$

(35)

再根据(17)式给出的充分条件可得参数p的数值

$5 - 4p = 2 \Rightarrow p = \frac{3}{4}$

(36)

并将该格式命名为WENO-NN3格式。

3. 数值试验

为了考察改进格式WENO-NN3格式的计算性能，选取线性精度测试、一维Sod激波管、激波与熵波相互作用、Rayleigh-Taylor不稳定性等经典算例进行自主编程计算，并将该格式计算结果与格式WENO-JS3、WENO-Z3和WENO-N3进行对比。

3.1 精度测试

该算例选自文献[8]，计算初始条件为

${u_0}\left( x \right) = \sin \left[ {{\rm{ \mathsf{ π} }}x - \frac{{\sin \left( {{\rm{ \mathsf{ π} }}x} \right)}}{{\rm{ \mathsf{ π} }}}} \right]$

(37)

其包含两个一阶极值点。表 1给出WENO-JS3、WENO-Z3、WENO-N3和WENO-NN3格式的L₁误差和精度。

表 1 针对初始条件(37)式在计算时间t=2时不同数值计算格式L₁误差和精度比较

Table 1. A comparative study of L₁ (error and order) for different schemes with initial condition Eq.(37) at t=2

N	WENO-JS3		WENO-Z3		WENO-N3		WENO-NN3
N	L₁ (Error)	L₁ (Order)	L₁ (Error)	L₁ (Order)	L₁ (Error)	L₁ (Order)	L₁ (Error)	L₁ (Order)
25	1.242 9×10^-1		5.727 6×10^-2		4.825 4×10^-2		2.595 4×10^-2
50	4.605 0×10^-2	1.432 4	1.511 6×10^-2	1.921 9	1.131 5×10^-2	2.092 4	3.796 1×10^-3	2.773 4
100	1.301 0×10^-2	1.823 6	3.464 5×10^-3	2.125 4	2.554 9×10^-3	2.146 9	4.660 5×10^-4	3.026 0
200	1.513 2×10^-3	3.103 9	8.372 1×10^-4	2.049 0	4.980 9×10^-4	2.358 8	4.865 7×10^-5	3.259 8
400	1.461 8×10^-4	3.371 8	1.682 3×10^-4	2.315 2	1.006 7×10^-4	2.306 8	6.294 1×10^-6	2.950 6

下载: 导出CSV

| 显示表格

根据表 1可知，本研究提出的改进格式WENO-NN3在极值点处能达到三阶设计精度，与2.3节中的理论分析保持一致。

3.2 Sod激波管

该算例初始条件^[22]为

$\left( {\rho ,u,p} \right) = \left\{ \begin{array}{l} \left( {1,0,1} \right)\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;0 \le x < 0.5\\ \left( {0.125,0,0.1} \right)\;\;\;\;\;\;\;0.5 \le x \le 1 \end{array} \right.$

(38)

网格数为200，计算结束时间为0.18。图 1给出该算例计算结束时刻密度曲线图及局部放大图。

图 1 Sod激波管计算结束后密度曲线及其局部放大图

Figure 1. Density curve and partially enlarged detail at the final time for the Sod problem

下载: 全尺寸图片幻灯片

3.3 激波与熵波相互作用

该算例初始条件^[20]为

$\left( {\rho ,u,p} \right) = \left\{ \begin{array}{l} \left( {3.857\;143,2.629\;369,10.333\;33} \right)\;\;\;\;\; - 5 \le x < - 4\\ \left( {1 + 0.2\sin \left( {5x} \right),0.1} \right)\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; - 4 \le x \le 5 \end{array} \right.$

(39)

网格数为800，计算结束时间为1.8。图 2给出了计算结束时刻密度曲线图及局部放大图。

图 2 激波与熵波相互作用计算结束后密度曲线及其局部放大图

Figure 2. Density curve and partially enlarged detail at the final time for the Shu-Osher problem

下载: 全尺寸图片幻灯片

根据上述两个一维经典算例计算结果可知，改进格式WENO-NN3计算性能最优，相较其他格式(WENO-JS3、WENO-Z3和WENO-N3)不仅具有更低耗散性，同时具有较高的精度。

3.4 Rayleigh-Taylor不稳定性

该问题主要描述重力场作用下，重流体加速进入轻流体界面失稳过程。文献[19, 23-24]也采用该算例探讨数值方法的分辨率特性。计算域设置为[0, 0.25]×[0, 1], 计算初始条件为

$\left( {\rho ,u,v,p} \right) = \left\{ \begin{array}{l} \left( {2,0, - 0.025\sqrt {\gamma p/\rho } \cos \left( {8{\rm{ \mathsf{ π} }}x} \right),2y + 1} \right)\;\;\;\;\;\;0 \le y < 0.5\\ \left( {1,0, - 0.025\sqrt {\gamma p/\rho } \cos \left( {8{\rm{ \mathsf{ π} }}x} \right),y + 1.5} \right)\;\;\;\;\;0.5 \le y < 1 \end{array} \right.$

(40)

该算例中，绝热指数γ取为5/3, 左右边界设置成反射边界条件，顶部和底部边界条件分别为(ρ, u, v, p)=(1, 0, 0, 2.5)，(ρ, u, v, p)=(2, 0, 0, 1)。计算结束时间为1.95。

图 3给出不同格式(WENO-JS3、WENO-Z3、WENO-N3和WENO-NN3)密度曲线图，网格数划分为240×960，共绘制15条等值线，其取值区间为[0.952 269，2.145 89]。从图 3接触间断附近精细结构的计算结果可以明显看出，改进格式WENO-NN3具有更好的分辨率特性。

图 3 Rayleigh-Taylor不稳定性问题不同格式(WENO-JS3、WENO-Z3、WENO-N3、WENO-NN3)密度曲线图

Figure 3. Density contours of the Rayleigh-Taylor instability computed using WENO-JS3, WENO-Z3, WENO-N3, and WENO-NN3 schemes

下载: 全尺寸图片幻灯片

3.5 Richtmyer-Meshkov不稳定性

该算例选自文献[25]，模拟激波管内部马赫数为2.0的平面激波冲击三模正弦形式空气/氦气界面产生的不稳定现象。界面左端是空气介质，界面右端是氦气。计算域取为[0, 0.6]×[0, 0.1]，各介质的初始参数为

$\left\{ \begin{array}{l} {\left( {\rho ,u,v,p} \right)_{\rm{L}}} = \left( {2.\;67,1.\;48,0,4.\;5} \right)\\ {\left( {\rho ,u,v,p} \right)_{\rm{M}}} = \left( {1,0,0,1} \right)\\ {\left( {\rho ,u,v,p} \right)_{\rm{R}}} = \left( {0.\;138,0,0,1} \right) \end{array} \right.$

(41)

平面激波初始位置位于x=0.06，三模界面的初始位置见(42)式，初场分布见图 4。

$x = 0.1 + 0.\;008\cos \left( {60{\rm{ \mathsf{ π} }}y} \right)$

(42)

图 4 三模Richtmyer-Meshkov不稳定性初始条件设置

Figure 4. Initial condition of the treble-mode Richtmyer-Meshkov instability

下载: 全尺寸图片幻灯片

上下边界条件设置为固壁边界，左边界设置为反射边界条件，右边界设置为流出边界条件。计算网格数为1 200×200，计算结束时间为17.22×10^-2。图 5给出了不同格式(WENO-JS3、WENO-Z3、WENO-N3和WENO-NN3)计算的密度曲线图，从图中尖钉结构可以明显看出，改进格式WENO-NN3格式具有更低的耗散性，描述出更精细的流场结构。

图 5 三模Richtmyer-Meshkov不稳定性问题不同格式(WENO-JS3、WENO-Z3、WENO-N3、WENO-NN3)密度曲线图

Figure 5. Density contours of the treble-mode Richtmyer-Meshkov instability computed using WENO-JS3, WENO-Z3, WENO-N3, and WENO-NN3 schemes

下载: 全尺寸图片幻灯片

4. 结论

在理论分析的基础上，提出一种提高极值点计算精度的改进三阶WENO-Z格式，并通过精度测试、Sod激波管等经典算例，从精度和耗散性两个角度考察了改进格式WENO-NN3的计算性能。

(1) 按照本研究的理论分析及数值精度测试，改进格式WENO-NN3格式能在极值点处达到收敛精度。推导的三阶WENO格式收敛精度的充分条件及采用泰勒展开的分析方法可以扩展应用到高阶格式。

(2) 改进格式WENO-NN3格式相较其他格式(WENO-JS3、WENO-Z3和WENO-N3)具有低耗散、高分辨率的特性。

图 1 Sod激波管计算结束后密度曲线及其局部放大图

Figure 1. Density curve and partially enlarged detail at the final time for the Sod problem

下载: 全尺寸图片幻灯片

图 2 激波与熵波相互作用计算结束后密度曲线及其局部放大图

Figure 2. Density curve and partially enlarged detail at the final time for the Shu-Osher problem

下载: 全尺寸图片幻灯片

图 3 Rayleigh-Taylor不稳定性问题不同格式(WENO-JS3、WENO-Z3、WENO-N3、WENO-NN3)密度曲线图

Figure 3. Density contours of the Rayleigh-Taylor instability computed using WENO-JS3, WENO-Z3, WENO-N3, and WENO-NN3 schemes

下载: 全尺寸图片幻灯片

图 4 三模Richtmyer-Meshkov不稳定性初始条件设置

Figure 4. Initial condition of the treble-mode Richtmyer-Meshkov instability

下载: 全尺寸图片幻灯片

图 5 三模Richtmyer-Meshkov不稳定性问题不同格式(WENO-JS3、WENO-Z3、WENO-N3、WENO-NN3)密度曲线图

Figure 5. Density contours of the treble-mode Richtmyer-Meshkov instability computed using WENO-JS3, WENO-Z3, WENO-N3, and WENO-NN3 schemes

下载: 全尺寸图片幻灯片

表 1 针对初始条件(37)式在计算时间t=2时不同数值计算格式L₁误差和精度比较

Table 1. A comparative study of L₁ (error and order) for different schemes with initial condition Eq.(37) at t=2

N	WENO-JS3		WENO-Z3		WENO-N3		WENO-NN3
N	L₁ (Error)	L₁ (Order)	L₁ (Error)	L₁ (Order)	L₁ (Error)	L₁ (Order)	L₁ (Error)	L₁ (Order)
25	1.242 9×10^-1		5.727 6×10^-2		4.825 4×10^-2		2.595 4×10^-2
50	4.605 0×10^-2	1.432 4	1.511 6×10^-2	1.921 9	1.131 5×10^-2	2.092 4	3.796 1×10^-3	2.773 4
100	1.301 0×10^-2	1.823 6	3.464 5×10^-3	2.125 4	2.554 9×10^-3	2.146 9	4.660 5×10^-4	3.026 0
200	1.513 2×10^-3	3.103 9	8.372 1×10^-4	2.049 0	4.980 9×10^-4	2.358 8	4.865 7×10^-5	3.259 8
400	1.461 8×10^-4	3.371 8	1.682 3×10^-4	2.315 2	1.006 7×10^-4	2.306 8	6.294 1×10^-6	2.950 6

下载: 导出CSV

参考文献(25)

[1]	LIU X D, OSHER S, CHAN T.Weighted essentially non-oscillatory schemes[J].Journal of Computational Physics, 1994, 115(1):200-212. doi: 10.1006/jcph.1994.1187
[2]	HARTEN A, ENGQUIST B, OSHER S, et al. Uniformly high order accurate essentially non-oscillatory schemes (Ⅲ)//HUSSAINI M Y, VAN LEER B, VAN ROSENDALE J. Upwind and high-resolution schemes. Berlin, Heidelberg: Springer, 1987: 231-303.
[3]	JIANG G S, SHU C W.Efficient implementation of weighted ENO schemes[J].Journal of Computational Physics, 1995, 126(1):202-228.
[4]	HSIEH T J, WANG C H, YANG J Y.Numerical experiments with several variant WENO schemes for the Euler equations[J].International Journal for Numerical Methods in Fluids, 2010, 58(9):1017-1039.
[5]	ZHAO S, LARDJANE N, FEDIOUN I.Comparison of improved finite-difference WENO schemes for the implicit large eddy simulation of turbulent non-reacting and reacting high-speed shear flows[J].Computers & Fluids, 2014, 95(3):74-87. https://www.sciencedirect.com/science/article/pii/S0045793014000802
[6]	WANG C, DING J X, SHU C W, LI T.Three-dimensional ghost-fluid large-scale numerical investigation on air explosion[J].Computers & Fluids, 2016, 137:70-79. https://www.sciencedirect.com/science/article/pii/S0045793016302353
[7]	ZAGHI S, MASCIO A D, FAVINI B.Application of WENO-positivity-preserving schemes to highly under-expanded jets[J].Journal of Scientific Computing, 2016, 69(3):1-25. doi: 10.1007/s10915-016-0226-5.pdf
[8]	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. doi: 10.1016/j.jcp.2005.01.023
[9]	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. doi: 10.1016/j.jcp.2007.11.038
[10]	YAMALEEV N K, CARPENTER M H.A systematic methodology for constructing high-order energy stable WENO schemes[J].Journal of Computational Physics, 2009, 228:4248-4272. doi: 10.1016/j.jcp.2009.03.002
[11]	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. doi: 10.1016/j.jcp.2010.11.028
[12]	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(1):68-86. doi: 10.1016/j.jcp.2012.06.016
[13]	SHEN Y Q, ZHA G C.Improvement of weighted essentially non-oscillatory schemes near discontinuities[J].Computers & Fluids, 2014, 96(12):1-9. https://www.sciencedirect.com/science/article/pii/S0045793014000656
[14]	CHANG H K, HA Y, YOON J.Modified non-linear weights for fifth-order weighted essentially non-oscillatory schemes[J].Journal of Computational Science, 2016, 67(1):299-323. doi: 10.1007/s10915-015-0079-3
[15]	YAMALEEV N K, CARPENTER M H.Third-order energy stable WENO scheme[J].Journal of Computational Physics, 2013, 228(8):3025-3047. https://www.sciencedirect.com/science/article/pii/S002199910900014X
[16]	WU X S, ZHAO Y X.A high-resolution hybrid scheme for hyperbolic conservation laws[J].International Journal for Numerical Methods in Fluids, 2015, 78(3):162-187. doi: 10.1002/FLD.v78.3
[17]	WU X, LIANG J, ZHAO Y.A new smoothness indicator for third-order WENO scheme[J].International Journal for Numerical Methods in Fluids, 2016, 81(7):451-459. doi: 10.1002/fld.v81.7
[18]	DON W S, BORGES R.Accuracy of the weighted essentially non-oscillatory conservative finite difference schemes[J].Journal of Computational Physics, 2013, 250(4):347-372. https://www.sciencedirect.com/science/article/pii/S0021999113003501
[19]	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. doi: 10.1016/j.jcp.2010.08.019
[20]	SHU C W, OSHER S.Efficient implementation of essentially non-oscillatory shock-capturing schemes, Ⅱ[J].Journal of Computational Physics, 1989, 77(2):439-471. doi: 10.1007/978-3-642-60543-7_14.pdf
[21]	GANDE N R, RATHOD Y, RATHAN S.Third-order WENO scheme with a new smoothness indicator[J].International Journal for Numerical Methods in Fluids, 2017, 85(2):171-185. https://www.researchgate.net/publication/314137535_Third_order_WENO_scheme_with_a_new_smoothness_indicator
[22]	SOD G A.A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws[J].Journal of Computational Physics, 1978, 27(1):1-31. https://www.sciencedirect.com/science/article/pii/0021999178900232
[23]	SHI J, ZHANG Y T, SHU C W.Resolution of high order WENO schemes for complicated flow structures[J].Journal of Computational Physics, 2003, 186(2):690-696. doi: 10.1016/S0021-9991(03)00094-9
[24]	ACKER F, BORGES R, COSTA B.An improved WENO-Z scheme[J].Journal of Computational Physics, 2016, 313:726-753. doi: 10.1016/j.jcp.2016.01.038
[25]	ZHANG P G, WANG J P.A newly improved WENO scheme and its application to the simulation of Richtmyer-Meshkov instability[J].Procedia Engineering, 2013, 61:325-332. doi: 10.1016/j.proeng.2013.08.023

施引文献

资源附件(0)

访问统计

点击查看大图

图(5) / 表(1)

计量

文章访问数: 7680
HTML全文浏览量: 3698
PDF下载量: 214

1. 控制方程
2. 数值方法
2.1 三阶WENO-JS格式
2.2 三阶WENO-NP3格式
2.3 改进的三阶WENO-Z格式
3. 数值试验
3.1 精度测试
3.2 Sod激波管
3.3 激波与熵波相互作用
3.4 Rayleigh-Taylor不稳定性
3.5 Richtmyer-Meshkov不稳定性
4. 结论

一种提高极值点处精度的三阶WENO-Z改进格式及应用

doi: 10.11858/gywlxb.20170696

作者简介: 徐维铮(1991-), 男, 博士, 主要从事爆炸波数值计算方法及程序开发研究.E-mail:xuweizheng@whut.edu.cn

通讯作者: 吴卫国(1960-), 男, 教授, 博士生导师, 主要从事结构动力学及舰船抗爆抗冲击研究.E-mail:mailjt@163.com

计量

出版历程

An Improved Third-Order WENO-Z Scheme for Achieving Optimal Order near Critical Points and Its Application

1. 控制方程

2. 数值方法

2.1 三阶WENO-JS格式

2.1.1 数值重构过程

2.1.2 收敛精度的充分条件

2.2 三阶WENO-NP3格式

2.3 改进的三阶WENO-Z格式

3. 数值试验

3.1 精度测试

3.2 Sod激波管

3.3 激波与熵波相互作用

3.4 Rayleigh-Taylor不稳定性

3.5 Richtmyer-Meshkov不稳定性

4. 结论

计量

出版历程

目录

1. 控制方程

2. 数值方法

2.1 三阶WENO-JS格式

2.2 三阶WENO-NP3格式

2.3 改进的三阶WENO-Z格式

3. 数值试验

3.1 精度测试

3.2 Sod激波管

3.3 激波与熵波相互作用

3.4 Rayleigh-Taylor不稳定性

3.5 Richtmyer-Meshkov不稳定性

4. 结论

作者简介:
徐维铮(1991-), 男, 博士, 主要从事爆炸波数值计算方法及程序开发研究.E-mail:xuweizheng@whut.edu.cn

通讯作者:
吴卫国(1960-), 男, 教授, 博士生导师, 主要从事结构动力学及舰船抗爆抗冲击研究.E-mail:mailjt@163.com