面向三维激波问题的装配方法

邹东阳，林敬周，黄洁，刘君

1. 中国空气动力研究与发展中心 超高速空气动力研究所，绵阳 621000 2. 大连理工大学 航空航天学院，大连 116024

强扰动在空气中以超过声速的速度进行传播，表现为流体中出现一个参数的间断面，称为激波[1]。激波间断的存在给描述无黏流动的Euler方程的求解带来了很大挑战。过去半个多世纪里，为了发展出能够精确地描述激波的数值方法，科研人员做了大量的工作。如今这些工作是计算流体力学(Computational Fluid Dynamics, CFD)的重要组成部分。根据对激波处理方式的不同，CFD可以分为激波捕捉方法[2-8]和激波装配方法[9-11]2种。激波捕捉方法是以弱解理论为基础，全场采用统一的计算方法，无需特殊处理，具有通用性的特点。激波装配方法则是利用激波关系式(Rankine-Hugoniot，R-H)对激波进行描述，理论上是精确的，但是需要找到匹配激波间断解与其他光滑区域数值解的方法，存在兼容性的问题。对于复杂激波流动问题，特别是三维问题，激波装配方法存在应用上的困难。随着激波捕捉方法的成熟，特别是TVD(Total Variation Diminishing)格式的诞生，对于激波的计算逐步满足了一定的要求。过去30多年里，激波捕捉方法得到了迅猛发展，在航空航天、汽车建筑、天文地理等领域发挥了巨大的作用。

虽然激波装配方法的研究在相当长一段时间内寥若晨星，但是其在计算激波时表现出来的优势还是给许多激波捕捉方法研究提供了思路。像激波装配方法那样，将计算网格面与激波面相匹配来减少激波计算误差的处理手段在很多方法中都得到了应用。在20世纪90年代中期，Parpia和Parikh[12]、van Rosendale[13]以及Camarero[14]等发展了不同的方法，全都是将非结构自适应三角形网格的棱边与流场中出现的间断进行匹配。类似的思想在新近发展的间断伽辽金有限元方法(DG-FEM)中也有一定的体现。两个完全独立的研究团队同时关注于解决网格节点(不包括流场参数)与间断面相匹配的问题[15-20]。

另外一方面，在过去的十几年里，CFD工作者对于激波装配技术的兴趣得到了急剧的提升，将目光从激波捕捉方法重新转向激波装配方法。造成这一现象的原因主要有2个方面：一是由于CFD工作者在解决激波捕捉方法本质性局限所遇到的困境；二是由于一些研究人员将激波装配与非结构网格相结合取得了显著成果，给解决激波捕捉方法所遇到的问题带来了希望。

意大利学者Paciorri和Bonfiglioli自2006年开始着力于研究基于浸没边界法的激波装配技术。他们方法有效地结合了浮动激波装配和边界激波装配的特点，激波的运动过程中允许其在背景网格上进行(该特点类似浮动激波装配)，计算过程中通过挖洞构造以激波为内边界的新的计算网格进行求解(该特点类似边界激波装配方法)。Paciorri和Bonfiglioli的方法较传统的激波装配方法在灵活性上有了显著的提高，利用发展的装配方法成功地解决了若干含有激波构型的问题[21-23]。在国内，刘君等基于非结构动网格技术对激波装配方法进行了发展，提出了新型非结构激波装配方法—嵌入式激波装配方法。嵌入式激波装配方法中，激波属于网格的一部分，激波节点的移动带动其他网格节点的运动，灵活性有了很大提高。该方法被成功用于求解许多定常/非定常流动问题[24-29]。

基于非结构网格的装配方法在诸多二维问题中的成功应用展示了其相比传统激波装配方法在灵活性上的优势。相比二维问题，三维问题虽然只是增加了一个维度，装配难度却是有了大幅度提升。三维空间的激波辨识难度大，给初始激波位置的选择带来了困难。描述激波面的网格点在运动过程中彼此要相互匹配，以保证激波收敛运动过程不发散。在三维空间中控制这些点的匹配关系要比二维更为复杂，对算法要求也较高。此外，三维激波相交等问题缺乏理论指导，使得装配方法在复杂问题中的应用难度较大。诸多问题的存在使得针对三维问题的装配方法应用研究还较少，目前只是在文献[30]中看到了相关报道。Bonfiglioli等[30]应用其发展的激波装配方法研究了一些超声速流动问题，证实了其方法在解决三维问题中的有效性。尽管如此，从文献给出的例子还是可以看出，装配的激波相对较为简单，流场内激波都是相互独立的，规避了三维激波相交带来的问题。本文基于前期关于嵌入式激波装配方法的研究工作，进一步从二维拓展应用到对三维问题的模拟中。与Bonfiglioli等[30]给出的三维装配不同，本文在对三维激波问题进行装配时，不只是针对独立激波，还考虑了更为复杂并且具有代表性的三维激波相交、反射等问题，结合给出的激波交点的计算方法，完成了对三类超声速/高超声速三维激波问题的装配。

1 三维激波装配方法

激波装配方法由2个要素组成：计算和追踪。所谓的“追踪”指的是对流场中出现各类间断位置进行追踪。所谓的“计算”指的是对间断边界包络计算域内流动利用CFD求解器对偏微分控制方程进行求解的过程。

无论使用什么样的CFD求解器，在非结构网格范围内装配激波是由一系列激波面所表示，这些激波面是由激波点自动连接而成。装配间断的位置是通过在激波点上求解R-H关系式获得。

根据求解器是格心型或者格点型的不同，追踪激波的方式也会有很大的差异，需要相应地去修改CFD求解器。本文的激波装配方使用的是格心型有限体积求解器。追踪激波时，在激波点上求解R-H关系式，以此来获得激波运动速度。通过激波面元的通量则是使用上游参数以及激波速度得到。一旦获得通过激波面元的通量，包含激波面元的单元参数的更新方式就和其他单元一样了。装配的激波按照求解的激波点速度在拉格朗日体系下移动。采用三维非结构激波装配方法进行数值分析的过程涉及4个主要内容，下面分别进行详细介绍。

1) 激波间断标识

激波装配方法不像激波捕捉方法那样对全场采用同样的处理方式，而是将激波间断单独标识出来，用一种解析的有别于光滑区域的方法进行求解。传统上的边界激波装配方法是将激波看作一种特殊的边界进行标识和计算的。这种方法需要将计算区域划分成若干个子区域，对于含有悬空激波以及拓扑结构变化的问题不易操作。本文中使用的激波装配方法是通过节点属性定义来对激波进行标识的。

如图1所示，采用网格生成程序(Tetgen)或软件(Pointwise)生成空间四面体网格。对离散后的网格节点进行属性定义，划分为激波节点(Shock nodes, S)和普通网格节点(Common nodes, C)2种。激波面元是由激波节点连接构成的，规定如果一个网格面元包含的节点全部为激波节点则该面元为激波面元。在计算过程中可以根据流场变化对网格节点属性进行调整，这就极大地提高了激波装配的灵活性。

图1 激波节点标识Fig.1 Labels of shock nodes

2) 网格节点参数确定

采用格心型有限体积方法时，流场参数是储存在单元格心的。网格节点参数V=[ρ,ux,uy,uz,p]T需要利用格心的流场参数U=[ρ,ux,uy,uz,p]T通过重构获得。其中，ρ、p分别表示流场密度和压力；ux、uy、uz为对应于3个方向的流动速度。按照步骤1)中的分类，对于激波节点和普通节点2种节点类型需要2种不同的计算方法。

普通网格节点上的参数主要通过格心参数加权平均得到

(1)

式中：Vi表示节点i上流场参数；N为与网格节点i相邻的单元数；αik和Uik分别为与之相邻的第k个单元影响的权系数和格心参数。αik采用节点与单元距离|Rik|的反比进行确定：

(2)

式中：xi和xk为节点i和与之相邻的第k个单元的格心坐标;Rik为第k个单元与节点i的距离矢量。

基于单元格心参数，利用式(1)、式(2)确定节点参数的方式可以看出，对于普通节点而言，其影响域包含了四周全部的单元。

使用装配方法对激波间断处理时需要用到上下游的两组参数。这个特点决定了被标识的激波节点与普通节点的一个重要区别：包含2组参数，一组代表上游低压区的影响，一组代表下游高压区的影响，如图2所示。对应于上下游，确定其激波节点2组参数所使用的网格单元或影响域和普通节点有所不同，上游节点参数Vu只使用上游低压区单元，下游节点参数Vd只使用下游高压区单元

(3)

式中：下标“u”“d”分别表示激波上下游；N1、N2为以激波面为分界包含激波节点的上下游单元个数。与式(1)中N不同的是，在求解式(3)时并非将所有与激波节点相邻的单元都纳入考虑范围，而是要求影响单元需要包含激波面元。此外，加权系数αik的确定也有所不同，采用式(4)确定:

图2 激波节点参数确定Fig.2 Determination of states on shock nodes

(4)

(5)

式中：mτk为第k个单元沿所包含激波面切向的马赫数，即

(6)

式中：定义u=[ux,uy,uz]T为流场速度矢量，则uk为第k单元的格心流动速度；xc为其所包含的激波面元运动速度；R′为对应激波面元中心到激波节点的空间位置矢量；定义a为声速，则ak为单元k个单元的格心声速。

利用光滑区域内单元格心参数确定激波节点上下游参数的过程建立了解析解和数值解之间的联系。对于激波间断而言，利用激波节点上下游参数通过求解R-H关系式可以得到激波节点的运动速度。由于R-H关系式是一维的，所以三维激波装配还涉及到三维问题和一维方程之间的匹配问题。通过建立基于激波节点法向的激波坐标系，可以将三维问题降为一维进行求解。

第i个激波节点法向ni是通过与之相邻的激波面元法向确定：

(7)

式中：Ns为与激波节点i相邻的激波面元的个数；nj为第j个激波面元的法向，并且定义激波方向由上游指向下游为正，加权系数α′ij的确定与式(4)类似：

(8)

式中：|R′ij|表示第j个激波面元中心到激波节点i之间的距离大小。

设定激波节点运动速度大小为w，方向为n。在激波坐标系下，激波节点上游流动的相对马赫数为

(9)

根据上游相对马赫数，利用R-H关系式并补充反向传播黎曼变量保持不变的条件可以确定下游参数修正值V′d以及激波运动速度

(10)

式中:前4式为补充了切向速度不变条件的R-H关系式，第5式是根据反向传播黎曼变量保持不变条件给出，γ表示比热比。求解上述非线性方程组可以得到激波节点运动速度w以及下游参数修正值V′d。一般情况下，V′d与Vd不相同，修正值V′d主要用于下游单元格心梯度的求解。

3) 激波动态跟踪

本文使用基于弹簧模型的非结构动网格技术实现对激波的显式跟踪。根据式(10)求解的激波节点速度w以及运动方向n，根据t时刻激波激波节点位置xt计算t+1时刻的激波节点位置

xt+1=xt+wnΔt

(11)

式中：Δt为时间步长。根据激波节点位移，利用动网格技术移动其他可移动的网格节点，以保证激波运动变形过程中计算网格的质量。图3给出了从初始时刻到收敛时刻激波面的运动示意图。运动过程中，网格质量较差时还需要配合网格重构的策略。

(4) 激波通量计算

采用格心型有限体积求解器对控制方程进行求解时需要计算流经单元面元的通量值。对于普通面元，可以采用各种成熟的迎风格式或中心差分格式进行计算。而对于激波面元来说，通量计算反而变得更加简单和清晰。在激波坐标系下，激波面元上游为超声速流动，按照超声速流动特征，面元通量Φ只和上游流动有关，因此可以根据式(12)进行确定：

图3 激波节点运动示意图Fig.3 Sketch of shock node motion

(12)

式中：φ1=ρ(u-xc)·n，φ2=xc·n，E=0.5·u·u+p/[ρ(γ-1)]。n=[nx,ny,nz]T为激波面元法向，激波面元运动速度xc可以根据激波节点运动速度wi和法向ni得到

(13)

式中：N为三维激波面元所包含的激波节点个数。

2 格心型有限体积方法

三维Euler方程的积分形式可以写为

(14)

式中：Q为守恒变量；Fc(Q)为无粘对流通量；Ω和∂Ω分别为控制体单元及其外边界。

在选定的控制体上，将式(14)中的对流项进行离散可以得到

(15)

式中：Fck、nk、Sk分别为控制单元的第k个表面边界的通量、外法向单位矢量和面积；Nf为控制体的表面个数。

采用格心型有限体积方法求解上述方程的过程涉及到空间物理量的重构。空间二阶格式通常假设物理量在控制单元内为线性分布，即控制单元内梯度为常数。控制单元表面物理量采用不同单元格心值重构得到的结果在表面边界处形成间断，采用van Leer等格式进行通量计算，可以获得通过单元表面的对流通量值，对各个控制单元边界遍历求和从而获得单元内参数变化量。然后对时间导数项进行离散，就可获得下一时间步的格心处流动参数。关于求解器的详细介绍可以参考文献[31]，在此不再赘述。

3 三维算例验证

3.1 三维弓形激波

高超声速球柱体绕流是考核计算方法对强间断问题计算能力的一个重要算例。本文首先使用三维激波装配方法对来流马赫数10条件下的该问题进行分析。球柱模型设计如图4所示，球头半径为L，柱体与圆柱相接，长度也为L。计算域选为6L×3L×6L的平行六面体。

不考虑攻角和侧滑角情况下，球头绕流为轴对称的，流动结构相对简单，有很多文献和数据可以用来对本文计算结果进行考核。

根据研究内容需要，设计了3套用于球柱体计算的网格，分别为细网格C-1、粗网格F-1和细网格F-2，网格节点和单元数如表1所示。细网格C-1使用均匀划分的方式，首先对物面和计算域外边界使用三角形进行离散，控制参数设为Δ=0.05L，根据边界网格生成内部计算网格，计算域内全部采用四面体进行划分。细网格C-1用于激波捕捉方法的计算。为了对激波点进行标记定义，在激波装配计算网格F-1和F-2中预先设置一个空间面，对该面进行离散后作为内部边界控制区域网格生成，并在计算中将面上的节点定义为激波点。采用与C-1同样的网格生成方式，生成粗网格F-1和细网格F-2。其中，粗网格F-1中物面和激波面上的网格离散控制参数设为Δ=0.1L，外边界控制参数设为Δ=0.3L；细网格F-2中对物面以及激波层内进行加密，网格离散控制参数设为Δ=0.05L，外边界控制参数同样设为Δ=0.3L。用于激波捕捉计算的网格C-1在计算过程中不发生网格变形，而网格F-1和F-2会随着激波的收敛过程发生变形重构。下面对使用3套网格和2种方法生成的计算结果进行对比分析。

图4 球柱体: 计算区域及网格Fig.4 Hemisphere-cylinder: computational domain and mesh

1) 激波捕捉和激波装配对比

本文对激波装配和激波捕捉2种结果对比时都是采用同样的求解器和同样的输入参数进行计算，不同之处在于激波装配方法中存在激波节点的定义，网格允许运动。

图5给出了对称面上的压力分布的计算结果。其中上侧为采用装配方法得到的结果，下侧为采用捕捉方法计算的结果，黑色线条为文献[30]中给出的压力分布。从对比结果可以明显看出，采用装配方法得到的结果和文献符合得非常好，压力分布完全一致；而采用捕捉方法的计算结果存在一个较厚的激波带，激波带对激波层内压力分布存一定影响，得到的压力分布明显前移。进一步，图6中对比了2种方法计算结果在物面上的压力分布。从对比可以明显看出，激波捕捉方法计算结果存在压力分布前移的现象，对应位置压力略低于激波装配的结果。此外，使用激波捕捉得到的压力分布较不规则，这主要是由于采用非结构网格离散时激波和网格面的斜交程度不同导致环向存在偏差。

表1 球柱体计算网格Table 1 Computational meshes for hemisphere-cylinder

图5 无量纲压力等值线云图Fig.5 Normalized pressure iso-contour levels

图6 球头表面的压力等值线图Fig.6 Pressure iso-contour levels on body

2) 激波装配中粗细网格结果对比

使用激波装配方法分别基于网格F-1和F-2进行了计算。图7给出了对称面上马赫数计算结果的对比。为了更加直观地显示F-1和F-2空间网格的分布情况，图7还给出了相应的背景网格，上侧为粗网格F-1，下侧为细网格F-2。可以看出，在激波层内细网格F-2要比粗网格F-1密得多，但基于2套网格得到的激波位置完全吻合，马赫数分布无较大差异，特别是在球头区域。对于马赫数大于2区域内计算结果的偏差主要是由于粗细网格对于膨胀波波系结构的描述以及流场显示时参数拟合的控制点存在差异造成的。

图7 激波装配马赫数等值线图Fig.7 Mach iso-contours for shock-fitting

对比网格F-1和F-2的结果可以得出结论，激波装配方法计算结果对网格粗细并不敏感，使用较少网格量同样可以获得足够高质量的计算结果。从这点来看，使用激波装配方法模拟时可以适当减少网格量，从而提高计算效率。

3.2 三维激波反射

球柱体绕流只包含一个弓形激波，装配难度相对较低。进一步，我们考虑一个较为复杂的包括一个正规反射的管内激波反射问题。设计模型如图8所示，在圆型通道内放置一个圆锥模型，以圆筒直径L为参考，圆筒半径为R=0.5L。圆锥模型底部到理论尖点的长度为L，圆锥半锥角δ=10°。为了更加贴近实际飞行器模型，圆锥头部选为球形，半径为0.02L。

考虑圆筒内来流马赫数Ma∞=2。超声速气流遇圆锥产生脱体圆锥激波，圆锥激波遇筒壁影响会发生反射。由于流动具有轴对称的特点，所以计算区域只选择1/4，即图8中实线框内区域。计算区域采用四面体单元进行离散，四面体单元数约为35万。

图8 激波反射：计算域和网格Fig.8 Shock reflection: computational domain and mesh

采用三维激波装配方法对入射激波和反射激波进行装配。由于激波在壁面发生反射，涉及到激波壁面交点运动参数的确定。三维激波交点运动参数的确定要比二维复杂得多，其运动速度求解不当会出现不相容问题，导致计算失败。所以激波交点运动参数的确定对于三维装配而言至关重要。

根据正规反射的特点，从几何构型的角度出发认为激波在筒壁上反射点的运动速度主要受入射激波影响，应与入射激波相匹配，即与壁面交点的运动速度wR应由入射激波决定，可以利用式(16)进行确定：

wR=(wI·τw)τw

(16)

式中：wI为按照入射激波确定交点参数得到的激波节点速度;τw为壁面的切向向量。

装配方法得到的激波型面如图9所示。为了更加直观地显示激波面形状将计算结果进行对称，可以看出激波面光滑合理，定性说明装配的正确性。对计算结果进行一步分析，给出了装配激波的角度。根据气体动力学相关知识可知，超声速轴对称绕流形成的附体圆锥激波从每个子午面来看都是斜激波。在马赫数Ma∞=2、半锥角δ=10°条件下，激波角β=31.2°。从图9的对比可以看出，由于圆锥头部为球形，激波会脱体，对下游流动造成影响。随着距离增加，激波角减小，最后收敛至锥形流理论激波角31.2°。入射激波壁面反射为正规反射，反射激波角度可以根据入射激波角31.2°和来流马赫数Ma∞=2进行确定。经计算反射激波角为β′=30.6°，与装配激波角度吻合较好。值得说明的是，远离洞壁时反射激波上游并非匀直流，反射激波角度会出现一定偏差。

图9 激波反射：装配得到的激波Fig.9 Shock reflection: fitted shock waves

3.3 三维激波-激波相交

最后考虑一个更为复杂的含有激波相交的流动问题—球柱锥组合体超声速绕流问题。针对该问题Houtmann[32]、Bonfiglioli[30]分别进行了实验和数值研究。选用他们研究的模型和参数，如图10所示，使用三维激波装配方法进行模拟计算。具体模型尺寸如图10所示，来流马赫数Ma∞=4.04，攻角为α=20°。该问题涉及到球头和拐角处形成的脱体激波和附体激波的装配以及对两个激波形成的三维激波相交的装配，装配难度较大。

根据流动特点，计算区域选择如图11所示。采用完全的四面体网格对计算区域进行离散，离散控制参数设为Δ=2 mm，离散网格量约为94万，用于激波捕捉方法使用的计算网格。如前所述，使用激波装配方法计算时为了对激波点进行标记定义，预先设置激波面，通过激波面和边界对计算区域进行网格划分，网格量约为110万。

图10 模型尺寸图Fig.10 Body geometry and dimension

图11 激波相交：计算域及网格Fig.11 Shock interaction: Computational domain and mesh

在来流条件作用下最外层形成的脱体弓形激波与锥裙和圆柱的结合处形成的内嵌激波相交形成VI类激波相交构型，如图12所示。如前文所述，对于该问题进行装配的难点和重点依然是如何确定激波交点运动速度。根据流动特点，本文给出了一种适用于此类三维激波交点运动速度的确定方法。通过对该类型激波结构的特点进行分析可知，激波交点T的运动速度由处于上游的2个激波I1和I2的运动决定。利用此关系，假设I1和I2的激波运动速度分别为w1和w2，运动方向分别为n1和n2，激波交点T的运动速度w可以通过求解式(17)确定：

(17)

式中：n=n1×n2为平面Π的法向向量。

利用式(17)对交点运动进行计算，得到收敛的激波装配结果，如图13所示。从结果来看，激波型面光滑合理，根据控制方程解唯一性特点可以在一定程度上反映装配结果的正确性。利用三维激波装配方法得到的激波型面使得对三维激波的认识更加直观。

图14给出了激波装配方法和激波捕捉方法在对称面上的马赫数云图对比，从对比结果可以看出：激波装配和激波捕捉得到的流场结果分布基本一致，也反映了装配结果的合理性，从而验证了三维激波装配方法适用于对复杂激波构型问题的装配。进一步对比可以看出，捕捉方法得到的激波存在一个较厚的激波带，随着激波强度的变弱厚度还有增大的趋势，而激波装配得到的激波为理论上的厚度为零。

图12 激波结构示意图Fig.12 Sketch of shock structure

图14 对称面上马赫数云图Fig.14 Mach number iso-contours on symmetric plane

4 结 论

本文将基于非结构动网格的激波装配法拓展应用于对三维问题的研究，数值模拟了若干含有不同激波个数以及不同激波结构的超声速/高超声速流动问题。计算结果显示该方法具有适用于三维问题的能力。通过对高超声速球柱绕流问题的模拟验证了该装配方法在三维问题中的合理性。在激波壁面反射和激波相交问题中，通过对激波结构的几何关系的分析发展了适用于描述这两类三维激波交点运动的方法。理论上，在其他类型的三维激波相交问题中利用激波构型的几何特征确定激波交点速度的方法也是适用的。通过本文工作，对传统装配方法进行了改进，拓展了装配方法的应用范围。三维激波装配可以获得合理的激波型面，增加了对三维激波形状的直观理解，对于相关理论方法的研究也具有一定参考意义。

激波装配方法发展的一个重要目的就是解决其通用性问题，实现其在更为复杂的工程问题中的应用价值。针对三维问题进行研究，会发现二维装配所遇不到的问题，以解决相关问题为导引有利于加速发展出更为通用的激波装配方法。本文关于三维激波装配方法的研究正是沿着这一思路而来。但是从目前来看，关于三维激波装配仍然还有很多工作要做。本文提出的激波交点的计算方法是利用激波结构的几何关系获得，背后的物理含义仍然需要作大量的探索研究。对于基于非结构动网格的激波装配方法，给定的初始激波和收敛时的激波位置距离较远时，激波运动幅度较大，需要不断进行网格重构，严重影响计算效率。精确的激波辨识可以提供相对准确的激波面，有助于提高计算效率。关于激波面的辨识依然是三维激波问题的一个重点。诸如此类的问题在一定程度上都还制约着激波装配方法的实际应用，有待后续深入研究。