﻿ Godunov型显式大时间步长格式研究进展
 文章快速检索 高级检索
Godunov型显式大时间步长格式研究进展

1. 中国航空工业空气动力研究院, 沈阳 110034;
2. 高速高雷诺数气动力航空科技重点实验室, 沈阳 110034

Research progress of Godunov type explicit large time step scheme
QIAN Zhansen1,2
1. AVIC Aerodynamics Research Institute, Shenyang 110034, China;
2. Aviation Key Laboratory of Science and Technology on High Speed and High Reynolds Number Aerodynamic Force Research, Shenyang 110034, China
Abstract: The research progress of the Godunov type explicit large time step scheme is reviewed. Firstly, the concept, classifications and advantages of explicit large time step scheme are introduced. Then, followed by its construction methods, higher order accuracy extension approaches, multi-dimension generalization methods, and characteristics analysis including stability, resolution and computational efficiency. Finally, further development direction of the Godunov type explicit large time step scheme is proposed.
Keywords: large time step schemes    Godunov type schemes    highly efficient algorithm    scheme stability    explicit schemes

1 显式大时间步长格式简介 1.1 数值格式的发展概述

1) 双曲型守恒律的时间相关方法

2) 激波捕捉格式的发展

1.2 显式大时间步长格式的概念

1.3 大时间步长格式的分类

1.4 大时间步长格式的优势

1) 定常问题计算效率高

2) 可直接应用于非定常问题计算

3) 对间断分辨率高

4) 大规模并行计算可扩展性好

5) 可与现有加速收敛的算法良好兼容

2 Godunov型大时间步长格式的提出 2.1 Godunov型格式

 $\left\{ {\begin{array}{*{20}{l}} {\frac{{\partial \mathit{\boldsymbol{u}}}}{{\partial t}} + \frac{{\partial f(\mathit{\boldsymbol{u}})}}{{\partial x}} = 0\quad - \infty < x < + \infty ,t \ge 0}\\ {\mathit{\boldsymbol{u}}(0,x) = {\mathit{\boldsymbol{u}}^0}(x)} \end{array}} \right.$ （1）

1) 第1种方法直接在区间[xj-1/2, xj+1/2]上求积分，如图 1所示，在CFL≤0.5的条件下，$\mathit{\boldsymbol{u}}_j^{n + 1}$可表达为

 $\begin{array}{*{20}{l}} {\mathit{\boldsymbol{u}}_j^{n + 1} = }\\ {\quad \frac{1}{{\Delta x}}\left[ {\int_{{x_{j - 1/2}}}^{{x_j}} \mathit{\boldsymbol{u}} (x,{t_{n + 1}}){\rm{d}}x + \int_{{x_j}}^{{x_{j + 1/2}}} \mathit{\boldsymbol{u}} (x,{t_{n + 1}}){\rm{d}}x} \right]} \end{array}$ （2）
 图 1 局部Godunov平均示意图 Fig. 1 Schematic of local Godunov averaging

2) 第2种方法采用守恒型格式：

 $\mathit{\boldsymbol{u}}_j^{n + 1} = \mathit{\boldsymbol{u}}_j^n - \frac{{\Delta t}}{{\Delta x}}({\mathit{\boldsymbol{\hat f}}_{j + 1/2}} - {\mathit{\boldsymbol{\hat f}}_{j - 1/2}})$ （3）

 ${\mathit{\boldsymbol{\hat f}}_{j + 1/2}} = \frac{1}{{\Delta t}}\int_{{t_n}}^{{t_{n + 1}}} {\mathit{\boldsymbol{f}}(\mathit{\boldsymbol{u}}({x_{j + 1/2}},t)} ){\rm{d}}t$ （4）

2.2 Godunov型大时间步长格式

 $\left\{ {\begin{array}{*{20}{l}} {\frac{{\partial u}}{{\partial t}} + \frac{{\partial f(u)}}{{\partial x}} = 0\quad {x_j} < x < {x_{j + 1}},t \ge 0}\\ {u_{j + 1/2}^n(x,{t_n}) = \left\{ {\begin{array}{*{20}{l}} {u_j^n}&{x < {x_{j + 1/2}}}\\ {u_{j + 1}^n}&{x \ge {x_{j + 1/2}}} \end{array}} \right.} \end{array}} \right.$ （5）

 $\left\{ {\begin{array}{*{20}{l}} {\frac{{\partial u}}{{\partial t}} + \frac{{\partial f(u)}}{{\partial x}} = 0}&{ - \infty < x < + \infty ,t \ge 0}\\ {{u^n}(x,{t_n}) = u_j^n}&{{x_{j - 1/2}} < x < {x_{j + 1/2}}} \end{array}} \right.$ （6）

 $\begin{array}{l} {{\tilde u}^n}(x,t) = {u^n}(x,{t_n}) + \sum\limits_j {(u_{j + 1/2}^n(} x,t) - \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} u_j^n(x,{t_n})) \end{array}$ （7）

 $\frac{{\varepsilon - {x_{j + 1/2}}}}{{\Delta x}}\Delta u$ （8）
 图 2 单波束在网格上的传播[84] Fig. 2 A single wave propagation on grid[84]

 图 3 多波束在网格上的传播[86] Fig. 3 Multi-wave propagation on grid[86]
 $\begin{array}{*{20}{l}} {u_j^{n + 1} = u_j^n - \bar \Delta {u_{j - 1/2}} - \cdots - \bar \Delta {u_{j - i + \frac{1}{2}}} - \cdots - }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \bar \Delta {u_{j - K + 1/2}} + \bar \Delta {u_{j + 1/2}} + \cdots + \bar \Delta {u_{j + i - 1/2}} + \cdots + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \bar \Delta {u_{j + K - 1/2}}} \end{array}$ （9）

 $\bar \Delta {u_{j - i + 1/2}} = \left\{ {\begin{array}{*{20}{l}} 0&{{s_{\rm{a}}} \le (i - 1)\Delta x/\Delta t}\\ {({s_{\rm{a}}}\Delta t/\Delta x - i + 1)({u_{j - i + 1}} - {u_{j - i}})}&{(i - 1)\Delta x/\Delta t < {s_{\rm{a}}} < i\Delta x/\Delta t}\\ {{u_{j - i + 1}} - {u_{j - i}}}&{i\Delta x/\Delta t \le {s_{\rm{a}}}} \end{array}} \right.$ （10）

LeVeque[80, 82]的数值实验表明，以上格式可以突破CFL≤1的限制，且在理论上CFL数可以取无限大。但实际上，由于线化假设等原因，CFL数只能取到有限值。

 图 4 Euler方程多波束在网格上的传播[86] Fig. 4 Multi-wave propagation on grid for Euler equations[86]

 $\begin{array}{l} u_j^{n + 1} = u_j^n - (\bar \Delta u_{j - 1/2}^1 + \cdots + \bar \Delta \mathit{\boldsymbol{u}}_{j - 1/2}^m + \cdots + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \bar \Delta u_{j - 1/2}^M) - \cdots - (\bar \Delta u_{j - i - 1/2}^1 + \cdots + \bar \Delta u_{j - i + 1/2}^m + \cdots + \\ \begin{array}{*{20}{l}} {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \bar \Delta u_{j - i + 1/2}^M) - \cdots - (\bar \Delta u_{j - K + 1/2}^1 + \cdots + \bar \Delta u_{j - K + 1/2}^m + \cdots + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \bar \Delta u_{j - K + 1/2}^M) + (\bar \Delta u_{j + 1/2}^1 + \cdots + \bar \Delta u_{j + 1/2}^m + \cdots + } \end{array}\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \begin{array}{*{20}{l}} {\bar \Delta u_{j + 1/2}^M) + \cdots + (\bar \Delta u_{j + i - 1/2}^1 + \cdots + \bar \Delta u_{j + i - 1/2}^m + \cdots + }\\ {\bar \Delta u_{j + i - 1/2}^M) + \cdots + (\bar \Delta u_{j + K - 1/2}^1 + \cdots + }\\ {\left. {\bar \Delta u_{j + K - 1/2}^m + \cdots + \bar \Delta u_{j + K - 1/2}^M} \right)} \end{array} \end{array}$ （11）

 $\bar \Delta u_{j - i + 1/2}^m = \left\{ {\begin{array}{*{20}{l}} 0&{{s_{{\rm{a}},m}} \le (i - 1)\Delta x/\Delta t}\\ {({s_{{\rm{a}},m}}\Delta t/\Delta x - i + 1) \cdot (\mathit{\boldsymbol{u}}_m^ + - \mathit{\boldsymbol{u}}_m^ - )}&{(i - 1)\Delta x/\Delta t < {{\left\| \mathit{\boldsymbol{u}} \right\|}_M} < i\Delta x/\Delta t}\\ {\mathit{\boldsymbol{u}}_m^ + - \mathit{\boldsymbol{u}}_m^ - }&{{s_{{\rm{a}},m}} \ge i\Delta x/\Delta t} \end{array}} \right.$ （12）

3 Godunov型大时间步长格式的改进

LeVeque[80-81]提出的大时间步长格式主要针对标量双曲型守恒律方程，自其提出该思路以来，研究者一直试图将其拓展到复杂问题的求解中，主要的应用包括气体动力学Euler方程和水动力学浅水波方程两个方面。主要的改进也涉及两个方面，一是对非线性方程的改进处理方法，二是从标量方程到方程组的推广。但目前仍存在一些问题有待解决，包括强非线性、膨胀波、同族波系干扰等带来的问题。

3.1 时间线插值技术

u(xj+1/2, t)为处于xj+1/2界面上的Riemann问题的解，则界面数值通量也可表达为

 ${\mathit{\boldsymbol{\hat f}}_{j + 1/2}} = \mathit{\boldsymbol{f}}({\mathit{\boldsymbol{\hat u}}_{{x_{j + 1/2}}}})$ （13）

 ${\mathit{\boldsymbol{\hat u}}_{j + 1/2}} = \frac{1}{{\Delta t}}\int_{{t_n}}^{{t_{n + 1}}} \mathit{\boldsymbol{u}} ({x_{j + 1/2}},t){\rm{d}}t$ （14）

 $\mathit{\boldsymbol{u}} = \sum\limits_m {{\alpha ^{(m)}}} {\mathit{\boldsymbol{K}}^{(m)}}$ （15）

 ${\mathit{\boldsymbol{\hat u}}_{j + 1/2}} = \sum\limits_m {\frac{1}{{\Delta t}}} \int_{{t_n}}^{{t_{n + 1}}} {{\alpha ^{(m)}}} K_{j + 1/2}^{(m)}(t){\rm{d}}t$ （16）

 $\begin{array}{l} {{\mathit{\boldsymbol{\hat u}}}_{j + 1/2}} = \sum\limits_{m,{\lambda ^{(m)}} \ge 0} {\frac{1}{{\Delta t}}} \int_{{t_n}}^{{t_{n + 1}}} {{\alpha ^{(m)}}} K_j^{(m)}(t){\rm{d}}t + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \sum\limits_{m,{\lambda ^{(m)}} < 0} {\frac{1}{{\Delta t}}} \int_{{t_n}}^{{t_{n + 1}}} {{\alpha ^{(m)}}} K_{j + 1}^{(m)}(t){\rm{d}}t \end{array}$ （17）

 $\begin{array}{*{20}{l}} {{{\mathit{\boldsymbol{\hat u}}}_{j + 1/2}} = \sum\limits_{m,{\lambda ^{(m)}} \ge 0} = \frac{1}{{\Delta t}}\int_{{t_n}}^{{t_{n + 1}}} {{\alpha ^{(m)}}} K_{j + 1/2,{\rm{L}}}^{(m)}(t){\rm{d}}t + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \sum\limits_{m,{\lambda ^{(m)}} < 0} {\frac{1}{{\Delta t}}} \int_{{t_n}}^{{t_{n + 1}}} {{\alpha ^{(m)}}} K_{j + 1/2,{\rm{R}}}^{(m)}(t){\rm{d}}t} \end{array}$ （18）

 ${\mathit{\boldsymbol{K}}^{(m)}} = {\bf{0}},\frac{{{\rm{d}}x}}{{{\rm{d}}t}} = {\lambda ^{(m)}}$ （19）

 $\begin{array}{l} \frac{1}{{\Delta t}}\int_{{t_n}}^{{t_{n + 1}}} {{\alpha ^{(m)}}} K_{j + 1/2}^{(m)}(t){\rm{d}}t = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{1}{{{x_{j + 1/2}} - {x_{\rm{A}}}}}\int_{{x_{\rm{A}}}}^{{x_{j + 1/2}}} {\alpha _n^{(m)}} K_n^{(m)}(x){\rm{d}}x \end{array}$ （20）

 $\begin{array}{l} \frac{1}{{\Delta t}}\int_{{t_n}}^{{t_{n + 1}}} {\alpha _{j + 1/2}^{(m)}} K_{j + 1/2}^{(m)}(t){\rm{d}}t = \frac{1}{{\Delta t}}\int_{{t_n}}^\tau {\alpha _{j + 1/2}^{(m)}} K_{j + 1/2}^{(m)}(t){\rm{d}}t + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{1}{{\Delta t}}\int_\tau ^{{t_{n + 1}}} {\alpha _{j + 1/2}^{(m)}} K_{j + 1/2}^{(m)}(t){\rm{d}}t \end{array}$ （21）

 图 5 时间线插值法示意图[83] Fig. 5 Schematic of time-line interpolation method
 $\begin{array}{l} \frac{1}{{\tau - {t_n}}}\int_{{t_n}}^\tau {\alpha _{j + 1/2}^{(m)}} K_{j + 1/2}^{(m)}(t){\rm{d}}t = \frac{1}{{\Delta {x_j}}}\int_{{x_{j - 1/2}}}^{{x_{j + 1/2}}} {\alpha _n^{(m)}} \cdot \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} K_n^{(m)}(x){\rm{d}}x = \alpha _{j,n}^{(m)}K_{j,n}^{(m)} \end{array}$ （22）
 $\begin{array}{l} \frac{1}{{{t_{n + 1}} - \tau }}\int_\tau ^{{t_{n + 1}}} {\alpha _{j + 1/2}^{(m)}} K_{j + 1/2}^{(m)}(t){\rm{d}}t = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{1}{{{\tau ^\prime } - {t_n}}}\int_{{t_n}}^{{\tau ^\prime }} {\alpha _{j - 1/2,n}^{(m)}} K_{j - 1/2,n}^{(m)}(t){\rm{d}}t \end{array}$ （23）

 $\begin{array}{l} \alpha _{j - 1/2,{\rm{L}}}^{(m)}K_{j - 1/2,{\rm{L}}}^{(m)}(t) = \alpha _{j - 1/2,n/2,{\rm{L}}}^{(m)}K_{j - 1/2,n/2,{\rm{L}}}^{(m)} + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left( {t - {t_n} - \frac{{\Delta t}}{2}} \right)s_{j - 1/2,{\rm{L}}}^{(m)} \end{array}$ （24）

 $\begin{array}{*{20}{c}} {\alpha _{j + 1/2,n/2,{\rm{L}}}^{(m)}K_{j + 1/2,n/2,{\rm{L}}}^{(m)} = \frac{{\alpha _{j,n}^{(m)}K_{j,n}^{(m)}}}{{{\rm{CFL}}}} + \left( {1 - \frac{1}{{{\rm{CFL}}}}} \right) \cdot }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left( {\alpha _{j - 1/2,n/2,{\rm{L}}}^{(m)}K_{j - 1/2,n/2,{\rm{L}}}^{(m)} - \frac{{\Delta t}}{{2{\rm{CFL}}}}s_{j - 1/2,{\rm{L}}}^{(m)}} \right)} \end{array}$ （25）
 $\begin{array}{l} s_{j + 1/2,{\rm{L}}}^{(m)} = - \frac{2}{{\Delta t}}\alpha _{j,n}^{(m)}K_{j,n}^{(m)} + \frac{2}{{\Delta t}}\alpha _{j - 1/2,n/2,{\rm{L}}}^{(m)} \cdot \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} K_{j - 1/2,n/2,{\rm{L}}}^{(m)} - \frac{1}{{{\rm{CFL}}}}s_{j - 1/2,{\rm{L}}}^{(m)} \end{array}$ （26）

3.2 膨胀波的处理技术

2.2节提到的叠波法只能直接用于波系均为间断(激波或接触间断)的情形，在非线性条件下会出现膨胀波，而膨胀波在tn+1时刻有可能跨越若干个区间。LeVeque[82]曾提出采用一束非物理膨胀间断来代替膨胀波，并采用此法求解了一维等熵Euler方程。该方法计算效率高，但存在违背熵条件的非物理解的可能性。在叠波法中，如果对Riemann问题分解后出现的膨胀波处理不当，则在遇到较强膨胀波时会导致数值解恶化，甚至出现包含膨胀间断的非物理解。

 图 6 膨胀波的多波束近似示意图[84] Fig. 6 Schematic of multi-wave approximation of expansion wave[84]

 图 7 膨胀波网格单元分解法示意图[86] Fig. 7 Schematic of grid cell decomposition method for expansion wave[86]
3.3 波系干扰的处理技术

 图 8 波系干扰示意图[84] Fig. 8 Schematic of interaction of wave systems[84]

 图 9 两波碰撞示意图[84] Fig. 9 Schematic of interaction of two discontinuities[84]
3.4 近似Riemann求解器的应用

1) LTS-Roe格式

Lindqvist等[97]针对标量双曲型方程，采用单参数(One Parameter)格式，描述了多种大时间步长格式，其讨论了数值黏性和通量差分分裂两种表达形式的大时间步长格式，并探讨了其TVD性质。基于数值黏性形式，标量双曲型方程的大时间步长格式的数值通量可写为

 ${f_{j + 1/2}} = \frac{1}{2}({f_j} + {f_{j + 1}}) - \frac{1}{2}\sum\limits_{i = - \infty }^\infty {Q_{j + 1/2 + i}^i} \Delta {u_{j + 1/2 + i}}$ （27）

 $\begin{array}{*{20}{l}} {u_j^{n + 1} = u_j^n - \frac{{\Delta t}}{{\Delta x}}\sum\limits_{i = 0}^\infty ( A_{j - 1/2 - i}^{i + }\Delta {u_{j - 1/2 - i}} + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} A_{j + 1/2 + i}^{i - }\Delta {u_{j + 1/2 + i}})} \end{array}$ （28）

Lindqvist[97]等指出数值黏性和通量差分分裂两种形式的大时间步长格式是一一对应的，具有如下关系：

 ${{A^{0 \pm }} = \frac{1}{2}(\lambda \pm {Q^0} \mp {Q^{ \mp 1}})}$ （29a）
 ${{A^{i \pm }} = \pm \frac{1}{2}({Q^{ \mp i}} - {Q^{ \mp (i + 1)}})}$ （29b）

 ${Q^i} = \left\{ {\begin{array}{*{20}{l}} {2\sum\limits_{p = - i}^\infty {{A^{p + }}} }&{i < 0}\\ {2\sum\limits_{p = 0}^\infty {({A^{p + }} - {A^{p - }})} }&{i = 0}\\ { - 2\sum\limits_{p = - i}^\infty {{A^{p - }}} }&{i > 0} \end{array}} \right.$ （30）

2K+1点的LTS-Lax-Friedrichs格式的数值黏性系数为

 ${Q_{{\rm{LxF}}}^0 = K\frac{{\Delta x}}{{\Delta t}}}$ （31a）
 ${Q_{{\rm{LxF}}}^{ \mp i} = \frac{{K - i}}{K}\left( { \pm \lambda + K\frac{{\Delta x}}{{\Delta t}}} \right)\quad i > 0}$ （31b）

2K+1点的LTS-Roe格式的数值黏性系数为

 ${Q_{{\rm{Roe}}}^0 = |\lambda |}$ （32a）
 ${Q_{{\rm{Roe}}}^{ \mp i} = 2{\rm{max}}\left( {0, \pm \lambda - i\frac{{\Delta x}}{{\Delta t}}} \right)\quad i > 0}$ （32b）

Lindqvist[97]等指出LTS-Roe格式和LTS-Lax-Friedrichs格式分别为耗散最小和最大的具有TVD性质的大时间步长格式。LTS-Lax-Friedrichs格式数值耗散过大，在工程计算中不实用；LTS-Roe格式数值耗散在某些情形下不足，导致数值解在间断波附近出现非物理振荡，并有可能出现违背熵条件的可能，得到非物理的解。Lindqvist等结合了两种格式的优劣，提出了混合型的LTS-RoeLxF(β)格式：

 ${Q^\alpha } = \beta Q_{{\rm{LxF}}}^\alpha + (1 - \beta )Q_{{\rm{Roe}}}^\alpha$ （33）

Lindqvist等[97]采用局部线性近似方法，将基于标量双曲型方程建立的大时间步长格式推广至方程组情形，但是由于方程组的复杂性，在标量情形所讨论的有关格式的各种性质将不再全部有效。尽管如此，Lindqvist等仍采用一维Euler方程激波管问题作为数值算例，初步验证了所建立的方法。

2) LTS-HLL格式

Prebeg等[98]在Lindqvist等[97]研究工作的基础上提出了基于HLL[57]格式和HLLC格式[58]的大时间步长格式(LTS-HLL/HLLC)。其采用更广义的双参数(Two Parameters)描述形式，表达了包含Lindqvist等[97]提出的多种格式在内的大时间步长格式，重点探讨了基于HLL和HLLC格式的大时间步长格式LTS-HLL和LTS-HLLC格式。

HLL格式是Harten等[57]提出的一个新颖的方法，其通过近似求解Riemann问题来给出数值通量，该格式假设Riemann问题的解由两道波组成，并且通过估算给出两道波各自的传播速度，可快速得到数值通量，不像Roe格式和Osher格式那样详细给出Riemann问题解的波系结构。该格式比Roe格式计算效率更高，但是由于对Riemann问题的间断分解后的物理解结构采用了两波近似，使得格式对接触间断的分辨率下降。针对这一问题，Toro等[58]通过添加接触间断的影响，扩展两波近似为三波形式，给出了改近的HLL格式，称为HLLC格式，提高了对接触间断的分辨率。

Prebeg等[98]通过对标量双曲型守恒律引入两参数形式的数值解，得到了更广义的方法。LTS-HLL格式写成数值黏性形式为

 $Q_{{\rm{HLL}}}^0 = \frac{{|{S_{\rm{R}}}|(\lambda - {S_{\rm{L}}}) + |{S_{\rm{L}}}|({S_{\rm{R}}} - \lambda )}}{{{S_{\rm{R}}} - {S_{\rm{L}}}}}$ （34a）
 $\begin{array}{*{20}{l}} {Q_{{\rm{HLL}}}^{ \mp i} = 2\frac{{\lambda - {S_{\rm{L}}}}}{{{S_{\rm{R}}} - {S_{\rm{L}}}}}{\rm{max}}\left( {0, \pm {S_{\rm{R}}} - i\frac{{\Delta x}}{{\Delta t}}} \right) + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 2\frac{{{S_{\rm{R}}} - \lambda }}{{{S_{\rm{R}}} - {S_{\rm{L}}}}}{\rm{max}}\left( {0, \pm {S_{\rm{L}}} - i\frac{{\Delta x}}{{\Delta t}}} \right)\quad i > 0} \end{array}$ （34b）

LTS-HLL格式也可以改写成叠波形式

 $\begin{array}{*{20}{l}} {u_j^{n + 1} = u_j^n - \frac{{\Delta t}}{{\Delta x}}\sum\limits_{i = 0}^\infty {\left( {\sum\limits_{p = 1}^2 {S_{j - 1/2 - i}^{p,i + }} w_{j - 1/2 - i}^p + } \right.} }\\ {\left. {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \sum\limits_{b = 1}^2 {S_{j + 1/2 + i}^{p,i - }} w_{j + 1/2 + i}^p} \right)} \end{array}$ （35）

 ${S^{p,i \pm }} = \pm {\rm{max}}\left( {0,{\rm{min}}\left( { \pm {S^p} - i\frac{{\Delta x}}{{\Delta t}},\frac{{\Delta x}}{{\Delta t}}} \right)} \right)$ （36）

LTS-HLLC格式直接针对方程组构造，以一维Euler方程为例，这里波的个数变为3个，写成叠波形式为

 $\begin{array}{*{20}{l}} {u_j^{n + 1} = u_j^n - \frac{{\Delta t}}{{\Delta x}}\sum\limits_{i = 0}^\infty {\left( {\sum\limits_{p = 1}^3 {S_{j - \frac{1}{2} - i}^{p,i + }} W_{j - \frac{1}{2} - i}^p + } \right.} }\\ {\left. {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \sum\limits_{p = 1}^3 {S_{j + 1/2 + i}^{p,i - }} W_{j + 1/2 + i}^p} \right)} \end{array}$ （37）

 ${S^1} = {S_{\rm{L}}},{S^2} = {S_{\rm{C}}},{S^3} = {S_{\rm{R}}}$ （38）

 $S_{{\rm{L,C,R}}}^{i \pm } = \pm {\rm{max}}\left( {0,{\rm{min}}\left( { \pm {S_{{\rm{L,C,R}}}} - i\frac{{\Delta x}}{{\Delta t}},\frac{{\Delta x}}{{\Delta t}}} \right)} \right)$ （39）

 $\left\{ {\begin{array}{*{20}{l}} {{S_{\rm{L}}} = (1 - \beta )S_{\rm{L}}^{\rm{E}} + \beta S_{\rm{L}}^{{\rm{LxF}}}}\\ {{S_{\rm{R}}} = (1 - \beta )S_{\rm{R}}^{\rm{E}} + \beta S_{\rm{R}}^{{\rm{LxF}}}} \end{array}} \right.$ （40）

4 Godunov型大时间步长格式的高阶精度推广方法

4.1 加权平均状态方法

WAS方法是Toro等[100]提出的用于求解双曲型守恒律的一种二阶精度的数值格式。其仍然按照数据重构、间断分解和通量计算等3步求解数值通量，在数据重构步仍然采用分段常数重构，间断分解仍采用精确或近似Riemann解进行，Toro方法的核心思想在通量计算步，其将原始的Godunov格式的通量计算方法进行了推广，即

 ${{{\mathit{\boldsymbol{\hat f}}}_{j + 1/2}} = \mathit{\boldsymbol{f}}({\mathit{\boldsymbol{U}}_{j + 1/2}})}$ （41）
 ${{\mathit{\boldsymbol{U}}_{j + 1/2}} = \frac{1}{{{t_2} - {t_1}}} \cdot \frac{1}{{{x_2} - {x_1}}}\int_{{t_1}}^{{t_2}} {\int_{{x_1}}^{{x_2}} {{\mathit{\boldsymbol{u}}^*}} } (x,t){\rm{d}}x{\rm{d}}t}$ （42）

 ${t_1} = {t_n},{t_2} = {t_{n + 1}},{x_1} = - \frac{1}{2}\Delta x,{x_2} = \frac{1}{2}\Delta x$ （43）

 ${\mathit{\boldsymbol{U}}_{j + 1/2}} = \frac{1}{{\Delta x}}\int_{ - \frac{1}{2}\Delta x}^{\frac{1}{2}\Delta x} {{\mathit{\boldsymbol{u}}_{j + 1/2}}} \left( {x,\frac{1}{2}\Delta t} \right)){\rm{d}}x$ （44）

 图 10 WAS方法通量计算示意图[11] Fig. 10 Schematic of flux calculation of WAS method[11]
 ${\mathit{\boldsymbol{U}}_{j + 1/2}} = \sum\limits_{k = 1}^{N + 1} {{\beta _k}} W_{j + 1/2}^{(k)}$ （45）

 $\left\{ {\begin{array}{*{20}{l}} {{\beta _k} = \frac{1}{2}({C_k} - {C_{k - 1}})}\\ {{C_k} = \frac{{\Delta t}}{{\Delta x}}{S_k},{C_0} = - 1,{C_{N + 1}} = 1} \end{array}} \right.$ （46）

 ${\mathit{\boldsymbol{U}}_{j + 1/2}} = \frac{1}{2}({W_j} + {W_{j + 1}}) - \frac{1}{2}\sum\limits_{k = 1}^N {{C_k}} \Delta W_{j + 1/2}^{(k)}$ （47）

Toro推导出了具有TVD性质的WAS格式来保证数值计算的稳定性。引入TVD约束条件后，WAS方法的数值通量式(47)可改写为

 $\begin{array}{l} {{\bar W}_{j + 1/2}} = \frac{1}{2}({W_j} + {W_{j + 1}}) - \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{1}{2}\sum\limits_{k = 1}^N {{\rm{ sign }}} ({C_k})\varPhi _{j + 1/2}^{(k)}\Delta W_{j + 1/2}^{(k)} \end{array}$ （48）

4.2 WAS型二阶精度大时间步长Godunov格式

 图 11 WAS方法中单束波的传播示意图[85] Fig. 11 Illustration of single wave propagation for WAS method[85]

 ${U_{(j,1)}} = {U_{(j,2)}} = {U_j}$ （49）

 ${U_{j + 1/2}} = \frac{1}{2}({U_{(j,2)}} + {U_{(j + 1,1)}})$ （50）

 $U_j^{n + 1} = U_j^n - \frac{{\Delta t}}{{\Delta x}}[F({U_{j + 1/2}}) - F({U_{j - 1/2}})]$ （51）

 $\Delta \hat q_{j + 1/2}^k = \varPhi (\gamma _{j + 1/2}^k)\Delta q_{j + 1/2}^k$ （52）

 $\gamma _{j + 1/2}^k = \frac{{\Delta q_{j + 1/2}^k}}{{\Delta q_{J + 1/2}^k}}$ （53）

Toro[11]指出，所选择的变量q只需满足穿过每一个波的跳跃不为0即可。通常可选择密度ρ或内能E。本文作者[85]选择了密度ρJ的取值与迎风方向有关，在此取

 $J = \left\{ {\begin{array}{*{20}{l}} {j - 1}&{C_{j + 1/2}^k \ge 0}\\ {j + 1}&{C_{j + 1/2}^k < 0} \end{array}} \right.$ （54）

5 多维问题推广

5.1 维数分裂

 $\begin{array}{l} \mathit{\boldsymbol{\hat U}}_{i,j,k}^{n + 1} = \mathit{\boldsymbol{\hat U}}_{i,j,k}^n - \lambda (\mathit{\boldsymbol{\hat F}}_{i + 1/2,j,k}^n - \mathit{\boldsymbol{\hat F}}_{i - 1/2,j,k}^n) - \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \lambda (\mathit{\boldsymbol{\hat G}}_{i,j + 1/2,k}^n - \mathit{\boldsymbol{\hat G}}_{i,j - 1/2,k}^n) - \lambda (\mathit{\boldsymbol{\hat H}}_{i,j,k + 1/2}^n - \mathit{\boldsymbol{\hat H}}_{i,j,k - 1/2}^n) \end{array}$ （55）

 $\mathit{\boldsymbol{\hat U}}_{i,j,k}^{n + 1} = (\mathit{\boldsymbol{I}} - \lambda {\delta _i}\mathit{\boldsymbol{\hat F}} - \lambda {\delta _j}\mathit{\boldsymbol{\hat G}} - \lambda {\delta _k}\mathit{\boldsymbol{\hat H}})\mathit{\boldsymbol{\hat U}}_{i,j,k}^n$ （56）

 $\begin{array}{l} ({\lambda ^2}{\delta _i}\mathit{\boldsymbol{\hat F}}{\delta _j}\mathit{\boldsymbol{\hat G}} + {\lambda ^2}{\delta _i}\mathit{\boldsymbol{\hat F}}{\delta _k}\mathit{\boldsymbol{\hat H}} + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\lambda ^2}{\delta _j}\mathit{\boldsymbol{\hat G}}{\delta _k}\mathit{\boldsymbol{\hat H}} - {\lambda ^3}{\delta _i}\mathit{\boldsymbol{\hat F}}{\delta _j}\mathit{\boldsymbol{\hat G}}{\delta _k}\mathit{\boldsymbol{\hat H}})\mathit{\boldsymbol{\hat U}}_{i,j,k}^n \end{array}$ （57）

 $\mathit{\boldsymbol{\hat U}}_{i,j,k}^{n + 1} = \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{I}} - \lambda {\delta _i}\mathit{\boldsymbol{\hat F}}}\\ {\mathit{\boldsymbol{I}} - \lambda {\delta _j}\mathit{\boldsymbol{\hat G}}}\\ {\mathit{\boldsymbol{I}} - \lambda {\delta _k}\mathit{\boldsymbol{\hat H\hat U}}_{i,j,k}^n} \end{array}} \right.$ （58）

5.2 对称维数分裂

 ${\mathit{\boldsymbol{\hat U}}^{n + 1}} = \mathit{\boldsymbol{T}}({\mathit{\boldsymbol{\hat U}}^n})$ （59）

 ${\mathit{\boldsymbol{\hat U}}^{n + 1}} = {\mathit{\boldsymbol{T}}_\xi }{\mathit{\boldsymbol{T}}_\eta }{\mathit{\boldsymbol{T}}_\zeta }({\mathit{\boldsymbol{\hat U}}^n})$ （60）

 ${\mathit{\boldsymbol{\hat U}}^{n + 1}} = \frac{{{\mathit{\boldsymbol{T}}_\xi }{\mathit{\boldsymbol{T}}_\eta }({{\mathit{\boldsymbol{\hat U}}}^n}) + {\mathit{\boldsymbol{T}}_\eta }{\mathit{\boldsymbol{T}}_\xi }({{\mathit{\boldsymbol{\hat U}}}^n})}}{2}$ （61）

 图 12 简单分裂和对称分裂的对比[86] Fig. 12 Comparison of simple split and symmetric split[86]
 ${\mathit{\boldsymbol{T}}_\xi }{\mathit{\boldsymbol{T}}_\eta }({\mathit{\boldsymbol{\hat U}}^n}) - {\mathit{\boldsymbol{T}}_\eta }{\mathit{\boldsymbol{T}}_\xi }({\mathit{\boldsymbol{\hat U}}^n}) = {\tau ^2}({\mathit{\boldsymbol{\hat F}}_\xi }{\mathit{\boldsymbol{\hat G}}_\eta } - {\mathit{\boldsymbol{\hat G}}_\eta }{\mathit{\boldsymbol{\hat F}}_\xi })({\mathit{\boldsymbol{\hat U}}^n})$ （62）

5.3 分片Riemann问题的引入

 $\mathit{\boldsymbol{T}} = \left[ {\begin{array}{*{20}{c}} 1&0&0&0&0\\ 0&{{{\hat \xi }_x}}&{{{\hat \xi }_y}}&{{{\hat \xi }_z}}&0\\ 0&{{{\hat \eta }_x}}&{{{\hat \eta }_y}}&{{{\hat \eta }_z}}&0\\ 0&{{{\hat \zeta }_x}}&{{{\hat \zeta }_y}}&{{{\hat \zeta }_z}}&0\\ 0&0&0&0&1 \end{array}} \right]$ （63）

 $\mathit{\boldsymbol{\bar U}} = \mathit{\boldsymbol{TU}} = \left[ {\begin{array}{*{20}{c}} 1&0&0&0&0\\ 0&{{{\hat \xi }_x}}&{{{\hat \xi }_y}}&{{{\hat \xi }_z}}&0\\ 0&{{{\hat \eta }_x}}&{{{\hat \eta }_y}}&{{{\hat \eta }_z}}&0\\ 0&{{{\hat \zeta }_x}}&{{{\hat \zeta }_y}}&{{{\hat \zeta }_z}}&0\\ 0&0&0&0&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} \rho \\ {\rho u}\\ {\rho v}\\ {\rho w}\\ E \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} \rho \\ {\rho \bar U}\\ {\rho \bar V}\\ {\rho \bar W}\\ E \end{array}} \right]$ （64）

 $\left\{ {\begin{array}{*{20}{l}} {\bar U = {{\hat \xi }_x}u + {{\hat \xi }_y}v + {{\hat \xi }_z}w}\\ {\bar V = {{\hat \eta }_x}u + {{\hat \eta }_y}v + {{\hat \eta }_z}w}\\ {\bar W = {{\hat \zeta }_x}u + {{\hat \zeta }_y}v + {{\hat \zeta }_z}w} \end{array}} \right.$ （65）
 $\left\{ {\begin{array}{*{20}{l}} {{{\hat \xi }_x} = \frac{{{\xi _x}}}{{|\nabla \xi |}},{{\hat \xi }_y} = \frac{{{\xi _y}}}{{|\nabla \xi |}},{{\hat \xi }_z} = \frac{{{\xi _z}}}{{|\nabla \xi |}}}\\ {{{\hat \eta }_x} = \frac{{{\eta _x}}}{{|\nabla \eta |}},{{\hat \eta }_y} = \frac{{{\eta _y}}}{{|\nabla \eta |}},{{\hat \eta }_z} = \frac{{{\eta _z}}}{{|\nabla \eta |}}}\\ {{{\hat \zeta }_x} = \frac{{{\zeta _x}}}{{|\nabla \zeta |}},{{\hat \zeta }_y} = \frac{{{\zeta _y}}}{{|\nabla \zeta |}},{{\hat \zeta }_z} = \frac{{{\zeta _z}}}{{|\nabla \zeta |}}} \end{array}} \right.$ （66）

 ${\left[ {\begin{array}{*{20}{c}} \rho \\ {\rho \bar U}\\ {\rho \bar V}\\ {\rho \bar W}\\ E \end{array}} \right]_t} + {\left[ {\begin{array}{*{20}{c}} {\bar U}\\ {\rho {{\bar U}^2} + p}\\ {\rho \bar U\bar V}\\ {\rho \bar U\bar W}\\ {\bar U(E + p)} \end{array}} \right]_\xi } = {\bf{0}}$ （67）

 $\mathit{\boldsymbol{\bar U}} = \left\{ {\begin{array}{*{20}{l}} {{{\mathit{\boldsymbol{\bar U}}}_{\rm{L}}}}&{\xi < 0}\\ {{{\mathit{\boldsymbol{\bar U}}}_{\rm{R}}}}&{\xi > 0} \end{array}} \right.$ （68）

ηζ方向的分析类同。完成坐标变换后即可直接进行计算，给出每一步的数值解后再采用当地旋转矩阵T的逆矩阵T-1:

 ${\mathit{\boldsymbol{T}}^{ - 1}} = \left[ {\begin{array}{*{20}{c}} 1&0&0&0&0\\ 0&{{{\hat \xi }_x}}&{{{\hat \eta }_x}}&{{{\hat \zeta }_x}}&0\\ 0&{{{\hat \xi }_y}}&{{{\hat \eta }_y}}&{{{\hat \zeta }_y}}&0\\ 0&{{{\hat \xi }_z}}&{{{\hat \eta }_z}}&{{{\hat \zeta }_z}}&0\\ 0&0&0&0&1 \end{array}} \right]$ （69）

 $\mathit{\boldsymbol{U}} = {\mathit{\boldsymbol{T}}^{ - 1}}\mathit{\boldsymbol{\bar U}}$ （70）
5.4 边界条件处理

 图 13 边界附近虚网格示意图[84] Fig. 13 Schematic of ghost cells near boundary[84]
6 Godunov型大时间步长格式的性能分析

6.1 收敛特性

1) 稳定性条件

2) TVD性质

TVD是标量双曲型守恒律方程的一个重要特性，根据Lax-Wendroff定理，只有具有TVD性质的守恒型相容格式才能收敛于弱解，并且可以保证数值解的无非物理振荡。对于3点格式的TVD特性，Harten等开展了充分的研究，取得了很大的成功。但是对于多点格式的TVD特性，长期以来开展的研究并不多，根据第5节的讨论，大时间步长格式就是典型的多点格式。

LeVeque[104]证明了对于标量双曲型守恒律方程，LTS-Godunov格式是具有TVD性质的，因而可以收敛到弱解。Jameson和Lax[105-106]最早研究了广义2K+1点格式的TVD特性，证明了该类多点格式具有TVD格式的充分条件。Lindqvist等[97]扩展了Jameson和Lax的结论，给出了标量双曲型方程的2K+1点格式具有TVD性质的充分必要条件，给出了2K+1点TVD格式的数值黏性系数需满足的上下边界条件，指出LTS-Roe格式和LTS-Lax-Friedrichs格式分别为耗散最小和最大的具有TVD性质的大时间步长格式，据此发展了3.4节的LTS-RoeLxF(β)大时间步长格式。

3) 熵稳定性

Lindqvist等[97]针对标量双曲型方程指出，LTS-Roe格式的数值耗散小于LTS-Godunov格式的耗散，因而可能出现不满足熵条件的情形，并且指出Riemann解中膨胀波的处理是违背熵条件的主要因素。因为在叠波法中，对膨胀波作了多波束近似，故必然存在出现膨胀间断的风险。针对膨胀波处理，Lindqvist等提出了两种方法：对于超声速膨胀波，其建议采用随机选取法来改变时间步长，进而竭力避免膨胀间断的出现；对于跨声速膨胀波，其沿用了Harten[34]提出的熵修正函数。采用这两种策略改进后的格式称为LTS-Roe*，数值实验表明其解的熵稳定性有了明显改善。

6.2 分辨率

1) 误差分析

 $D = \frac{{(|C| + 1 - K)(K - |C|)}}{{2|C|}}$ （71）

 图 14 误差系数随CFL数变化[86] Fig. 14 Error coefficients vs CFL numbers[86]

2) 数值耗散机制分析

Toro[11]提出典型的Godunov型格式的求解过程可以分为3个步骤，即重构、演化和平均。重构主要实现网格节点或单元内平均值分布的高阶近似，故一般不会引入耗散。演化时，耗散大小与所采用的Riemann求解器相关，特别是采用精确Riemann解时，可以认为不引入耗散。平均是为了实现了间断分解后的复杂波系解在下一时刻向网格上的投影，该步是引入耗散的主要环节。文献[85]指出，大时间步长格式的分辨率随CFL数的增大而逐渐提高的原因在于，Godunov型格式求平均步会产生较大的耗散，导致间断面的抹平，对激波和接触间断分辨率不高，而大时间步长格式与普通Godunov型格式相比，由于采取了大的时间步长，在相同的物理计算时间内，减少了求平均的次数，从而降低了耗散，提高了对间断的分辨率。CFL数越大，相同时间内求平均的次数就越少，所以格式的分辨率也就越高。

3) 数值实验验证

6.3 计算效率

 图 15 文献[84]计算效率随CFL数变化 Fig. 15 Computational efficiency vs CFL numbers in Ref.[84]

 图 16 文献[86]计算效率随CFL数变化 Fig. 16 Computational efficiency vs CFL numbers in Ref.[86]

 图 17 文献[115]计算效率随CFL数变化 Fig. 17 Computational efficiency vs CFL numbers in Ref.[115]
7 典型问题应用示例

Godunov型显式大时间步长格式的主要应用对象为以双曲型守恒律方程组为控制方程的问题，主要包含空气动力学领域的Euler方程[82, 84, 86, 97-98]、水力学领域的浅水波方程[111-115]和电磁学领域的Maxwell方程[116]等典型问题。本文主要关注其在空气动力学领域的应用，给出了其在激波管、翼型、机翼等典型一维、二维和三维流动的应用示例。

7.1 一维Sod激波管问题

Sod激波管问题[117]是CFD计算格式验证的经典算例，几乎所有的格式研究都采用其作为验证算例。文献[84, 86, 97-98]均采用所发展的大时间步长格式对该问题进行了计算。文献[84]的计算结果如图 18所示，其通过采用多波束对膨胀波进行近似处理，可以很好地避免非物理解的产生，并且随着CFL数的增大，格式对间断波的分辨率显著提高。文献[86]的计算结果如图 19所示，同样也表明随着CFL数的增大，格式的分辨率提高。文献[97]等的计算结果如图 20所示，其通过熵修正技术可以使得LTS-Roe*格式在CFL数6.0范围内获得较为满意的结果。文献[98]的结果如图 21所示，其发展的LTS-HLL格式基本呈现无数值振荡，但耗散相对较大，LTS-HLLC格式耗散明显减小，呈现略微数值振荡。

 图 18 文献[84]Sod激波管问题计算结果 Fig. 18 Numerical results of Sod shock tube problem in Ref.[84]
 图 19 文献[86]Sod激波管问题计算结果 Fig. 19 Numerical results of Sod shock tube problem in Ref.[86]
 图 20 文献[97]Sod激波管问题计算结果 Fig. 20 Numerical results of Sod shock tube problem in Ref.[97]
 图 21 文献[98]Sod激波管问题计算结果 Fig. 21 Numerical results of Sod shock tube problem in Ref.[98]
7.2 二维翼型绕流

 图 22 NACA0012翼型绕流计算结果(Ma=0.8、α=1.25°)[84] Fig. 22 Numerical results of flow field around NACA0012 airfoil(Ma=0.8, α=1.25°)[84]
7.3 三维机翼绕流

 图 23 文献[84]M6机翼绕流计算结果(Ma=0.839 5、α=3.06°) Fig. 23 Numerical results of flow field around M6 wing(Ma=0.839 5, α=3.06°) in Ref.[84]
 图 24 文献[86]M6机翼绕流计算结果(Ma=0.839 5、α=3.06°) Fig. 24 Numerical results of flow field around M6 wing(Ma=0.839 5, α=3.06°) in Ref.[86]
7.4 高超声速双椭球绕流

 图 25 双椭球高超声速绕流计算结果(Ma=4.0、α=0°)[85] Fig. 25 Numerical results of hypersonic flow field around double ellipsoid (Ma=4.0, α=0°)[85]
8 大时间步长格式研究的发展方向

1) 高阶精度Godunov型大时间步长格式

6.2节的分析表明Godunov型大时间步长格式所取的时间步长越大，则格式的分辨率越高，其在大时间步长下的分辨率甚至高过了常规的高阶精度格式。但是发展高阶精度的大时间步长格式仍是必要的，特别是在提高黏性分辨率方面，有更大的潜力。但是构造高阶格式并不容易，如果采用比分段常数更高精度的高阶重构，如分段线性重构，则会面临所谓的广义Riemann问题的求解，目前对于广义Riemann问题尚未有有效的精确或近似解法，这将引起叠波法的使用困难。为了避免该难题，继续采用分段常数重构相对容易，本文作者[85]曾结合加权平均状态方法和叠波法，构造了具有二阶精度的Godunov型大时间步长格式，但目前CFL数仅能达到2.0。如何克服广义Riemann问题带来的叠波法应用困难，应该是下一步工作要探索的重点。

2) 引入自适应参数的大时间步长格式

3) 方程源项的大时间步长处理

4) 与自适应网格技术结合的真正多维大时间步长格式

http://dx.doi.org/10.7527/S1000-6893.2020.23575

0

#### 文章信息

QIAN Zhansen
Godunov型显式大时间步长格式研究进展
Research progress of Godunov type explicit large time step scheme

Acta Aeronautica et Astronautica Sinica, 2020, 41(7): 023575.
http://dx.doi.org/10.7527/S1000-6893.2020.23575