Godunov型显式大时间步长格式研究进展

钱战森

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

自20世纪七八十年代高性能计算机出现以来，数值计算方法、网格生成技术、并行算法及图像处理等得到了蓬勃发展，数值模拟的应用范围不断拓宽，深度和复杂程度不断加大，解决实际问题的能力也不断提高，使计算流体力学(CFD)逐步成为与理论分析和实验相并列的研究流体力学问题的三大手段之一[1-12]。计算流体力学的应用遍及航空航天、天气预报、海洋与船舶、汽车和列车外型的空气动力学优化、风扇设计、降噪技术、换热器设计等基础工业的各个方面，以及生物工程、土木工程、水利工程、环境保护及风沙治理等工程应用领域[13-16]。特别是在航空航天领域，一方面，计算流体力学的应用可以拓宽和深化流体力学的基础研究，揭示经典方法无法获取的清晰图像，例如飞机在大攻角飞行或机动飞行的情况下，传统的风洞试验手段一般难以取得相关复杂现象的清晰信息；另一方面，相对于风洞试验存在的试验周期长、价格昂贵的缺点，计算流体力学快捷、低成本的优点更为突显。

然而航空航天工程问题数值模拟的计算量是十分巨大的，特别是在大规模工程计算中。例如对复杂飞行器气动特性数值模拟，需要在边界层和分离区等流动结构复杂的区域内密布网格,而在三维计算中，网格尺度每增加一倍的密度，计算量便要增加约16倍。在湍流问题的直接数值模拟中，需要真实雷诺数下的计算，其计算量将和雷诺数的9/4次方成比例。在飞行器气动布局优化过程中，更是需要计算成百上千个样本状态，才能获得最终的优化气动布局形式。根据目前工程使用经验，计算速度仍是制约当前计算流体力学应用的一大瓶颈，在计算机性能和资源相对有限的情况下，建立高效的数值算法以提高计算速度对开展大规模科学与工程计算具有重大的现实意义。

提高计算速度的主要方法包括优化程序设计、并行计算、区域分解、改善数值格式等。优化程序设计主要考虑在空间循环中减少耗费CPU时间的运算过程，但对已开发的格式，通过优化程序设计所能节约的计算量是有限的。随着计算机资源的不断发展，并行计算已成为助力复杂飞行器数值模拟分析的重要手段，但在一定时间内计算机资源总是有限的，模拟规模和数量往往受限。区域分解方法[17]可以根据局部流场的主控尺度特征对不同区域采用不同主控方程进行模拟，也可在一定程度上明显提高数值模拟的效率。改善数值格式，即从数值格式构造方面提高计算速度的方法，是提高数值模拟效率最直接的途径。格式速度每提高一倍，相当于计算机CPU主频提高一倍或并行计算机数目减少一半，这是相当可观的收益。该方法也可以和前述方法耦合使用，获得更大的收益。因此，提高流场解算器数值格式的计算速度对于大规模计算具有特殊的意义。

放宽显式计算格式的时间步长限制，提高计算的CFL(Courant-Friedrichs-Lewy)数，从而提高计算效率，是一种有意义的探索，即发展所谓的显式大时间步长格式。这里CFL数的定义一般可写为CFL=aΔt/Δx，或Δt=CFLΔx/a，其中a为双曲型守恒律方程的Jacobi矩阵的最大特征值(通常表征特征波的最大传播速度)，Δt为离散的时间步长，Δx为离散的空间网格尺度。对于显式格式而言，一般情况下要满足CFL≤1(称为CFL条件)，这导致时间步长受限，从而计算效率受限。这里的显式大时间步长(Large Time Step， LTS)格式是指一种CFL数突破1限制的全离散的、单步的、显式的计算格式，仍满足CFL条件，但CFL条件需要推广，这将在下文中详述。

本文首先简要介绍了显式大时间步长格式的概念、分类和优势，然后重点回顾了针对Godunov型显式大时间步长格式的发展，并给出了该类格式今后的发展方向。本项工作旨在为促进大时间步长格式领域的交流提供参考。

1 显式大时间步长格式简介

1.1 数值格式的发展概述

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

在航空航天计算流体力学中广泛采用时间相关方法[18-20]数值求解Euler或Navier-Stokes(N-S)方程。按照时间推进方式，求解流体力学方程的时间相关格式可以归纳为显式格式和隐式格式；按照离散方式可以分为全离散格式和半离散格式。

全离散格式指的是在离散控制方程时同时考虑时间项和空间项的导数近似，以取得时间和空间方向所具有的设定精度。全离散格式基本都是显式的，一般可通过原偏微分方程以空间导数来表达其时间导数。对于线性常系数方程，更高阶的场变量时间导数离散可容易地以其等价空间导数表达，但对于非线性方程，该变换实现较为复杂，甚至无法显式表达。

半离散格式是指对控制方程进行离散时将时间方向和空间方向离散解耦的方法。该思路可通过分别离散空间与时间导数较容易地导出设定精度的数值格式。通常可先对空间方向作离散，保持时间方向的连续性，导出关于时间变量的常微分方程，称为原控制方程的半离散方程。而导出的常微分方程的离散可以借助已有的标准方法给出。

显式格式在工程实际中，特别是大规模计算中被广泛应用。该类格式的优点是每推进一个时间步长计算量和存储量较少，程序简单，缺点是时间步长受到稳定性条件的限制，导致总体计算步数往往较多。目前已发展了大量可以求解常微分方程组初值问题的显式方法[21-22]，如Euler折线法、Adams外插法、Adams内插法、Milne方法、Hamming算法、Adams三步四阶预估-校正算法及Runge-Kutta方法等。航空航天流体力学计算中使用较多的显式格式为Runge-Kutta方法，其时间精度可达到二或更高阶，时间步长可以适当放宽。但Runge-Kutta方法属于多步法，所有多步法均以增加计算量为代价，每多一层子循环，计算量便增加一倍，不符合高效计算的要求。

隐式格式在模型方程的线性稳定性分析中一般是呈无条件稳定的，在数值求解实际问题时也可取较大的时间步长，目前也得到了较为广泛的应用。航空航天流体力学计算中具代表性的隐式方法有Beam-Warming隐式格式[23]、MacCormack隐式方法[24]、近似因式分解方法[25]、改进的对角化近似因式分解方法[26]、LU-SGS方法[27]以及GMRES方法[28-29]等。这些隐式格式在定常流场计算中得到了广泛应用。然而，一般提高计算效率的近似处理均会导致迭代过程中时间精度的降阶,而低阶精度特别是一阶精度的时间格式，不适用于非定常流的模拟。为了克服这一难题，20世纪90年代发展了模拟非定常问题的双时间步迭代方法[30]。但是目前工程实践经验表明，双时间步迭代方法由于其内迭代常常无法在短时间内收敛或者甚至不收敛，需要通过设定内迭代的次数来降低计算量，该措施必然引入附加误差，甚至可能导致非定常流动的模拟失真。

2) 激波捕捉格式的发展

航空航天领域广泛关注可压缩流动的求解问题，可压缩流动中的激波模拟给数值格式带来了强非线性的重要特征，常用的激波模拟方法包括激波装配法和激波捕捉法。相对于激波装配法，激波捕捉法因其操作简单、实用而成为计算流体力学广泛采用的方法。

Godunov于1959年提出的间断分解法被认为是最早的激波捕捉格式[31]，即通过Riemann问题间断分解的计算来构造差分格式，这类格式统称作Godunov型格式，这个思路直到今天依然是CFD差分格式研究的热点，并且一直处于不断完善和发展之中。近30年来，Godunov型格式的研究和发展主要有两个方向：一是对格式的高阶推广，提高格式的精度以提高对各种间断的分辨率，代表性工作包括Boris和Book[32]、Harten等[33-34]通过对一阶差分格式的数值通量进行修正而获得二阶精度的格式，Jameson等[35]在中心差分格式(JST格式)基础上采用人工黏性法获得二阶精度的格式，van Leer[36-39]及Collela和Woodward[40]通过对初始时刻的解采用分段线性或抛物线重构函数代替原Gudonov格式中的分段常数重构以获得二阶或三阶精度的格式，Harten等[41]提出的基本无震荡(ENO)格式及Liu[42]、Shu[43]等提出并改进的高效WENO格式，通过采用隐式差分可以在较小的模板点上获得更高的插值精度的紧致格式及其改进形式[44-53]；二是对近似Riemann解的研究，用各种近似Riemann求解器来代替原始Godunov格式采用的Riemann问题精确解以减少计量，代表性工作包括Roe[54]、Osher[55-56]、Harten[57]、Toro[58]等分别发展出了各种通量差分分裂形式的近似Riemann解算器；Steger[59]、van Leer[60]、Liou[61-65]、Jameson[66]、张涵信[67-68]等分别发展出了各种矢通量分裂形式的近似Riemann解算器；Lax[69]、MacCormack[70]、Jameson[35]、Tadmor[71-72]、Levy[73-74]等发展的各种以Lax-Friedrichs格式[75]为基础的中心型Riemann解算器。总体来看，对近似Riemann解的研究近年来活跃度不如高阶推广高，其中较有代表性的研究工作是针对多维问题发展的混合型近似Riemann解[76-79]。

当前广泛采用的各种CFD求解方法，都可以看作是Godunov格式的高阶推广或间断分解方法的改进，即都属于Godunov型格式，因而都需要联合采用高阶精度的插值方法和类似Godunov格式的精确或近似的Riemann问题解的数值通量，两者共同组成了可压缩流动的激波捕捉格式。

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

航空航天领域空气动力学数值模拟的主要控制方程为N-S方程，一个完整的N-S方程可以通过算子分裂[18]分解成无黏部分、有黏部分以及源项部分。其中无黏部分即Euler方程，属双曲型；有黏部分属抛物型；源项部分可以是线性或非线性代数函数、常微分或偏微分形式。抛物型方程由于不存在占主导地位的信息传播方向，最佳离散格式为中心差分格式。常微分方程的离散可以借助已有的线性多步法等标准方法处理。

目前CFD领域有关差分格式的研究主要集中在Euler方程相关项，其为典型的双曲型守恒律方程，显式格式和隐式格式均有所发展，增大数值格式的时间步长是提高格式计算速度最有效的途径之一。但是对于显式格式而言，一般情况下时间步长要受到稳定性要求(CFL条件)的限制。隐式格式因其理论上没有时间步长限制，在实际模拟中也可以取较大的时间步长，在工程计算领域得到了更广泛的应用。然而，Euler方程由于其双曲型性质，即解的依赖域是有限的，理论上不宜采用其数值解的依赖域，往往是全场相关的隐式格式求解。

发展显式大时间步长格式，突破显式格式的时间步长限制是另一种提高计算效率的途径。如前文所述，这里的显式大时间步长格式其CFL数可突破1的限制，但仍保持单步、显式、全离散的特征，并且满足稳定性和相容性等基本特性。本文重点关注该类格式的研究进展。

1.3 大时间步长格式的分类

根据目前就双曲型守恒律方程的求解格式开展的研究工作来看，显式大时间步长格式按照构造方法可以分为两类：第1类为几何型构造，第2类为代数型构造。第1类主要基于Godunov型格式构造所得，其中代表性的为LeVeque提出的行波法[80-82]，针对该格式的不足，Guinot[83]提出了时间线插值法，本文作者[84-85]提出了膨胀波多波近似法和近距双波干扰法，Dong和Liu[86]提出了膨胀波单元分解法，分别对其进行了改进，增强了格式的分辨率和鲁棒性；第2类主要采用代数运算来实现，其中具有代表性的为Harten[87]提出的TVD大时间步长格式和董海涛、李椿萱[88]提出的熵条件大时间步长格式。Harten通过算子合并将原TVD格式构造为单步大时间格式，本文作者[89]对其作了改进并将其推广至多维流动的求解。董海涛和李椿萱[88]利用Hamilton-Jacobi方程精确解构造了双曲型守恒律方程的大时间步长格式。第2类方法与第1类方法得到的格式性质差异较大，本文主要关注第1类大时间步长格式，即Godunov型大时间步长格式。

1.4 大时间步长格式的优势

与常规显式格式和隐式格式相比，Godunov型显式大时间步长格式具有以下几个方面的突出优点：

1) 定常问题计算效率高

时间相关法是求解双曲型守恒律方程的典型方法，但时间推进的步长直接决定着流场信息的推进速度。常规显式格式的CFL数不能超过1，大时间步长格式则可以突破这一限制，实际计算结果[82-84]表明，后者计算效率也确实得到了充分的提高。隐式格式的CFL数虽然理论上不受限，但在实际复杂问题应用中往往不能取得过大。

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

显式大时间步长格式可直接应用于非定常问题计算，且因其属于全离散格式，具有时空一致的代数精度。与常规显式格式和以Runge-Kutta方法为代表的多步法相比计算效率明显提高；与隐式格式相比，其不需要双时间步迭代，具有精度更高、鲁棒性更强的优势。

3) 对间断分辨率高

多个研究结果[82-84]表明，在同等代数精度下，随着CFL数的增大，Godunov型大时间步长格式分辨率逐渐提高。特别是对于激波等强非线性间断和接触间断等线性或拟线性间断均具有很高的分辨率，适合于空气动力学的激波干扰、混合层、涡结构、声波等复杂流动现象的求解。

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

大时间步长格式因是显式格式，对计算机内存要求较低。同时因其信息依赖域有限，故而在分区并行中具有较大优势，尤其适合于当前CFD领域普遍采用的几何分区并行计算，一方面不同几何区域界面之间的数据通讯量较小，另一方面收敛效率基本不受几何分区的影响。故大时间步长格式有着非常优越的大规模并行计算可扩展性。与之相比较，隐式格式则在大规模并行计算方面的扩展性稍欠。由于其计算往往是全场相关的，如果要严格按照原始方式推进，则需要在不同几何分区之间进行大量的数据传递，这将占用庞大的时间；若简单地在每一个分区内单独使用诸如LU-SGS等形式的隐式推进格式，则将导致算法的鲁棒性和收敛效率大大降低。

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

长期以来已发展了不少行之有效的基于时间相关法求解定常问题的快速收敛算法，如当地时间步长[90]、残值光顺[91-93]、焓阻尼技术[94]和多重网格技术[30，95]等。大时间步长格式因其显式特点，可与这些加速收敛的算法很好地兼容，大大改善对于定常问题的求解效率，这对于大时间步长格式的工程应用推广具有十分重要的意义。

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

2.1 Godunov型格式

考虑一维双曲型守恒律组的初值问题：

(1)

图1 局部Godunov平均示意图Fig.1 Schematic of local Godunov averaging

(2)

由于函数u(x,tn+1)在区间[xj-1/2,xj]和[xj,xj+1/2]上由不同的间断分解产生，需要进行分段积分。这种方法的实现较烦琐，且时间步长受到CFL≤0.5的严格限制，以避免左右两相邻间断波发生相互干扰。因此，在实际计算中很少采用此方法。

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

(3)

式中：

(4)

显然，以上两种方法的CFL数都要受到限制，即不能超过1，从而给时间推进步长带来了严格的限制。

2.2 Godunov型大时间步长格式

沿着原始Godunov型格式的思想，LeVeque最早提出了一种用于计算双曲型守恒律的大时间步长格式[80-82]。该格式的CFL数可突破1的限制，其核心思想是：沿用分段常数重构函数来代替网格点上的离散值，形成一系列的Riemann问题。在对该系列的Riemann问题作间断分解后，将得到一系列的波(包括激波、膨胀波、接触间断等)。于是可对其进行时间积分求出数值解。当CFL≤0.5时，波系间显然不会发生相互干扰，可直接采用原始Godunov型格式的第1种方法或第2种方法求解。当CFL≤1.0时，采用第2种平均方法也不致出现实质性困难。为了使CFL数进一步提高，LeVeque对波系之间的相互干扰作如下的线性假设：假设两束波相交时互不影响，即两束波的传播速度和强度均保持不变，也无新的波产生，只是简单的相互穿越。

由于大时间步长格式允许CFL>1，因此需面对间断分解后的强波与强波之间的干扰问题：强波相交问题。对于线性方程组，异族强波相互穿越并且穿越后波速与波强不变(直线穿越)，同族强波则变为相互平行的直线，故LeVeque的上述假设在常系数线性双曲型守恒律是精确成立的；对于非线性方程，异族强波相互穿越并且穿越后波速与波强会有改变，同族强波一般不易相交，如相交则合并成一个强波。可见对非线性方程，LeVeque假设会带来一定误差，但是实际计算结果表明，其仍可获得较好的计算结果。

(5)

则式(5)的系列分片常数Riemann问题

(6)

的解可表达为

(7)

首先考虑单束波的简单问题，即由标量双曲型守恒律主控的问题，如图2所示[84]。设在点(xj+1/2,tn)处有Δu=uj-uj+1的间断，sa为该间断的传播速度。当sa>0时，对于x∈[xj+1/2,xj+1/2+saΔt]，un(x,tn+1)的值相对于un(x,tn)存在Δu的变化，即un+1(x,tn+1)=un(x,tn)+Δu。若该束波的传播已越过[xj+1/2,xj+3/2]的整个区间，则对区间[xj+1/2,xj+3/2]内的每一点x，un(x,tn+1)在时间段[tn,tn+1]均增加了Δu，于是整个区间的平均值也增加了Δu。若在时间段[tn,tn+1]内波束未穿过整个区间，仅到达区间[xj+1/2,xj+3/2]内某一点ε(xj+1/2<ε

图2 单波束在网格上的传播[84]Fig.2 A single wave propagation on grid[84]

(8)

同样，对于sa<0的情况也可作类似的处理。

采用显式格式数值求解该问题，由稳定性条件(CFL条件)知微分方程数值解的依赖域必须包含精确解的依赖域，因此CFL数取K时，至少需2K+1个模板单元，此时推进时间步长可以达到Δt=KΔx/sa，即可能有2K个波束(左右各K个)对目标单元有贡献，图3仅给出了特征波速为正的波束示意图[86]。根据叠波法思路，CFL=K的大时间步长格式可写成

图3 多波束在网格上的传播[86]Fig.3 Multi-wave propagation on grid[86]

(9)

(10)

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

对于守恒律系统，与单方程不同的是间断分解比较复杂，每一间断可以分解出多个波，如图1所示。可仍假设不同波之间相互自由穿越，互不干扰。这个假设对非线性方程也属于线化近似，在一定CFL数范围内该近似足够准确，当CFL数非常大时，它会导致解的扭曲，此时可以通过加密网格来弥补(可以证明[80]，当网格加密时，该方法对所有CFL数都是收敛的)。图4给出了Euler方程(每个间断包含三波束)的多波束叠加示意图[86]，这时波束虽然变得复杂，但在标量情形下发展的扫描方法仍可直接推广应用。

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

在方程组的情形下，式(10)可以推广为

(11)

(12)

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

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

3.1 时间线插值技术

时间线(Time-line)插值技术是Guinot[83]提出的一种改进方法，旨在减小数值通量计算的模板点数，从而进一步降低计算量，其主要思路如下。

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

(13)

式中：

(14)

变量u可以写为如下的线性组合：

(15)

(16)

按照特征值的正负，可将式(16)改写为

(17)

考虑界面上的等价Riemann问题，有

(18)

在特征变量空间内，按照特征线法，有

(19)

考虑λ(m)为正值的情形，在CFL数不超过1的条件下，根据特征线法可将时间积分转换为空间积分：

(20)

式中：xA点是特征线的起始点，在CFL≤1条件约束下，一般落于xj或xj+1单元。

在CFL数超过1的条件下，xA点将落在xj单元左侧或xj+1单元右侧。因为界面xj-1/2或xj+3/2处可能有激波存在，特征线将是不可逆的，基于逆向追踪法的式(20)将无法直接使用。此时，界面通量可以改写成

(21)

式中：τ表示从界面xj-1/2发出的特征线与界面xj+1/2的交点，如图5所示。仍考虑λ(m)为正值的情形，进一步可有

图5 时间线插值法示意图[83]Fig.5 Schematic of time-line interpolation method

(22)

(23)

式中：τ′表示从界面xj+1/2逆向发出的特征线与界面xj-1/2的交点。式(23)的表达相当于把特征线不可跨域的激波区采用相邻的界面值来代替。

对于界面xj-1/2的变量采用线性重构，则有

(24)

式中：中点值和斜率的表达式分别为

(25)

(26)

其中：s为线性重构的斜率，可以对s引入限制器以防数值解振荡，具体可参见文献[83]。

该格式需要从左向右依次推进，以得到所有波速为正的界面通量贡献，从右向左依次推进，以得到所有波速为负的界面通量贡献。所以，数值通量表现出很强的边界条件依赖性，需要妥善处理边界条件，详情可见文献[83]。本质上来讲，因为无法实现真正意义上的局部推进，该格式仅是一种半显式格式，可以看作是对隐式格式的一种简化。

3.2 膨胀波的处理技术

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

针对Euler 方程的特点，本文作者与Lee[84]提出采用若干个膨胀间断波来近似间断分解中出现的膨胀波，可称之为LTS-Godunov-MW格式。如图6所示，各个膨胀间断波之间的状态由膨胀波波头和波尾的状态通过线性插值得到。这里的膨胀间断波严格意义上也是违反熵条件的，但是因为与压力相比，熵变随着近似波数目的增多衰减要快得多，所以使得最终数值解出现违反熵条件的概率急剧下降。尝试了分别采取一束膨胀间断波(其本质就是一束膨胀激波，简称单波近似)、两束膨胀间断波(简称双波近似)和三束膨胀间断波(简称三波近似)的处理。理论上可以采用更多波束(多于三波)来代替膨胀波，且采用越多波越逼近于原膨胀波，但也相应加大了计算量。数值实验表明，在绝大数情况下，两波近似已可提供足够的精度，既能避免非物理解的出现，又不致过大地增加计算量。

图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 波系干扰的处理技术

对于非线性双曲型守恒律组，如一维Euler方程，在CFL>0.5时，如图8所示[84]，两个相邻间断发出的波即可能发生干扰。对于诸如Euler方程的非线性方程，异族强波相互穿越并且穿越后波速与波强均会有改变，同族强波不易相交，如相交则合并成为一个强波。在实际问题中，非线性波系的干扰非常复杂。相邻的两个间断分解后，各发出3个波——左传波(可以是左向激波或左向膨胀波)、接触间断、右传波(可以是右向激波或右向膨胀波)，其中激波和膨胀波为非线性波，接触间断为线性波。

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

从图8可以看到，在给定的时间步Δt内仅两个相邻间断就可能出现多个波对干扰问题，对于N个间断的情形，将出现更多个波对的相互干扰。文献[96]对两族波之间的相互作用，包括激波、膨胀波和接触间断的两两相互作用进行了详细论述，给出了解析解。但是，在实际应用中，考虑完整的波系干扰(包括不相邻的两个或多个间断)，在编程和计算等实际操作上是难以实现的。为简化计算，文献[84]提出仅考虑相邻的两个间断发出的相邻波——左间断发出的右传波(图8中波a)和右间断发出的左传波(图8中波b)——之间的碰撞问题。经以上简化后，只需处理两个异族波的对撞问题，如图9所示。其对撞存在4种可能性：两个激波的对撞、左膨胀波和右激波的对撞、左激波和右膨胀波的对撞、两个膨胀波的对撞，可采用文献[96]给出的解析解来实现。波a和波b碰撞后必然产生两个新的非线性波和一个接触间断，即两波之间的状态由单一状态M变为两个由接触间断所间隔的状态M*(若波a和波b均为膨胀波时则不存在接触间断)。数值计算结果表明，通过采用相邻间断干扰精确处理的方法，使得格式的鲁棒性大大提高，稳定计算的CFL数范围得到有效扩展。

图9 两波碰撞示意图[84]Fig.9 Schematic of interaction of two discontinuities[84]

3.4 近似Riemann求解器的应用

在大时间步长格式的研究中，文献[80, 84, 86]的工作均采用了Riemann问题的精确解来进行间断分解。近期Lindqvist[97]、Prebeg[98]等尝试了采用近似Riemann求解器来进行间断分解。

1) LTS-Roe格式

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

(27)

基于通量差分分裂形式，标量双曲型方程的大时间步长格式可写为

(28)

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

(29a)

(29b)

式中:i∈{1,2,…,K-1}。

(30)

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

(31a)

(31b)

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

(32a)

(32b)

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

(33)

该格式属于LTS-Roe格式和LTS-Lax-Friedrichs格式的线性组合，通过LTS-Lax-Friedrichs格式引入了更多的数值黏性，消弱了数值解在激波附近的非物理振荡，并使得数值解朝熵解靠拢。其在数值算例中探索了参数β对计算结果的影响。

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格式写成数值黏性形式为

(34a)

(34b)

式中：SL和SR分别为HLL格式中近似Riemann问题的解中两道波的波速。

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

(35)

式中：w为波强，即分解后所得的波两侧状态的差量。其中S1=SL和S2=SR，并且有

(36)

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

(37)

式中：

S1=SL,S2=SC,S3=SR

(38)

其中：SC为HLLC格式中接触间断的传播速度，进一步有

(39)

容易理解，LTS-HLL格式耗散较大，LTS-HLLC格式耗散较小，分辨率有明显提高。但是从数值实验也发现LTS-HLLC格式在间断附近容易滋生非物理的数值振荡。为改善该类格式的数值特性，参考HLLE格式[99]的波速选取，Prebeg建立了LTS-HLLE格式，为进一步改进格式，Prebeg等[98]提出了LTS-HLLE格式与LTS-Lax-Friedrichs格式的混合格式LTS-HLLELxF(β)格式：

(40)

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

自LeVeque[80-81]提出大时间步长格式以来，所发展的绝大多数格式均只具有一阶代数精度。虽然在同等代数精度下，随着CFL数的增大，Godunov型大时间步长格式分辨率逐渐提高，甚至超过了常规二价代数精度的格式，但是构造更高阶的格式对于提高格式分辨率还是大有裨益的。本文作者[85]首次尝试了构造二阶精度大时间步长的工作，其综合采用了Billett和Toro[100]提出的加权平均状态(Weighted Average State，WAS)方法和LeVeque的叠波法[80]，建立了具有二阶精度的Godunov型大时间步长格式，但目前CFL数仅能达到2.0。

4.1 加权平均状态方法

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

(41)

(42)

式中：u*(x,t)为tn时刻间断分解后的数值解。式(42)是一种广义的数值通量计算方法，本节的讨论仅局限于

(43)

的局域。与原始Godunov型格式仅有时间方向的积分不同，改进的数值通量的计算引入了空间方向的积分。为获得二阶精度格式，这里对时间方向采用中点积分方式，则有

(44)

如图10所示[11]，对区域[-Δx/2,Δx/2]进行积分，可得

图10 WAS方法通量计算示意图[11]Fig.10 Schematic of flux calculation of WAS method[11]

(45)

(46)

式中：Ck为波速为Sk的第k个波对应的CFL数。式(45)也可以改写为

(47)

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

(48)

式中：Φ为限制器函数，具体形式见文献[100]。

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

采用WAS方法的起点是引入状态变量的加权平均式(42)，然后采用式(41)给出数值通量。对于Euler方程，式(42)中的变量状态可以选择守恒变量、原始变量或其他状态参数，均不影响格式的守恒性。为了同原有的一阶Godunov型大时间步长格式相容，文献[85]仍采用守恒变量，具体构造思路如下。

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

步骤1设置

U(j,1)=U(j,2)=Uj

(49)

步骤3由间断xj+1/2处两侧的两个半区间值(U(j,1)和U(j,2))的平均值给出计算数值通量Fj+1/2的Uj+1/2，即

(50)

步骤4采用守恒格式更新变量值，即

(51)

(52)

(53)

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

(54)

由于在时间方向采用了中点积分，故而该格式的CFL数理论上仅能达到2.0，可以考虑构造使用更多点的积分进一步扩展CFL数，但算法的复杂度将大幅增加。

5 多维问题推广

目前提出的大时间步长格式都是针对一维问题的，向多维情形的推广主要采用维数分裂法[101]。维数分裂法虽然会引入一定误差，但是在一定精度范围内还是可以满足要求。

5.1 维数分裂

考虑多维Euler方程[85]，其离散形式可写成

(55)

式中：λ=Δt/Δξ=Δt/Δη=Δt/Δζ，这里ξ、η、ζ为曲线贴体坐标系下的广义坐标。记算子δi(·)=δi+1/2(·)-δi-1/2(·)，则式(55)可以写成算子形式

(56)

若给式(56)右端添加一个二阶小量:

(57)

则可得到

(58)

以上形式的全离散型显式差分形式可以在i、j、k3个方向上分别进行准一维的计算，程序设计也大大简化。

5.2 对称维数分裂

式(58)可以改写成

(59)

式中：T为推进算子。将T分解为3个方向的子推进算子，则有

(60)

常用的方法是如式(60)所示的简单的维数分裂，但对于某些情形该分裂方法可能会导致一定的数值振荡，尤其在大CFL数时问题可能更突出。Dong和Liu[86]发现这是由于非对称的维数分裂造成的，他们建议使用所谓的对称算子分裂。以二维情形为例：

(61)

图12[86]给出了算子分裂的示意，使用3个间断的解为例来分析，3个间断形成一个三波点，可用于考察简单分裂和对称分裂对于斜间断数值解的影响。从图中可以看出，对于斜间断解，简单分裂与分裂的顺序有关，一个在左上，一个在右下，两者间存在一个如下的非对称小量差异：

图12 简单分裂和对称分裂的对比[86]Fig.12 Comparison of simple split and symmetric split[86]

(62)

式中：Tξ=I-τTξ；Tη=I-τTη。即便是这个小量的存在，也会改变斜间断的形状，使之发生畸变。在多维情形，这个畸变可能会导致间断附近出现非物理的振荡。使用对称分裂方法可以去除该非对称差异，可以抑制斜间断的变形及其附近的数值解振荡。

5.3 分片Riemann问题的引入

一维情形下，Euler方程含有3个变量，但对于三维情形，每一方向均需要处理5个变量的问题。文献[85]参考了文献[102]的方法，采用当地旋转矩阵 将其变换至局部笛卡尔坐标系。该变换矩阵的形式为

(63)

(64)

式中：

(65)

(66)

经过坐标变换，物理域中位于(x,y,z)点处ξ、η或ζ方向的间断分解问题可转换为局部笛卡尔坐标系下的一维Riemann问题。该问题和原始笛卡尔坐标系下的一维问题具有相同的形式，但多出了两个速度分量，Toro[11]称之为分片Riemann问题。以沿ξ方向的三维分片Riemann问题为例，其变换后的形式可写为

(67)

其初值记为

(68)

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

(69)

(70)

5.4 边界条件处理

在绝大多数算法研究的论文中仅考虑了一维激波管算例，在该算例中，边界条件均采取一阶外推给定，处理较为容易。在数值模拟多维问题时，边界条件的设置则相对复杂。本文作者[84-85]通过在计算域边界外添加虚网格的方式简化边界条件的设置，时间步长越大则所需的虚网格数量越多，如图13所示。多维计算中，入流远场边界虚网格上的值采用入流无反射条件给定；出流边界采用出流无反射外推边界条件；壁面边界虚网格上的值由镜像反射条件[103]给定。对于Euler方程，其壁面切向流动不受约束，可使用镜像反射法给出边界外的虚流场。对于N-S方程，壁面满足无滑移条件，虚网格上速度可由反对称镜像法给定。

图13 边界附近虚网格示意图[84]Fig.13 Schematic of ghost cells near boundary[84]

6 Godunov型大时间步长格式的性能分析

自LeVeque提出大时间步长格式以来，有多个工作研究了该类格式的性能，但绝大多数工作集中在针对标量双曲型守恒律的稳定性条件、TVD特性、熵稳定性、分辨率等性能的探讨，对于方程组的情形，目前开展的理论分析仍很少见。当然，计算效率作为大时间步长格式的主要研究出发点，也得到了较为广泛的关注。本节简要介绍有关的代表性工作。

6.1 收敛特性

1) 稳定性条件

直观上看，该类格式似乎违背了CFL条件，无法实现稳定计算，其实即使在CFL>1时，其依然满足稳定性条件。CFL条件指出数值解的依赖域不得小于物理条件下的依赖域，其实该条件并未限定具体的CFL数。根据CFL条件，当格式的模版点为2K+1时，其CFL数的上限可达到K。文献[85]指出Godunov型大时间步长格式正是通过增加CFL数来自动增加模版点的数目以保证满足CFL条件。该格式采用波的传播速度来完成模版点选择的自适应性，即CFL越大，每一时间步长中涉及的模版点数就越多。

为何普通多点格式如二阶TVD格式、MUSCL格式、WENO格式等在时间方向显式离散时，其CFL数仍不能超过1。这是因为在二阶TVD、MUSCL、WENO等格式中，除了3 个主要模版点之外，其他的模版点仅存在于系数之中(主要体现在限制器中)，处于次要地位，因而不能获得更大的CFL数。但该处理可带来光滑区格式精度的提高。

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) 熵稳定性

在标量双曲型守恒律方程情形下，具有TVD性质的守恒型格式收敛于弱解，但是弱解不一定是熵解，只有具有熵稳定性的格式才能收敛至熵解。对于常规格式的熵稳定性，目前已有较多的结论，但是对于大时间步长格式的熵稳定性的研究尚很不充分。LeVeque[104]曾推测，只要间断分解过程中使用的Riemann求解器可提供熵解，则采用叠波法的LTS-Godunov格式就可以收敛到熵解，但时至今日，该推测依然没有得到严格证明。Wang和Warnecke[107-108]曾证明：在CFL≤2的条件下，如果通量函数是单调的，则LTS-Godunov格式是熵稳定的；如果初始值是单调的，则LTS-Godunov格式对于任意CFL数都是熵稳定的。Wang等[109]证明了在某些特定的初值条件下，LTS-Godunov格式对于任意CFL数都是熵稳定的。Tang和Warnecke[110]研究了LTS-Lax-Friedrichs格式，证明了因其为单调格式，故而是熵稳定的。

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

6.2 分辨率

多个研究均表明Godunov型大时间步长格式的分辨率与时间步长相关，且时间步长越大，则格式的分辨率越高，即引入的耗散越小。本节从误差分析、数值耗散机制分析和数值实验验证3个 角度简述对该特性的研究。

1) 误差分析

文献[86]给出标量双曲型守恒律方程的大时间步长格式的截断误差系数为

(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) 数值实验验证

自大时间步长格式提出以来，多个研究者均对其进行了数值算例验证。文献[80，82]针对标量方程和等熵Euler方程，文献[83，111-115]针对浅水波方程，文献[84，86，97]针对Euler方程均进行了大量数值实验，计算结果均表明Godunov型大时间步长格式具有随着CFL数增大而分辨率提高的性质。这与上述的误差分析和数值耗散机制分析的结论是一致的。

6.3 计算效率

在针对大时间步长格式的研究中，绝大多数工作仅针对格式本身的特性，并未直接给出计算CPU时间的直接对比。文献[84,86,115]等各自给出了LTS-Godunov格式的CPU时间，评估了所发展的各种大时间步长格式的计算效率，给出了其在典型模型问题求解中的加速特性。

文献[84]给出了采用双波近似的LTS-Godunov-MW格式对一维Sod激波管和二维NACA0012翼型的流场计算的效率对比,如图15所示[84]。在CFL数5.0以内，对一维Sod激波管问题格式加速比随CFL数增加而提高，在CFL数5.0时能达到23%；对二维NACA0012翼型问题加速特性有所下降，在CFL数5.0时加速比能达到42%左右。但该格式效率还存在优化的空间。

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

文献[86]给出了改进后的LTS-WA格式对一维Burges方程、Sod激波管问题和二维Ringleb流动的流场计算的效率对比,如图16所示[86]。对于一维Burges方程，LTS-WA格式在CFL数9条件下的加速比约为53%；对于一维Sod激波管问题，LTS-WA格式在CFL数9条件下的加速比能达到21%；对于二维Ringleb流动，LTS-WA格式在CFL数8条件下的加速比能达到25%。

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

文献[115]给出了LTS-Roe格式在一维浅水波方程求解中的效率对比，并研究了膨胀波多波束近似和CFL数之间的关系，多波束处理会增加计算时间，但总体上来看增大CFL数还是可明显提高计算效率，如图17所示[115]。在双波近似条件下，以原始Roe格式在CFL数0.8条件下的计算时间作为参考，LTS-Roe格式在CFL数4条件下的加速比能达到45%。

图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 二维翼型绕流

目前大部分的大时间步长格式研究工作主要集中在一维格式的构造上，对二维问题进行计算的工作不多。Dong和Liu[86]曾采用LTS-WA格式计算了二维Riemann问题，但仍采用了等距网格。文献[84]首次将LTS-Godunov格式用于二维NACA0012翼型绕流计算。来流条件分别为马赫数Ma=0.8、攻角α=1.25°和马赫数Ma=1.2、攻角α=7.0°。文中两个绕流情形均给出了CFL数从1.0～5.0的计算结果，可以发现随着CFL数的增大，LTS-Godunov格式的分辨率逐渐提高，该趋势尤其明显地体现在对下翼面强度较弱的激波的分辨能力上，如图22所示。

图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 三维机翼绕流

目前仅文献[84]和Dong[86]等给出了采用LTS-Godunov格式计算M6机翼三维绕流的结果，计算的来流马赫数Ma=0.839 5、攻角α=3.06°。此算例为超临界状态，流场结构较复杂，上翼面共有3道激波：前一道为弱激波，激波后仍为超声速流动；后一道激波为强激波，波后为亚声速流动；此两道激波在翼梢附近相交，汇合成λ型第3道最强激波。下翼面均为亚声速流动。图23 给出了文献[84]的部分计算结果，显示了CFL数1.0～5.0时的上壁面压力云图。从图中可以看出，大时间步长Godunov型格式能正确捕捉到上翼面由于两道激波相交而形成的λ激波，且和一维及二维情形相似，随着CFL数的提高，该格式的分辨率也逐步提高，主要表现在λ激波的捕捉上。文献[84]还给出了相应CFL数下的压力分布及其与风洞试验结果的对比，验证了该类格式的计算可靠性。文献[86]给出了CFL数最大到8.0 的计算结果，得到的流场分布如图24所示，结论与文献[84]基本一致，且文献[86]针对80% 展长处的压力分布捕捉精度不足的问题，通过对网格的改进得到了更优化的计算结果。

图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 高超声速双椭球绕流

本文作者[85]还给出了LTS-Godunov格式对高超声速双椭球绕流的计算结果，来流马赫数Ma为4.0、攻角α为0°，计算结果如图25所示。其指出大时间步长Godunov格式在不同CFL数下的计算结果和Haten-Yee TVD格式[118]所得到的结果相符，都能准确地捕捉到椭球头部的脱体弓形激波和两椭球镶嵌部位附近的二次激波，图中给出了CFL数为1.0，3.0和5.0时双椭球物面及流场对称面上的等压线分布和对称面上的上下壁面压力分布。

图25 双椭球高超声速绕流计算结果(Ma=4.0、α=0°)[85]Fig.25 Numerical results of hypersonic flow field around double ellipsoid (Ma=4.0,α=0°)[85]

8 大时间步长格式研究的发展方向

自LeVeque提出大时间步长格式的概念以来，经过30多年的发展，Godunov型大时间步长格式引起了多个研究者的兴趣。特别是在近年来随着超级计算机系统的迅速发展，该类格式因其优异的并行特性，得到了更广泛的关注。经过逐渐探索和完善，愈来愈接近推向工程应用计算的阶段，但是仍有一些问题有待解决。根据作者观点，以下的发展方向值得进一步探索。

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

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

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

目前基于叠波法所发展的大时间步长格式具有无自由参数的优点，但是从数值算例还是可以发现该格式在大CFL数下，在激波波后仍有微幅数值振荡，这可能是数值耗散不足造成的。Lindqvist等[97]结合了LTS-Roe和LTS-LxF两种格式的优劣，提出了混合型的LTS-RoeLxF(β)格式；Prebeg等[98]提出了LTS-HLLE格式与LTS-Lax-Friedrichs格式的混合型LTS-HLLELxF(β)格式。其目的都是通过LTS-LxF引入更多的耗散，但是其引入的混合参数β为固定值，对于复杂流动计算难以得到表现优异的全场统一的混合参数。故而，发展自适应形式的混合参数，使其根据流场结构的特点自适应地调节，是一种优化的策略。对于叠波法而言，该思路也是可以尝试的，在叠波过程引入自适应调节参数，也可方便调节格式的耗散特性。

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

在航空飞行器绕流的数值模拟中，黏性是不可忽略的重要因素。当前飞行器空气动力学研究普遍采用RANS湍流模型来模拟黏性湍流流动，湍流模型中一般存在描述湍流生成和耗散的源项。在考虑化学非平衡的高超声速流动或燃烧效应的流动中也存在描述化学反应效应的源项。因此，开展带源项的双曲型守恒律的大时间步长格式研究具有重要的实际工程意义。Murillo[111]、Hernandez[113]等在浅水波方程的求解中发展了处理渠底坡度和摩擦损失等带来的源项处理方法，可为空气动力学领域的控制方程求解所借鉴。但空气动力学领域控制方程的源项一般非线性更强，故处理难度更大，是今后需要重点探讨的方向之一。

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

自适应网格技术是空气动力学数值模拟的重要手段，可随着流场结构的变化动态调整网格的分布，是高精度求解复杂流场问题的高效手段。目前的大时间步长格式均基于一维情形构造，然后借助维数分裂法来实现多维问题推广。但实际上Godunov型大时间步长格式的叠波法具备推广至多维情形的潜力，并且叠波法具有先天追踪波系信号的特点，在采用自适应网格技术方面具有独特优势。如能将叠波法和自适应网格技术相结合，则可能构造出显式的自适应网格大时间步长格式，具备更高效求解空气动力学问题的能力，也是值得探索的方向。

致 谢

感谢清华大学陈海昕教授，本文是在他的鼓励下完成的。