← 回到首页 Skip to main content

3.1 空间离散

第 3 章 控制方程的求解原理

上一章给出了完整的 Navier–Stokes/Euler 方程组,补充了完全气体的热力学关系,并为化学反应气体定义了附加输运方程。现在,我们可以求解整个控制方程组,以确定各流动变量。可以想见,求解方法多种多样。若暂不考虑仅适用于简化流动问题的解析方法,几乎所有求解策略都遵循同一基本过程。

首先,将需要计算流动的空间,即物理空间,划分为大量几何单元,称为网格单元。这个过程称为网格生成;有些作者使用 mesh 一词,含义相同。也可以将其理解为:先在物理空间中布置网格点,也称节点或顶点,再用直线,即网格线,将这些点连接起来。二维(2-D)网格通常由三角形和/或四边形组成;三维(3-D)网格通常由四面体、六面体、棱柱或棱锥组成。网格生成工具必须满足的最重要要求,是网格单元之间既不存在空隙,也不相互重叠。此外,网格应当光顺,即网格单元的体积或拉伸比不能发生突变,单元形状应尽可能规则。对于四边形或六面体网格,网格线还应避免出现剧烈折转,否则数值误差将显著增大。

一种做法是使网格紧密贴合物理空间的边界,这种网格称为贴体网格(图 3.1a)。其主要优点,是能够在边界处十分精确地解析流动。

Image
图 3.1 固体附近的二维贴体网格(a)、采用浸入边界的笛卡尔网格(b),以及采用多边形边界单元的笛卡尔网格(c)。图(c)中以粗线标出边界单元。图内 Solid body:固体。

这种边界解析能力对于固体表面附近的剪切层至关重要。代价是网格生成工具会变得十分复杂,尤其是在处理实际工程几何形状时。另一种做法,是使用所谓的笛卡尔网格 $ [1-8] $,其单元边与笛卡尔坐标轴平行,因此生成起来非常容易。相较于贴体网格,这种网格上式(2.19)的通量计算也简单得多。不过,从图 3.1b、c 可以看出,要以通用且精确的方式处理边界,并不容易 [3]。因此,通常优先采用贴体网格或组合方法(见第 11.2.4 节),尤其是在模拟对象几何形状往往十分复杂的工业环境中。

如今,求解 Euler 和 Navier–Stokes 方程的数值方法,绝大多数都分别进行空间离散和时间离散,即采用所谓的线方法(method of lines)$ [9] $。根据所选算法,网格可用于构造控制体并计算通量积分,也可用于近似流动物理量的空间导数。随后,从一个已知初始解出发,使用适当方法对所得时间相关方程进行时间推进。如果流动变量不随时间变化,还可以通过迭代过程求得控制方程的定常解。

按照前面推导控制方程(2.19)的方式,连续性方程(2.3)中包含密度对时间的导数。密度作为独立变量,用于计算压力(式(2.29)),因此,在动量方程中,密度与压力的时间演化相互耦合。基于这一原因,对控制方程(2.19)进行离散的求解方法称为基于密度的格式。这一表述的困难在于:对于不可压缩流体,连续性方程中的密度时间导数为零,压力便不再由某个独立变量驱动。另一个困难是,随着 Mach 数降低,声波传播速度与对流波传播速度之间的差异越来越大,使控制方程的刚性不断增强,因而更难求解 $ [10] $。为处理这一问题,主要发展了三种方法。第一种是求解从动量方程推导出的压力 Poisson 方程 $ [11-17] $,这类方法称为基于压力的格式。第二种称为人工可压缩性方法,其基本思想是用压力的时间导数替代连续性方程中密度的时间导数 $ [18, 19] $,从而直接耦合速度场和压力场。第三种也是最通用的方法,是对控制方程进行预处理 $ [20-31] $。这一方法允许对极低 Mach 数和高 Mach 数流动使用同一套数值格式。第 9.5 节将作更详细的讨论。

下面将介绍不同求解方法的基本原理,涉及控制方程在空间和时间上的数值近似、湍流建模,以及边界处理。

3.1 空间离散

首先考察求解过程的第一步:Navier–Stokes 方程的空间离散,也就是对流通量、黏性通量以及源项的数值近似。过去已为此提出许多方法,相关发展仍在继续。为加以归类,可先将空间离散格式分为三大类:有限差分格式、有限体积格式和有限元格式。这些方法都依赖某种网格来离散控制方程(2.19)。基本上,网格可分为两种类型:

  • 结构网格(图 3.2):每个网格点,即顶点或节点,都由索引 i、j、k 唯一确定,相应的笛卡尔坐标分别为 $ x_{i,j,k} $、$ y_{i,j,k} $ 和 $ z_{i,j,k} $。二维网格单元为四边形,三维为六面体。如果网格像图 3.1a 那样贴合物体边界,也可称为曲线网格。
Image
图 3.2 结构贴体网格方法的二维示例:(a)为物理空间;(b)为计算空间;$ \xi $ 和 $ \eta $ 表示曲线坐标系。图内 Solid body:固体。
Image
图 3.3 二维非结构混合网格;数字标示各个单元。图内 Solid body:固体。
  • 非结构网格(图 3.3):网格单元与网格点没有特定的排列次序,因此,不能仅凭索引直接确定相邻单元或相邻网格点,例如单元 6 可能与单元 119 相邻。过去,二维单元为三角形,三维单元为四面体。如今,为充分解析边界层,非结构网格在二维中通常由四边形和三角形混合组成,在三维中则由六面体、四面体、棱柱和棱锥混合组成。这类网格称为混合网格。

结构网格的主要优点,来自索引 i、j、k 构成线性地址空间这一性质。该空间也称为计算空间,因为它直接对应流动变量在计算机内存中的存储方式。只需对相应索引加减一个整数,就能快速、方便地访问某网格节点的邻点,例如 $ (i + 1) $、$ (k - 3) $ 等(见图 3.2)。这一性质显著简化了梯度、通量的计算以及边界条件的处理。由于通量 Jacobian 矩阵具有规则的带状结构,隐式格式的实现也因此得到简化。但结构网格也存在不足:对于复杂几何形状,网格生成十分困难。如图 3.4 所示,一种办法是将物理空间划分为若干拓扑结构较简单的部分,即网格块,以便分别生成网格。这就是多块方法 $ [32–36] $。当然,流动求解器会随之变得更加复杂,因为需要专门的处理逻辑,在网格块之间交换物理量或通量。如果界面两侧的网格点可以彼此独立地布置,即不要求网格线在块边界处相接(如图 3.4 中 C 或 F 内部所示),还可以获得更大的灵活性。只存在于块界面一侧的网格点,称为悬挂节点。这种方法的优点显而易见:各网格块的网格线数量

Image
图 3.4 块间界面匹配/非匹配的结构多块网格;粗线表示网格块边界。

可以根据需要分别选择。灵活性提高的代价,是对悬挂节点进行守恒处理所需的附加开销增加 $ [37, 38] $。多块方法还为通过区域分解,在并行计算机上实现流动求解器提供了有益途径。然而,对于复杂外形,网格生成仍需较长时间,往往达到数周甚至数月。

另一种与分块结构网格有关的方法,是 Chimera 技术 $ [39–46] $。其基本思想是:先分别为计算域中的各个几何组成部分生成网格,再将这些网格组合起来,使其在交接处相互重叠。图 3.5 给出了一个简单外形的示意。关键操作是,在重叠区域准确地传递不同网格之间的物理量。因此,重叠区域的范围需按所要求的

Image
图 3.5 Chimera 技术示意图(二维示例)。

插值阶数作相应调整。与多块方法相比,Chimera 技术的优点是各套网格可完全独立地生成,无需顾及网格间的界面。其问题则在于,重叠区域内不能满足控制方程的守恒性质。文献 $ ^{[47, 48]} $ 近期报道了这一方法在非结构网格上的应用。

第二类是非结构网格,它在处理复杂几何形状方面具有最大的灵活性 $ [49] $。其主要优点是:原则上,无论计算域多么复杂,二维三角形网格或三维四面体网格都可自动生成。当然,实际应用中仍需适当设置若干参数,以获得质量良好的网格。此外,为精确解析边界层,二维情况下宜在固壁附近采用矩形单元,三维情况下宜采用棱柱或六面体单元 $ [50–60] $。这种混合网格的另一优点,是可以减少网格单元、边、面,以及可能的网格点数量。不过,对于几何结构要求较高的算例,混合网格的生成并非易事。即便如此,为复杂外形构造非结构混合网格所需的时间,仍显著少于生成结构多块网格所需的时间。如今,流动模拟对几何形状的保真要求迅速提高,能够快速生成网格并尽量减少用户干预,因而变得越来越重要,尤其是在工业环境中。非结构网格的另一个优点,是能够较为自然、顺畅地实现依赖流动解的网格加密与粗化。其不足之一,是流动求解器必须采用复杂的数据结构。这些结构使用间接寻址,因计算机硬件不同,计算效率会受到不同程度的影响。与结构格式相比,内存需求通常也更高。不过,能够在较短的完整计算周期内处理复杂几何问题,仍具有更大的价值。因此,几乎所有商业 CFD 软件厂商都转向非结构流动求解器,也就不足为奇了。文献 $ [61] $ 对非结构网格上的多种空间和时间离散方法作了详细评述。

生成网格后,接下来的问题是如何具体离散控制方程。如前所述,基本上可在有限差分、有限体积和有限元三类方法中选择。下面分别作简要介绍。

3.1.1 有限差分法

有限差分法是最早用于微分方程数值求解的方法之一,Euler 可能早在 1768 年就已首次使用。该方法直接作用于控制方程的微分形式,其基本原理是利用 Taylor 级数展开,对流动变量的导数进行离散。以下用一个例子说明。

假设需要计算标量函数 $ U(x) $ 在某一点 $ x_{0} $ 的一阶导数。将 $ U(x_{0} + \Delta x) $ 关于 x 展开为 Taylor 级数,得到

$$ U(x_{0}+\Delta x)=U(x_{0})+\Delta x\frac{\partial U}{\partial x}\Big|_{x_{0}}+\frac{\Delta x^{2}}{2}\frac{\partial^{2}U}{\partial x^{2}}\Big|_{x_{0}}+\cdots.\tag{3.1} $$

由此,U 的一阶导数可近似为

$$ \frac{\partial U}{\partial x}\bigg|_{x_{0}}=\frac{U(x_{0}+\Delta x)-U(x_{0})}{\Delta x}+\mathcal{O}(\Delta x).\tag{3.2} $$

这一近似是一阶的,因为与余项中最大项成正比的截断误差(记为 $ \mathcal{O}(\Delta x) $),随 $ \Delta x $ 的一次幂趋于零。有关精度阶数的讨论,见第 10 章。采用相同过程,也可以推导出精度更高的有限差分公式,以及更高阶导数的近似。

有限差分法的一个重要优点是简单。另一个优点,是容易获得高阶近似,从而实现高阶空间离散精度。但由于该方法要求使用结构网格,其适用范围受到明显限制。此外,有限差分法不能直接用于贴体曲线坐标,必须先将控制方程变换到笛卡尔坐标系;换言之,就是从物理空间变换到计算空间(图 3.2)。问题在于,坐标变换的 Jacobian 会出现在流动方程中(例如见第 A.1 节),必须以相容的方式进行离散,以避免引入额外的数值误差 $ [62] $。因此,有限差分法只能直接应用于较为简单的几何形状。目前,它用于湍流研究,也结合浸入边界单元(图 3.1b)用于生物学研究。更多细节可参见文献 $ [63] $,或有关偏微分方程求解的教科书。

3.1.2 有限体积法

有限体积法直接利用守恒定律,也就是 Navier–Stokes/Euler 方程的积分形式。McDonald $ [64] $ 首次将其用于二维无黏流动模拟。该方法先将物理空间划分为若干任意形状的多面体控制体,再对控制方程进行离散。随后,将式(2.19)中的面积分,近似为穿过控制体各个面的通量之和。空间离散的精度,取决于计算通量所采用的具体格式。

控制体相对于网格的形状和位置,有多种定义方式。基本上可区分为两类:

  • 单元中心格式(图 3.6a):将流动物理量存储于网格单元的形心,因此控制体与网格单元相同。

  • 单元顶点格式(图 3.6b):将流动变量存储于网格点。控制体可以是所有共享该网格点的单元的并集,也可以是以该网格点为中心构造的某个体积。前者称为重叠控制体,后者称为对偶控制体。

第 4、5 章讨论空间离散时,将介绍单元中心表述和单元顶点表述各自的优缺点。

有限体积法的主要优点,是直接在物理空间中进行空间离散,因此不会出现有限差分法中物理坐标系与计算坐标系之间变换所带来的问题。相较于有限差分法,它还有一个优点:灵活性很高,既容易用于结构网格,也容易用于非结构网格。因此,有限体积法尤其适合模拟复杂几何内部或周围的流动。

(a)
Image
Image
(b)
图 3.6 单元中心格式(a)与单元顶点格式(对偶控制体)(b)的控制体。

有限体积法直接离散守恒定律,因此其数值格式也保持质量、动量和能量守恒。由此产生了该方法的另一个重要特点:能够正确计算控制方程的弱解。不过,对于 Euler 方程,还必须满足一个附加条件,即熵条件。这是因为弱解并不唯一。熵条件可以防止膨胀激波等非物理现象出现;这类现象会使熵减少,违反热力学第二定律。守恒离散的另一个结果,是能够直接满足在解的间断处,例如激波或接触间断处,必须成立的 Rankine–Hugoniot 关系。

在某些条件下,可以证明有限体积法等价于有限差分法,或低阶有限元法。由于具备这些有利性质,有限体积法如今十分流行,应用广泛。后续章节将详细介绍。

3.1.3 有限元法

有限元法最初仅用于结构分析,由 Turner 等人 $ ^{[65]} $ 于 1956 年首次提出。约十年后,研究人员开始将其用于连续介质场方程的数值求解。不过,直到 20 世纪 90 年代初,有限元法才在 Euler 和 Navier–Stokes 方程的求解中逐渐普及。文献 $ ^{[66]} $ 对经典有限元法作了很好的入门介绍,文献 $ ^{[67-70]} $ 则讨论了其在流动问题中的应用。数年前,研究重点转向一种称为间断 Galerkin 格式的方法变体 $ [71-98] $。该格式的高阶表述尤其适用于气动声学问题。

通常用于求解 Euler/Navier–Stokes 方程的有限元法,首先将物理空间划分为离散单元:二维中多为三角形或四边形,三维中多为四面体或六面体,因此需要生成非结构网格。根据单元类型和所需精度,在单元边界和/或内部规定一定数量的点,并在这些点求解流动问题。点的总数乘以未知量的数量,即确定了自由度数。还需定义形函数,用以表示各单元内部解的变化。在实际实现中,采用线性单元或高阶单元;线性单元仅使用网格节点。相应的形函数为线性分布或高阶分布,

其在所属单元之外的值为零。由此,在规则网格上可获得二阶或更高阶精度的解表示。

在有限元法中,必须将控制方程从微分形式变换为等价的积分形式。可采用两种方式。第一种基于变分原理,即寻找使某个泛函取得极值的物理解。第二种称为加权余量法或弱表述,要求残差在物理计算域上的加权平均恒等于零。残差可视为解的近似误差。弱表述具有与守恒定律有限体积离散相同的优点:能够处理激波等间断。因此,通常优先采用弱表述,而非变分方法。间断 Galerkin 格式与经典有限元法的区别之一,在于其质量矩阵仅定义于生成它的局部单元。由此,可以用简单的显式方法,对未知解系数进行时间推进。

有限元法采用积分表述和非结构网格,而这两点对复杂几何内部或周围的流动均十分有利,因此具有吸引力。除气动声学外,该方法也特别适合处理非 Newton 流体。有限元法具有十分严谨的数学基础,尤其是在椭圆型和抛物型问题中。虽然在某些情况下,可以证明它与有限体积离散在数学上等价,但所需数值计算量明显更大。这可能解释了为什么有限体积法在一般流动问题中更受欢迎。不过,两种方法有时也会结合使用,尤其是在非结构网格上。例如,边界处理和黏性通量离散通常会“借用”有限元法。

3.1.4 其他离散方法

还有少数数值格式,在某些情况下优于前述方法。下面简要介绍其中三种。

谱元法

第一个例子是谱元法 $ [99–107] $。这种方法将有限元技术在几何处理上的灵活性,与谱格式的高阶空间精度(例如十阶)和快速收敛特性结合起来 $ [101] $。该方法用高阶多项式(通常为拉格朗日插值多项式)表示解,并结合标准的 Galerkin 有限元法或加权余量法。谱元法适用于能够保证高阶正则性的特定问题,或高阶正则性普遍存在的问题,例如不可压缩流体力学问题;它尤其适用于涡流。谱元法的主要优势在于,对对流算子的近似不引入扩散和色散,并能很好地近似对流—扩散边界层。只要满足高阶正则性的条件,该方法就能处理几何和物理上复杂的问题。除适用范围较窄外,谱元法的主要缺点是计算量很大,例如远高于有限体积法。

格子 Boltzmann 方法

格子 Boltzmann 方法(LBM)$ [108–129] $ 于 20 世纪 90 年代末作为格子气方法的后继方法提出。随着高性能大规模并行计算系统日益普及,该方法也逐渐受到重视。前面讨论的方法以宏观物理量(即质量、动量和能量)的守恒为基础,而 LBM 将流体建模为由虚拟粒子组成的系统。这些粒子在由网格单元构成的离散格子上反复进行微观迁移和碰撞。该方法基于动理学理论,可视为连续 Boltzmann 方程的一种特殊离散形式。通常,LBM 首先在局部以松弛过程计算碰撞的作用(碰撞步),然后将计算结果——离散分布函数——沿特征方向传递到相应的相邻节点(迁移步)。宏观场通过求分布函数的离散矩得到,压力场也可直接由密度分布获得。应当指出,根据 Chapman–Enskog 理论,可以从 LBM 算法恢复连续性方程和 Navier–Stokes 方程。近期的一篇文献还揭示了 LBM 与一种有限差分方法之间值得关注的相似性,并比较了两者的计算效率 $ [123] $。

由于具有粒子动理学机制和局部动力学特性,LBM 相比传统 CFD 方法具有多项优势,尤其是在处理复杂边界或引入微观相互作用时。例如,多相/多组分流动的模拟非常困难,因为其中存在运动和变形的界面,以及更为复杂的相间过渡边界。然而,LBM 的形式允许通过修改碰撞算子,以相对简单且一致的方式处理这类问题。相分离由粒子动力学自动产生,无须对界面进行特殊处理。因此,LBM 已成功用于模拟界面不稳定性、气泡/液滴动力学、固体表面的润湿以及液滴的电流体动力学变形。它还应用于多孔介质中的多相流、气动声学和微流体领域。

“碰撞—迁移”算法使 LBM 非常适合采用区域分解技术进行高效并行计算。另一方面,传统 LBM 本质上采用时间推进(物理时间或伪时间),因而往往需要较慢的收敛过程才能得到所求的解(可能是定常解)。此外,迄今尚无适用于强可压缩性效应流动的 LBM 表述。这意味着,LBM 不能模拟跨声速或超声速流动。

无网格方法

另一种近来受到关注的离散格式是无网格方法 $ [130-136] $。该方法仅使用点云进行空间离散,不要求像传统结构或非结构网格格式那样将点连接成网格。无网格方法以笛卡尔坐标系中控制方程的微分形式为基础,利用某一点周围指定数量的邻点,通过最小二乘重构确定流动变量的梯度。由于无须计算坐标变换、面面积或体积,无网格方法既不是有限差分法,也不是有限体积法或有限元法;可以将其视为有限差分法与有限元法的混合。其主要优势在于,处理复杂构型绕流时具有类似非结构网格方法的灵活性,而且可以在适当位置布置或加密点(或点云)。例如,计算梯度时,很容易只选择特征方向上的邻点。然而,仍有一个尚未解决的问题:尽管无网格方法求解的是 Euler 方程或 Navier–Stokes 方程的守恒形式,却还不完全清楚它是否真正保证质量、动量和能量守恒。

无论选择哪一种空间离散格式,都必须确保格式具有相容性,即当网格充分加密时,解能够收敛到离散方程的解。因此,检查网格加密(例如将网格点数加倍)后解的变化幅度非常重要。如果解仅有微小改善,就称其为网格收敛解。另一个相当明显的要求是:离散格式应具有与待求流动问题相适应的精度阶次。有时,为获得更快的收敛,会放弃这一要求,尤其是在工业应用中,理由是“较差的解总比没有解好”。这当然是一种十分危险的做法。第 10 章将再次讨论精度、稳定性和相容性问题。

3.1.5 中心格式与迎风格式

到目前为止,我们只讨论了空间离散的基本选择。然而,在上述三种主要方法——有限差分法、有限体积法和有限元法——之中,都存在多种具体的空间离散格式。讨论这些格式时,宜区分式(2.19)中的对流通量 $ \vec{F}_{mathrmc} $ 与黏性通量 $ \vec{F}_{v} $ 的离散。由黏性通量的物理性质决定,唯一合理的做法是采用中心差分(中心平均)离散。因此,在结构网格上离散黏性通量十分直接。在非结构三角形或四面体网格上,即使采用有限体积格式,也最好用 Galerkin 有限元方法近似黏性通量 [137]。对于非结构混合网格,情况更复杂,此时采用修正的梯度平均更为合适 [138–142]。

不过,真正丰富多样的是对流通量的离散方法。为对各类方法进行分类,下面只讨论为有限体积法构造的格式;其中的大多数概念同样可以直接应用于有限差分法或有限元法。

中心格式

第一类格式完全基于中心差分公式或中心平均,称为中心格式。其基本思想是,对控制体某一侧面左右两侧的守恒变量取平均,以计算该面上的通量。由于中心格式无法识别和抑制解的奇偶解耦(即离散方程产生两个相互独立的解),必须加入所谓的人工耗散以稳定计算;这一名称源于它与黏性项的相似性。最著名的实现由 Jameson 等人提出 $ ^{[143]} $。在结构网格上,它将二阶差分与四阶差分混合,并用对流通量 Jacobian 的最大特征值进行尺度调整。在非结构网格上,则采用不除以网格间距的 Laplace 算子与双调和算子的组合 $ [144] $。对每个方程采用不同的尺度因子,可以显著改善该格式;这种方法称为矩阵耗散格式 $ [145] $。应当指出,在非结构混合单元网格上,显式 Runge–Kutta 时间推进与传统中心格式结合时可能失稳 $ [146] $。

迎风格式

另一方面,还存在一些更先进的空间离散格式,它们根据 Euler 方程的物理性质构造。由于这些格式区分上游和下游的影响(即波的传播方向),因此称为迎风格式。大致可将它们分为四类:

  • 通量矢量分裂格式;

  • 通量差分分裂格式;

  • 总变差减小(TVD)格式;

  • 波动分裂格式。

下面分别作简要介绍。

通量矢量分裂格式

一类通量矢量分裂格式按照某些特征变量的符号,将对流通量矢量分成两部分。这些特征变量通常与对流通量 Jacobian 的特征值相似,但并不完全相同。随后,用偏迎风差分分别离散这两部分通量。此类最早的通量矢量分裂格式由 Steger 和 Warming $ [147] $ 以及 Van Leer $ [148] $ 于 20 世纪 80 年代初分别提出。第二类通量矢量分裂格式将通量矢量分为对流部分和压力(声学)部分。Liou 等人的 AUSM(对流迎风分裂法)$ [149, 150] $ 以及 Jameson 的 CUSP(对流迎风分裂压力)格式 $ [151, 152] $ 都采用了这一思想。类似方法还包括 Edwards 提出的低扩散通量分裂格式(LDFSS)$ [153] $,以及 Rossow 的基于 Mach 数的对流—压力分裂(MAPS)格式 $ [154, 155] $。第二类通量矢量分裂格式近来更加流行,主要原因是它能改善剪切层的分辨能力,而计算量仍较适中。通量矢量分裂格式的另一项优点是,相较于通量差分分裂格式和 TVD 格式,它更容易推广到真实气体流动。下面还将讨论真实气体模拟。

通量差分分裂格式

第二类即通量差分分裂格式,以界面两侧不连续状态下局部一维 Euler 方程的解为基础。这对应于 Riemann(激波管)问题;界面两侧的值通常称为左状态和右状态。在两个控制体之间的界面上求解 Riemann 问题的思想,最早由 Godunov $ [156] $ 于 1959 年提出。为降低精确求解 Riemann 问题所需的计算量,研究者构造了近似 Riemann 求解器,例如 Osher 和 Solomon $ [157] $ 以及 Roe $ [158] $ 的方法。Roe 求解器因能很好地分辨边界层并清晰捕捉激波,如今得到广泛应用。它容易在结构和非结构网格上实现 $ [159] $。

其他基于 Riemann 问题近似解的方法称为 Harten–Lax–van Leer(HLL)求解器。其中,Einfeldt 提出的 HLLE 变体 $ [160] $ 以流体间断的稳定性理论为基础,考虑界面处最小和最大的波速。Toro 等人提出了另一种称为 HLLC(HLL contact,即含接触间断的 HLL)的变体 $ [161] $,旨在恢复 HLL 格式中的接触间断;随后,其他研究者对它进行了研究和修改 $ [162–164] $。Nishikawa 和 Kitamura 又提出了旋转混合 Riemann 求解器 $ [165] $,以同时克服原始 Roe 格式的 carbuncle(激波鼓包)现象 1 和 HLL 求解器过大的扩散。应当指出,对 Roe 格式作一项简单修改——Harten 熵修正——即可较容易地避免激波鼓包问题。第 4.3.3 节讨论 Roe 格式时将介绍这一修正。考虑到这一点,很难看出用某种 HLL 变体取代 Roe 迎风格式能带来什么实际收益。

1 激波鼓包现象是指强弓形激波前方沿驻点线增长的扰动,其原因是格式未能识别声速点。

TVD 格式

Harten $ [166] $ 于 1983 年首次提出 TVD 格式的思想。其基本概念是防止流动解中产生新的极值。TVD 格式的主要条件为:最大值不得增加,最小值不得减小,且不得产生新的局部极值;满足这些条件的格式称为保单调格式。因此,具有 TVD 性质的离散方法能够分辨激波,而不产生虚假的解振荡。TVD 格式通常实现为对流通量的平均值加上一个额外耗散项。耗散项可以依赖特征速度的符号,也可以不依赖;前者称为迎风 TVD 格式 $ [167] $,后者称为对称 TVD 格式 $ [168] $。经验表明,应优先选用迎风 TVD 格式,因为它对激波和边界层的分辨能力优于对称 TVD 格式。TVD 格式的一项缺点是,不容易推广到二阶以上的空间精度。

波动分裂格式

最后一类——波动分裂格式——实现了真正的多维迎风处理。其目标是准确分辨那些与网格不对齐的流动特征。这是它相对于上述全部迎风格式的一项显著优势,因为后者只依据网格单元的方向分裂方程。在波动分裂方法中,流动变量定义在网格节点上。首先,在二维三角形或三维四面体单元上,根据通量平衡计算中间残差;然后以偏迎风的方式将单元残差分配给节点,再使用节点值更新解。对于方程组(Euler 或 Navier–Stokes 方程),必须将单元残差分解为标量波。由于二维和三维情况下分解并不唯一,研究者提出了多种方法,从 Roe 的波模型 [169, 170]、Sidilkover 的代数格式 [171],到更先进的特征分解方法 [172–176]。尽管相较于按维分裂的 Riemann、TVD 等求解器具有上述优势,波动分裂方法迄今仍仅用于研究程序;原因在于其复杂性、较大的计算量以及收敛问题。

解的重构

上述全部格式都需要知道控制体边界上的密度、速度分量和压力等流动变量;有时还需要分别给出界面两侧的左状态和右状态。由于每个控制体内已知的仅是流动变量的平均值,必须采用某种插值过程,这一过程称为重构。显然,解的重构精度决定了离散格式的空间精度阶次。

一阶与二阶格式

对于中心格式,对公共界面左右两侧控制体的平均值作简单算术平均,即可达到二阶精度(参见第 4.3.1 节)。对于迎风格式,若将左、右状态分别设为对应控制体的平均值,则得到一阶精度格式。要在结构网格上达到二阶空间精度,需要采用单边差分(见第 4.3 节式(4.46))。在非结构网格上实现二阶迎风格式,则首先需要计算流动变量的梯度,再利用梯度和平均值,通过方向导数插值得到左、右状态(第 5.3.3 节)。

ENO/WENO 格式

采用 ENO(本质无振荡)或 WENO(加权 ENO)离散格式 $ [177–197] $,可以获得高于二阶的精度。这些方法采用方向可变的离散模板,计算控制体面上变量的不同重构(插值)多项式。ENO 格式选择振荡最小的多项式;WENO 方法则对各多项式加权组合,权重与振荡程度成反比。这样,在光滑区域内可获得高阶空间精度,同时又不会过度抹平真实间断。文献 $ [185] $ 对 WENO 格式作了很好的介绍。

中心格式与迎风格式的比较

读者可能会问,各种空间离散方法的优缺点是什么?一般而言,中心格式的计算量远小于迎风格式,因此每次计算所需的 CPU 时间也显著较少。另一方面,迎风格式能够比中心格式更准确地捕捉间断。由于数值扩散较低,迎风格式还可以用更少的网格点分辨边界层;特别是 Roe 通量差分分裂格式和 CUSP 通量矢量分裂格式,能够十分准确地计算边界层。迎风格式的缺点在达到二阶或更高阶空间精度时显现:为了防止强间断附近产生虚假振荡,必须采用所谓的限制器函数(简称限制器)。限制器在光滑流动区域内的意外切换,会使迭代格式的收敛停滞。Venkatakrishnan $ [198-200] $ 提出了一种补救方法,对大多数实际问题效果令人满意,不过必须允许解中存在微小扰动。限制器的另一项缺点是计算量较大,尤其是在非结构网格上。后面详细讨论迎风离散格式时,将再次讨论限制器。

真实气体流动的迎风格式

针对真实气体模拟,尤其是化学反应流,研究者提出了多种迎风离散格式的扩展。对于处于热力学和化学平衡的流体,Van Leer 通量矢量分裂格式 $ [148] $ 和 Roe 近似 Riemann 求解器 $ [158] $ 的修改形式见文献 $ [201–203] $。对于更复杂的化学与热力学非平衡流,文献 $ [204–209] $ 及其所引文献给出了这两种迎风方法的表述。其中,文献 $ [203, 207, 208] $ 对迎风离散格式所采用的方法作了很好的概述。文献 $ [210] $ 汇总了化学反应流的控制方程,以及迎风格式可能需要的 Jacobian 矩阵和变换矩阵;文献 $ [211] $ 则给出了平衡真实气体的相应结果。