0 引言

近几十年来，计算流体力学(CFD)越来越广泛地应用于航空航天，可以说CFD改变了航空航天的设计过程。计算机性能的提高和CFD技术的发展，是CFD大规模应用的“助推器”。据统计，高性能计算机(HPC)计算速度每十年将提高三个量级^[1]。然而，现今的CFD方法在建立过程中很多都没有考虑到并行计算的特点，并行水平较低，已经不能适应HPC发展的要求。NASA在2030年的CFD展望^[2]中指出，目前鲁棒性的CFD流动解算器非常缺乏可扩展性，极少有可以有效利用超过O(1000)个内核的应用能力，为此需要开发并行程度很高的CFD算法，用以平衡计算与通讯。

在CFD并行计算中，普遍采用基于网格分区的并行策略。在CFD时间推进中，为了不受稳定性条件的限制以及出于加速流动模拟收敛的目的，一般采用各种隐式计算方法。然而，传统的以LU-SGS(Lower-Upper Symmetric Gauss-Seidel)方法^[3]为代表的隐式近似因子分解(AF)类方法应用于并行计算时，存在两个问题。一方面，隐式算法中相邻网格单元之间较强的数据依赖使得并行求解较为困难^[4]。另一方面，如果直接在各网格分区中采用隐式算法，则随着分区数的增加算法收敛性能将恶化^[5]。

线松弛方法是一类隐式求解的迭代算法。MacCormack(1985)^[6]提出求解Navier-Stokes方程的GSLR(Gauss-Seidel Line Relaxation)方法。Candler(1988)^[7]将其应用到热化学非平衡流动的数值模拟中，它的特点是收敛很快。GSLR方法的思想是：在黏性流动中，相较于流向和周向，流场梯度在物面法向更大，流动在物面法向耦合更紧密，数值求解时在物面法向把隐式部分的对角项和非对角项耦合求解才能反映流动的物理本质。GSLR方法在流向和周向采用Gauss-Seidel扫描，不利于并行计算。NASA Ames研究中心的Wright(1997)^[8]，以GSLR为出发点，把Gauss-Seidel扫描改为松弛过程，消除相邻物面法向网格线流场变量之间的数据依赖，既保留了GSLR收敛快的特点，又使它有利于并行计算，这就是DPLR(Data-Parallel Line Relaxation)方法。

Wright发明并在单块结构网格上应用DPLR方法^[8]，后以此为基础研制著名的气动热力学软件DPLR^[9-10]，这是一款基于多块结构网格的CFD软件。但是DPLR方法在多块网格上的性能如何，网格分区对算法的收敛性能会产生多大影响，这些问题未在公开文献中报道。由于DPLR需要结合隐式边界条件使用以及多块对接结构网格中存在大量的对接边界，本文首先提出一种隐式对接边界条件的处理方法，为把DPLR应用于多块结构网格奠定基础。然后，以二维和三维球头绕流为算例，分析比较了DPLR与LU-SGS的收敛速度，研究了网格分块方式对算法稳定性和收敛性能的影响。最后进行了返回舱外形绕流的数值模拟，给出返回舱外形在马赫数9.96、雷诺数2.04×10⁵、迎角10°~30°的气动力特性，并与风洞试验结果对比，考察算法在复杂外形飞行器绕流数值模拟中的性能。

1 数值计算方法 1.1 控制方程

为叙述基于多块对接结构网格的DPLR隐格式实现途径，以二维无量纲Navier-Stokes方程为例，说明本文的计算方法，三维控制方程的计算方法与之类似。控制方程^[11]可写为：

$ \frac{{\partial \mathit{\boldsymbol{Q}}}}{{\partial t}} + \frac{{\partial \mathit{\boldsymbol{E}}}}{{\partial x}} + \frac{{\partial \mathit{\boldsymbol{F}}}}{{\partial y}} = \frac{{\partial {\mathit{\boldsymbol{E}}_v}}}{{\partial x}} + \frac{{\partial {\mathit{\boldsymbol{F}}_v}}}{{\partial y}} $

(1)

式中，Q为守恒变量矢量，E、F为无黏通量矢量，E_v、F_v为黏性通量矢量，它们的表达式详见文献[11]。

1.2 空间项离散

应用有限体积法，以网格单元为控制体，对控制方程进行离散。在有限体积法中，流场变量储存在单元格心，在单元内部均匀分布。本文记流向为ξ，下标索引为i；法向为η，下标索引为j。把控制方程(1)在控制体(i, j)积分，并应用Gauss-Green定理，则有：

$ \frac{{\partial \mathit{\boldsymbol{Q}}}}{{\partial t}} = - \frac{1}{{{V_{i,j}}}}\mathit{\boldsymbol{RH}}{\mathit{\boldsymbol{S}}_{i,j}} $

(2)

$ \begin{array}{l} \mathit{\boldsymbol{RH}}{\mathit{\boldsymbol{S}}_{i,j}} = {\mathit{\boldsymbol{E}}_{i + \frac{1}{2},j}}{S_{i + \frac{1}{2},j}} - {\mathit{\boldsymbol{E}}_{i - \frac{1}{2},j}}{S_{i - \frac{1}{2},j}} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;{\mathit{\boldsymbol{F}}_{i,j + \frac{1}{2}}}{S_{i,j + \frac{1}{2}}} - {\mathit{\boldsymbol{F}}_{i,j - \frac{1}{2}}}{S_{i,j - \frac{1}{2}}} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;{\mathit{\boldsymbol{E}}_{v,i + \frac{1}{2},j}}{S_{i + \frac{1}{2},j}} - {\mathit{\boldsymbol{E}}_{v,i - \frac{1}{2},j}}{S_{i - \frac{1}{2},j}} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;{\mathit{\boldsymbol{F}}_{v,i,j + \frac{1}{2}}}{S_{i,j + \frac{1}{2}}} - {\mathit{\boldsymbol{F}}_{v,i,j - \frac{1}{2}}}{S_{i,j - \frac{1}{2}}} \end{array} $

(3)

式中，RHS表示残差，E、F分别为通过垂直于ξ、η方向单元面的无黏通量，E_v、F_v分别为相应的黏性通量。V是控制体单元的体积，S是单元面的面积。

在本文中，无黏通量采用Steger-Warming矢通量分裂^[12]以及NND格式^[13]计算；黏性通量采用中心差分格式计算。

1.3 时间项离散

对方程(2)取n+1时间层，左侧的时间导数进行后向差分，在离散中需要对无黏和黏性通量线性化。无黏通量E的线性化为:

$ {\mathit{\boldsymbol{E}}^{n + 1}} \approx {\mathit{\boldsymbol{E}}^n} + {\mathit{\boldsymbol{A}}^n}\delta {\mathit{\boldsymbol{Q}}^{n + 1}} $

(4)

黏性通量E_v根据Tysinger和Caughey的方法^[14]线性化为:

$ \mathit{\boldsymbol{E}}_v^{n + 1} \approx \mathit{\boldsymbol{E}}_v^n + \frac{\partial }{{\partial \xi }}{\left( {\mathit{\boldsymbol{L}}\delta \mathit{\boldsymbol{Q}}} \right)^n} $

(5)

式中，A是无黏通量的Jacobian矩阵，$ \mathit{\boldsymbol{A}}{\rm{ = }}\frac{{\partial \mathit{\boldsymbol{E}}}}{{\partial \mathit{\boldsymbol{Q}}}}$; L是黏性通量的线性化矩阵，L=$ \frac{{\partial {\mathit{\boldsymbol{E}}_\mathit{v}}}}{{\partial {\mathit{\boldsymbol{Q}}_\xi }}}$; δQⁿ⁺¹=Qⁿ⁺¹-Qⁿ。则隐式离散方程为:

$ \begin{array}{l} {{\mathit{\boldsymbol{\hat A}}}_{i,j}}\delta \mathit{\boldsymbol{Q}}_{i,j}^{n + 1} + {{\mathit{\boldsymbol{\hat B}}}_{i,j}}\delta \mathit{\boldsymbol{Q}}_{i,j + 1}^{n + 1} + {{\mathit{\boldsymbol{\hat C}}}_{i,j}}\delta \mathit{\boldsymbol{Q}}_{i,j - 1}^{n + 1}\\ \;\;\;\;\;\; + {{\mathit{\boldsymbol{\hat D}}}_{i,j}}\delta \mathit{\boldsymbol{Q}}_{i + 1,j}^{n + 1} + {{\mathit{\boldsymbol{\hat E}}}_{i,j}}\delta \mathit{\boldsymbol{Q}}_{i - 1,j}^{n + 1} = - \frac{{\Delta t}}{{{V_{i,j}}}} \cdot \mathit{\boldsymbol{RHS}}_{i,j}^n \end{array} $

(6)

$ \begin{array}{l} {{\mathit{\boldsymbol{\hat A}}}_{i,j}} = I + \frac{{\Delta t}}{{{V_{i,j}}}}\left\{ {\mathit{\boldsymbol{\tilde A}}_{i + \frac{1}{2},j}^{n + }{S_{i + \frac{1}{2},j}} + \mathit{\boldsymbol{\tilde B}}_{i,j + \frac{1}{2}}^{n + }{S_{i,j + \frac{1}{2}}} - } \right.\\ \;\;\;\;\;\;\;\;\;\left. {\mathit{\boldsymbol{\tilde A}}_{i - \frac{1}{2},j}^{n - }{S_{i - \frac{1}{2},j}} - \mathit{\boldsymbol{\tilde B}}_{i,j - \frac{1}{2}}^{n - }{S_{i,j - \frac{1}{2}}}} \right\}; \end{array} $

$ {{\mathit{\boldsymbol{\hat B}}}_{i,j}} = \frac{{\Delta t}}{{{V_{i,j}}}}{S_{i,j + \frac{1}{2}}}\mathit{\boldsymbol{\tilde B}}_{i,j + \frac{1}{2}}^{n - };{{\mathit{\boldsymbol{\hat C}}}_{i,j}} = - \frac{{\Delta t}}{{{V_{i,j}}}}{S_{i,j - \frac{1}{2}}}\mathit{\boldsymbol{\tilde B}}_{i,j - \frac{1}{2}}^{n + } $

$ {{\mathit{\boldsymbol{\hat D}}}_{i,j}} = \frac{{\Delta t}}{{{V_{i,j}}}}{S_{i + \frac{1}{2},j}}\mathit{\boldsymbol{\tilde A}}_{i + \frac{1}{2},j}^{n - };{{\mathit{\boldsymbol{\hat E}}}_{i,j}} = - \frac{{\Delta t}}{{{V_{i,j}}}}{S_{i - \frac{1}{2},j}}\mathit{\boldsymbol{\tilde A}}_{i - \frac{1}{2},j}^{n + } $

(7)

其中:

$ \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\tilde A}}}^ \pm } = {\mathit{\boldsymbol{A}}^ \pm } \pm L}&{{{\mathit{\boldsymbol{\tilde B}}}^ \pm } = {\mathit{\boldsymbol{B}}^ \pm } \pm N} \end{array} $

(8)

A^±、B^±分别是ξ、η方向的经过分裂的无黏Jacobian矩阵，L和N分别是ξ、η方向黏性通量的线性化矩阵，Δt是时间步长。

1.4 DPLR格式

在隐式离散方程(6)中，$ {{\mathit{\boldsymbol{\hat A}}}_{i, j}}$是对角项系数矩阵，$ {{\mathit{\boldsymbol{\hat B}}}_{i, j}}、{{\mathit{\boldsymbol{\hat C}}}_{i, j}}、{{\mathit{\boldsymbol{\hat D}}}_{i, j}}、{{\mathit{\boldsymbol{\hat E}}}_{i, j}}$是非对角项系数矩阵。在LU-SGS等AF类隐格式中，对这些矩阵采用最大特征值分裂算法实现系数矩阵的对角化以避免矩阵求逆。Wright指出^[8]，这种做法虽然计算量较小但会导致流场信息的缺失，从而影响算法的收敛速度。为此，在DPLR格式中，这些矩阵的计算要使用精确的经过Steger-Warming分裂的无黏Jacobian矩阵以及黏性线性化矩阵。

设η方向为物面法向，则DPLR格式^[8]可写为: 对于m=1，m_max

$ \begin{array}{l} {{\mathit{\boldsymbol{\hat A}}}_{i,j}}\delta \mathit{\boldsymbol{Q}}_{i,j}^{\left( m \right)} + {{\mathit{\boldsymbol{\hat B}}}_{i,j}}\sigma \mathit{\boldsymbol{Q}}_{i,j + 1}^{\left( m \right)} + {{\mathit{\boldsymbol{\hat C}}}_{i,j}}\delta \mathit{\boldsymbol{Q}}_{i,j - 1}^{\left( m \right)}\\ \;\; = - {{\mathit{\boldsymbol{\hat D}}}_{i,j}}\delta \mathit{\boldsymbol{Q}}_{i + 1,j}^{\left( {m - 1} \right)} - {{\mathit{\boldsymbol{\hat E}}}_{i,j}}\delta \mathit{\boldsymbol{Q}}_{i - 1,j}^{\left( {m - 1} \right)} - \frac{{\Delta t}}{{{V_{i,j}}}} \cdot \mathit{\boldsymbol{RHS}}_{i,j}^n \end{array} $

(9)

则：δQ_{i, j}ⁿ⁺¹=δQ_{i, j}^(m_max)

其中，松弛步初值Q_{i, j}⁽⁰⁾=0，m_max是松弛次数。文献[8]认为，m_max=4是最佳松弛次数。

从式(9)看出，在时间方向每推进一步，DPLR格式需要求解m_max次三对角块矩阵方程组，在本文中采用追赶法^[15]求解。根据追赶法的特点，在m_max次松弛中只需要求逆左端系数矩阵$ {{\mathit{\boldsymbol{\hat A}}}_{i, j}}、{{\mathit{\boldsymbol{\hat B}}}_{i, j}}、{{\mathit{\boldsymbol{\hat C}}}_{i, j}}$一次。DPLR格式需要在每个网格存储5个系数矩阵(三维问题是7个系数矩阵)，存储量较大。

1.5 边界条件

DPLR需要和隐式边界条件结合起来使用^[8-9]。本文采用虚拟网格方法^{[6-7, 16]}处理边界条件。虚拟网格不由网格软件生成，而是在代码中根据不同边界类型逻辑生成。图 1和图 2分别是物面边界和对接边界虚拟网格示意图。

文献[16]认为，虚拟网格上流场变量设置的原则是边界面上的无黏和黏性通量能够正确计算。在GSLR方法中，MacCormack认为^[6]可以在隐式计算中把边界条件包含进来，这就是隐式边界条件。因此，DPLR格式的边界条件分为显式和隐式两部分，显式部分是为了正确计算边界面上的通量，隐式部分则体现了边界条件对隐式求解的影响。这和LU-SGS方法是不同的，LU-SGS隐式部分通常不包含边界条件。

参照图 1，设边界面两侧网格中心的某流场变量分别为φ₁和φ₀, 设边界面上的物理量为φ_b，下标“1”表示内点，“0”表示虚拟点，“b”表示边界面。显式边界条件可写为通式:

$ {\varphi _0} = 2{\varphi _b} - {\varphi _1} $

(10)

对于整个求解区域的外边界，显示边界条件可应用式(10)。例如，对于超声速入口，直接给出φ₀=φ_∞，其中“∞”表示自由来流；对于超声速出口，φ_b=φ₁；对于绝热壁无滑移物面，P_b=P₁，u_b=0，v_b=0，T_b=T₁；对于等温壁无滑移物面，P_b=P₁，u_b=0，v_b=0，T_b=T_w。其中，P、u、v、T分别表示流场压力、x方向速度分量、y方向速度分量和温度。对于相邻网格块之间的对接边界(见图 2)，虚拟网格上的流场变量φ₀从与之对接的网格块内点直接得到，φ₀和φ₁之间没有类似式(10)的函数关系。

对式(10)微分，得到dφ₀和dφ₁的函数关系，经过进一步推导可得到边界两侧守恒变量微分之间的函数关系:

$ {\rm{d}}{\mathit{\boldsymbol{Q}}_0} = \mathit{\boldsymbol{H}}{\rm{d}}{\mathit{\boldsymbol{Q}}_1} $

(11)

其中H在文献[16]中称为折叠矩阵，式(11)就是隐式边界条件。把式(11)代入离散方程(9)，就可得到靠近边界内点的离散方程。例如，设某边界位于j=1/2位置，则在控制体(i, 1)处的离散方程为:

$ \begin{array}{l} \left( {{{\mathit{\boldsymbol{\hat A}}}_{i,1}} + {{\mathit{\boldsymbol{\hat C}}}_{i,1}} \cdot {\mathit{\boldsymbol{H}}_{i,\frac{1}{2}}}} \right)\delta \mathit{\boldsymbol{Q}}_{i,1}^{\left( m \right)} + {{\mathit{\boldsymbol{\hat B}}}_{i,1}}\delta \mathit{\boldsymbol{Q}}_{i,2}^{\left( m \right)}\\ = - {{\mathit{\boldsymbol{\hat D}}}_{i,1}}\delta \mathit{\boldsymbol{Q}}_{i + 1,1}^{\left( {m - 1} \right)} - {{\mathit{\boldsymbol{\hat E}}}_{i,1}}\delta \mathit{\boldsymbol{Q}}_{i - 1,1}^{\left( {m - 1} \right)} - \frac{{\Delta t}}{{{V_{i,1}}}} \cdot \mathit{\boldsymbol{RHS}}_{i,1}^n \end{array} $

(12)

从式(12)可看出，隐式边界条件把边界条件的影响包含进离散方程的隐式部分。

但是对于对接边界，边界面两侧的流场变量φ₀和φ₁不存在显式的函数关系式(10)，也无法推导出隐式边界条件式(11)，那么在对接边界处怎样把边界条件对隐式求解的影响考虑进来呢？本文提出一种处理方法，还是以j=1/2边界(假定是对接边界)为例：

$ \begin{array}{l} {{\mathit{\boldsymbol{\hat A}}}_{i,1}}\delta \mathit{\boldsymbol{Q}}_{i,1}^{\left( m \right)} + {{\mathit{\boldsymbol{\hat B}}}_{i,1}}\delta \mathit{\boldsymbol{Q}}_{i,2}^{\left( m \right)}\\ \;\; = - {{\mathit{\boldsymbol{\hat C}}}_{i,1}}\delta \mathit{\boldsymbol{Q}}_{i,0}^{\left( {m - 1} \right)} - {{\mathit{\boldsymbol{\hat D}}}_{i,1}}\delta \mathit{\boldsymbol{Q}}_{i + 1,1}^{\left( {m - 1} \right)} - {{\mathit{\boldsymbol{\hat E}}}_{i,1}}\delta \mathit{\boldsymbol{Q}}_{i - 1,1}^{\left( {m - 1} \right)}\\ \;\;\;\; - \frac{{\Delta t}}{{{V_{i,1}}}} \cdot \mathit{\boldsymbol{RHS}}_{i,1}^n \end{array} $

(13)

从(13)式看出，涉及到虚拟点上的项$ {{\mathit{\boldsymbol{\hat C}}}_{i, 1}}$δQ_{i, 0}^(m-1)采用上一松弛步的变量增量δQ_{i, 0}^(m-1)计算。而δQ_{i, 0}^(m-1)可通过对接边界处网格的对应关系直接从与之对接的另一个块的内点得到，对于本松弛步(m步)是已知量。这样就把对接边界条件对隐式计算的影响考虑进来。

2 算例验证及讨论 2.1 计算条件及网格

为了验证本文的计算方法和程序、比较DPLR和LU-SGS的收敛性能以及研究网格分区方式对算法稳定性和收敛性能的影响，本文计算两个算例。算例1是二维球头绕流，算例2是三维球头绕流，计算条件如表 1所示。算例1的网格如图 3(a)所示，网格数为121×81(流向×法向)。算例2的网格如图 3(b)所示，总共3个网格块，每个块的网格数都是16×101×16(流向×法向×周向)。两个算例法向最小网格间距都设置为使得网格雷诺数Re_Δn=2。

2.2 计算结果正确性验证

对于二维球头绕流，本文把计算得到的物面无量纲压力和热流(与驻点值的比值)分布与Wieting试验结果^[17]比较(见图 4)，两者吻合较好。对于三维球头绕流，本文把计算结果和Coleman的试验结果^[18]比较，并如图 5所示。由图可见，本文计算的三维球头物面压力(与自由来流值的比值)、热流分布与试验吻合较好。

2.3 DPLR稳定性和收敛性能以及网格分区方式影响的分析与讨论

收敛快是DPLR算法的特点之一，Wright在这方面已经做了详尽的研究^[8]。对于二维球头和三维球头算例，本文把DPLR与LU-SGS算法的收敛曲线进行对比，分别如图 6和图 7所示，其中横坐标是CPU时间(Intel Core i7-4790单核)。由于DPLR不是无条件稳定格式^[8]，在应用DPLR格式中这两个算例的CFL都取为能够使计算正常进行并保证收敛的最大CFL数，其中二维球头绕流算例CFL数取为3000，三维球头绕流算例CFL数取为3500。LU-SGS是无条件稳定格式，其CFL在两个算例中都取为1.0×10⁶。此外，基于Wright的经验^{[8, 10]}在DPLR计算时采用全局时间步长法；而在LU-SGS计算时应用当地时间步长法以加速收敛。从DPLR和LU-SGS收敛性能的对比看，DPLR格式收敛更快，二维球头绕流收敛到机器零的时间相当于LU-SGS的35%，三维球头绕流收敛到机器零的时间相当于LU-SGS的22%。从残差曲线的变化历程看，在迭代的初始阶段，LU-SGS方法由于单步计算量小而收敛较快，然而在迭代后期，DPLR方法就显示出优势。这启发我们，可以在迭代初期使用LU-SGS而在后期使用DPLR方法，这样既能利用LU-SGS方法单步计算量小以及稳定性好的特点，在时间推进的初始阶段更好地避免强烈非定常特性给收敛带来的不利影响并使残差快速下降，又能发挥DPLR格式在迭代后期的优势。

为了考察多块网格下DPLR方法的性能。本文把二维球头绕流网格和三维球头绕流网格进行分块。有两种分块方式，一是在流向和周向分块；二是在法向分块。

图 8是对初始二维球头绕流网格在流向进行等间距分块、得到的不同网格分块数流场的残差收敛曲线。由图可见，不同网格分块数流场的残差收敛曲线基本重合，在流向进行网格分块对DPLR算法的收敛性能几乎没有影响。图 9是对三维球头绕流网格在流向和周向进行等间距分块、得到的不同网格分块数流场的残差收敛曲线，由图可以得到和二维球头绕流同样的结论。这很容易理解，参看式(9)，在第m松弛步，DPLR算法在每条法向网格线上求解，而流向网格上的δQ已经在m-1松弛步得到计算，这就是说每条法向网格线上的隐式求解都是独立的、不会相互影响，因此在流向和周向进行网格分块不会影响算法的收敛过程。

下面考察在法向进行网格分块对算法性能的影响。图 10是对二维球头绕流网格在法向进行分块的分块位置示意图，本文分别考察在两个分块位置分块后形成的多块网格流场的收敛性能，其中分块位置1(splitting 1)远离边界层，分块位置2(splitting 2)位于边界层内。图 11展示了二维球头绕流网格在法向分块后得到的多块网格流场收敛曲线，由图可见，在分块位置1进行网格分块对收敛性能影响较小，但在分块位置2进行网格分块降低了算法的收敛性能。分析原因，相较于分块位置1，分块位置2在边界层内，流场法向梯度大，相邻法向网格单元之间流场变量耦合紧密，而DPLR算法在分块后的两个块上分别求解使得两个网格块流场变量联系得不那么紧密，收敛速度下降是必然的。需要说明的是，DPLR隐格式不是无条件稳定的^[8]，在初始网格求解时CFL数可以取到3000(图 6，图 8)，在分块位置2分块后CFL数最大只能取到1800(图 11)，这说明分块后格式稳定性下降。图 11中网格分块后流场的DPLR计算，采用本文提出的隐式对接边界条件(见式(13)，图中标注为“implicit interface boundary condition”)，而不采用隐式对接边界条件即直接舍弃与虚拟网格点有关的项，算法CFL数最大可取到100，收敛很慢，在图中未给出。

图 12是三维球头绕流网格在法向分块后得到的多块网格流场收敛曲线，图中的两个分块位置(splitting 1和splitting 2)与二维球头绕流网格法向分块位置类似，即分块位置1远离边界层，分块位置2在边界层内。图 12中在分块位置2分块得到的多块网格的流场在计算时CFL数取不到初始网格和分块位置1计算时的3500，最大可取到2000，说明三维计算时在法向进行网格分块导致格式稳定性下降。从是否采用隐式对接边界条件的收敛曲线对比可以看出，采用隐式对接边界条件能够改善DPLR算法在法向分块后收敛速率下降的问题。

3 返回舱绕流计算分析

返回舱外形复杂，包含凸起物、耳片等局部结构。计算网格分区如图 13所示，共77个分区，约116万网格数。在这样的网格上计算，DPLR算法将会出现相邻网格求解方向交叉的现象，这是与上文球头绕流计算的不同之处，也是对本文提出的隐式对接边界条件模型的重要检验。为了与风洞试验结果对比，本文选择风洞试验状态^[19]作为计算状态：模型缩比6%，马赫数9.96，基于模型长度的雷诺数2.04×10⁵，总温1005 K，总压1.976 MPa，迎角10°~30°。风洞试验在FD-17A Φ1 m高超声速低密度风洞^[20]进行，这是一座是高压下吹-真空抽吸暂冲式风洞。分析试验条件，自由来流努森数为7.08×10^-5，流动处于连续流范围；雷诺数较低，流动保持层流的可能性大；总温较低，流动的真实气体效应很弱。因而在计算中，本文选择层流和完全气体模型。

图 14给出迎角22°状态的残差收敛曲线，图中包含了DPLR和LU-SGS两种方法的对比，其中DPLR算法的CFL数为5000，而LU-SGS算法的CFL数取到了1.0×10⁶，图中的横坐标是CPU时间。两种隐式算法都基于原始网格分区在Intel(R) Xeon(R) CPU E5-2680共8核上进行并行计算，未考虑负载均衡。由图可见，LU-SGS算法残差最多只能下降3个量级，而DPLR算法的残差可以下降到机器零。这不仅说明了DPLR算法和本文提出的隐式对接边界条件能够成功地应用于复杂外形，也再次说明了DPLR相对于LU-SGS在收敛速度上的优势。

图 15给出返回舱在迎角22°对称面和物面的压力、马赫数等值线云图以及对称面的流线结构，可以看到返回舱头部的弓形激波和表面凸起物附近模拟得到的流场结构较为清晰。从图 15(b)的对称面流线结构可以看到，在返回舱背部的凸起物前部，由于凸起物的遮挡，流动被局部滞止，马赫数减小并发生了流动分离；而在凸起物的后部，流动分离形成旋涡，部分分离的流线折向返回舱底部并与下部流线汇集。图 16给出返回舱的阻力、升力和俯仰力矩系数随迎角变化关系。可看出，阻力和升力系数随迎角的增大而减小；俯仰力矩系数随迎角的增大而减小，说明返回舱具有纵向静稳定性。计算和试验符合较好，阻力系数相对偏差不超过为3%，升力系数相对偏差不超过2%；计算与试验预测的配平迎角分别为22.05°和21.63°，相对偏差为1.9%。图 15、图 16较好地揭示了返回舱绕流面貌。

4 结束语

由于DPLR算法需要和隐式边界条件相结合，为把DPLR应用到返回舱再入近空间连续流区气动力/热复杂绕流多块对接结构网格的计算，本文首先提出一种隐式对接边界条件的处理方法，通过在对接边界处应用隐式边界条件把边界条件对隐式求解的影响考虑进来。然后研究了网格分块方式对算法稳定性和收敛性能的影响。最后进行了返回舱绕流的数值模拟，给出气动力特性，并与风洞试验结果对比，考察算法在复杂外形数值模拟中的性能。可得出如下结论：

(1) DPLR算法采用松弛方法把求解方向(一般是物面法向)和其它方向解耦，在非求解方向上进行网格分块不会对算法的收敛性能产生影响。而在求解方向上对网格分块，特别是分块位置在边界层内，会降低算法的稳定性和收敛性能。

(2) 本文提出的隐式对接边界条件处理方法，能够改善在求解方向进行网格分块时造成的算法稳定性和收敛性能下降的问题。

(3) 如果应用多块结构网格和DPLR算法进行数值模拟，在用网格软件生成网格或基于并行计算的目的剖分网格时，本文的研究给出一些注意事项：可以在流向和周向进行网格分区，尽量不要在物面法向进行网格分区，更不要在边界层内沿物面法向进行网格分区。这为相关方面的研究提供指导原则。

(4) 返回舱连续流区的复杂绕流计算结果表明，DPLR算法结合隐式对接边界条件能够成功地应用于复杂外形结构网格数值模拟，并以较快的速度收敛，大大缩短了计算所需时间。

本文的研究是初步的阶段性工作下，一步将把DPLR算法应用到湍流以及热化学非平衡流动的数值模拟。