傅 林，高正红，张晓东

（西北工业大学航空学院 翼型、叶栅空气动力学国防科技重点实验室，陕西 西安 710072）

0 引 言

Godunov首次将Riemann问题及其相关解引入到计算流体力学领域，提出了Godunov方法，将问题转化为求解网格单元界面的Godunov型通量。目前比较成熟和应用广泛的为一阶通量数值计算方法，具有较强的数值稳定性并且可以正确的收敛到熵解。主要有基于通量差分思想和基于矢通量分裂思想的两类格式，包括著名的ROE，Van Leer等格式。ROE格式在跨声速和超声速计算中经常由于不满足熵条件，而要引入熵修正，增加了数值耗散；Van Leer格式本身耗散较大，适合超声速流场计算。Harten，Lax，Van Leer等人首先提出了一种求解近似黎曼解的方法－HLL格式［1］。采用HLL格式可以直接求解近似黎曼通量，计算效率较高且稳定，适应范围广。然而由于HLL格式在间断面具有较大的数值耗散，计算结果容易将接触间断抹平，在跨声速复杂流场计算中误差较大。Toro等人提出了一种改进的方法，称为HLLC格式［2］。HLLC近似黎曼求解器能够准确的捕捉激波、接触间断和稀疏波。与精确求解黎曼问题的方法相比具有计算量小、低耗散、高分辨率的特点，且严格满足熵条件。在低速及跨声速可压缩流场中，得到了成功且广泛的应用。

经数值研究发现在可压缩流中，当马赫数较低时，HLLC格式可以获得比ROE相对于实验更好的结果［1］。然而随着马赫数继续增加，HLLC格式在强激波区域会出现激波不稳定的现象［3］。相反，HLL格式不会出现激波不稳定现象，但是由于HLL格式的耗散过大，得到的解略差于ROE格式。本文研究了HLL，HLLC格式的数学原理，构建了混合的HLL－HLLC格式。其在低马赫数下仍然能获得较低的耗散，在强激波流场中，可以克服HLLC格式的激波不稳定现象。为此，优先设计了激波探测器，用来准确捕捉流场中的强激波区域。在强激波区域附近，可实现耗散较大的具有HLL性质的格式，以克服激波不稳定现象；而在远离强激波区域采用低耗散的具有HLLC性质的格式。从而使HLL－HLLC格式在各种马赫数条件下都能得到较好的解。

文中计算采用的求解器为作者自行研究开发的具有定常、非定常数值模拟能力的三维有限体积代码－LMNS3D LMNS3D拥有基于CPU计算的MPI并行技术和基于GPU计算的CUDA并行技术，能够进行大规模并行计算。通量计算可选择不同格式，如ROE，HLL，HLLC［4－6］，采用 MUSCL插值获得二阶空间离散精度，同时具有四种通量限制器可供选择。湍流模型为在航空界应用广泛的SA一方程模型和K－Omega SST两方程模型。粘性项采用二阶中心差分离散，拥有薄层近似及全NS计算能力。时间推进技术采用LU－SGS隐式方法，同时采用低速预处理技术进行低速流场数值模拟。

1 数值方法

1.1 控制方程

在物理坐标系下NS方程为：

经过雅克比转换，NS方程转化为计算系下的守恒形式，为本文求解的模型方程：Q′代表有限体体心值，等代表界面上的无粘数值通量，由一阶单调近似黎曼求解器给出；R′、S′、T′代表粘性数值通量。要得到高精度格式，要求提高界面上守恒变量值的插值精度。

1.2 HLL及HLLC近似黎曼求解器

HLLC近似黎曼求解器基于HLL发展而来，HLLC格式将黎曼问题求解区域由向左传最快波系、接触间断、向右传最快波系分成四个部分（如图1所示）。最左边和最右边的波系可能为激波或稀疏波，用符号表示。经过这些波系，流场变量发生变化，特别是激波使流场解变得不连续。中间的波系为接触间断，左右两边的流场变量为。求解时假设在各自的扇形区域里是常量。Harten等人通过对每一个扇形域中守恒变量做积分平均，提出了一种方法来简化黎曼问题解的结构。Harten、Lax、Van Leer等人首先提出了一种近似黎曼求解器HLL［1］，忽略了中间接触间断的影响，即认为接触间断两边的是相同的。理论研究表明［1］：如果解收敛，那么可以收敛到满足熵条件的物理解。HLL格式是一种非常有效和鲁棒的格式，可以直接应用于Euler／NS方程。

图1 黎曼问题示意图Fig.1 Riemann problem schematic

其解向量的表达式为：

然而HLL求解器具有容易抹平接触间断的特点。为了克服这一弱点，Toro等人提出了HLLC（C表示contact）格式来精确捕捉接触间断［2］。其解的结构为：

HLLC格式需要预先估计最慢波速、中间波系波速和最快波速。本文采用基于ROE平均的直接波速估计方法，具有很好的鲁棒性和精度［8］。

则HLLC格式的通量为：

HLLC格式具有以下优点［9］：

（1）精确捕捉激波和接触间断。

（2）保证计算过程中的标量衡正。

（3）增强的熵条件，不需要熵修正，适合跨声速和低速计算。

（4）不需要知道通量雅克比矩阵的细节。

1.3 HLL－HLLC格式的构造

1.3.1 激波不稳定现象的产生机理

Liou［10］分析了不同数值通量计算格式的耗散特性。在分析中，把每种格式的耗散项展开，用主变量ρ、u和p的差分形式表达，定义密度差和压力差的系数为算子D（ρ）、D（p）。研究不同格式（AUSM ＋，AUSMDV，HLLE和ROE）耗散项中的。研究结果表明：激波不稳定现象出现的根源是在求解质量通量的算法中包含了D（p）的影响；当一个格式求解质量通量算法中，具有D（p）＝0的特性，则是一个激波稳定的格式。Pandolfi和D′Ambrosio通过进一步研究发现，当一种格式显式的处理了接触间断，就会显示出明显的激波不稳定现象［11］。从已有的研究中可以看出，HLLC格式之所以在强激波区域出现激波不稳定现象，是因为显式的表达了接触间断。

具体分析HLLC格式发现：连续性方程并不显式的包含压力差的贡献，但是由公式（19）预估波速的时候，求解中间波系波速的公式中包含压力差项的影响。

所以HLLC格式无法避免出现激波不稳定现象。而同样的分析可以发现，由于HLL格式并没有显式的表达中间波系，所以不会出现激波不稳定现象。

1.3.2 克服HLLC格式激波不稳定现象的方法

本文试图在HLLC的框架下调整通过黎曼解中间波系（接触间断）的通量差分来消除HLLC的激波不稳定现象。减小通过中间波系的解向量的变化，三波系的黎曼问题会退化为两波系黎曼问题，因此新格式将具有HLL格式的性质（如图2所示）。

图2 状态解构造示意图Fig.2 Construction schematic of state solution

f为本文定义的激波捕捉函数［12］：

其中pL和pR是作用在网格面上左右两边的压力。在激波以外的区域，压力连续的变化，所以混合函数f的值为1.0，此时格式具体表现为HLLC的性质；在强激波附近，压力差pR－pL很大，混合函数f为一个小于1的值，格式表现为HLLC和HLL的混合型；当f很小时，格式表现为完全的HLL性质，从而抑制激波不稳定现象的产生。在三维应用中，三个方向的混合函数f的定义为：

2 HLL－HLLC格式在可压缩速流动中的应用

2.1 高超声速双楔流动数值模拟

Neuenhahm等人为了研究前缘半径Rn大小和壁面温度对高超声速双楔流动的影响做了一系列实验［13］。超声速双楔流动包含强激波附面层干扰和小分离现象，是考察数值计算格式和求解器的典型算例。本算例计算模型如图3。

图3 高超声速双楔流动模型Fig.3 Hypersonic flow on the model of double wedge

计算条件为：来流马赫数为8.1，静压为520Pa，静温为106K，单位长度雷诺数为3.8×106。壁面为等温无滑移边界条件，Tw为300K。如图4所示，计算网格为三块点对接网格，y＋的最大值为0.9。无粘通量分别采用四种数值格式－ROE（不加熵修正）、HLL、HLLC、HLL－HLLC。计算采用层流模型。

图4 超声速双楔模型计算网格Fig.4 Computational grid of double wedge

图5为四种格式计算压力云图的比较，可以看出HLLC格式在这种强激波流动中出现了激波不稳定现象，激波后出现非物理解。ROE、HLL和本文构造的HLL－HLLC格式计算结果较好。图6为不同格式计算压力系数分布与实验对比，由于实验洞壁干扰等因素，实验中取了三次实验值分别为EXP＿1、EXP＿2、EXP＿3。HLLC格式出现激波不稳定现象，压力分布出现非物理震荡；同时可以看出，HLL格式耗散过大，分离点靠后与实验结果差距较大。本文构造的格式HLL－HLLC，计算稳定，结果优于HLL格式，接近于ROE格式。从图7可以看出，文中设计的激波感受因子能精确捕捉强激波。图8给出了HLL－HLLC格式计算的双楔交界处的小分离区域流线。

图5 压力系数分布云图Fig.5 Comparison of Cp distribution

图6 压力系数分布与实验比较Fig.6 Comparison of aerodynamic coefficient

图7 激波感受因子捕捉到的激波示意图Fig.7 Shock captured by shock feeling factor

图8 小分离区域流线示意图Fig.8 Streamline schematic of small separation region

2.2 超声速后台阶流动数值模拟

超声速后台阶流动是典型的复杂流动问题（如图9），流动现象包括：湍流边界层、膨胀波、流动分离和再附着、再附激波等，从而成为考察CFD求解器和计算格式的典型算例之一。Smith［14］等人对此问题进行了一系列实验，分别研究了不同台阶高度（0in、0.433in、0.75in），马赫数从2.5到5.0，雷诺数从2.5×105到1.8×106等温壁面下的流动情况。

图9 超声速后台阶流动模型Fig.9 Supersonic flow over the model of backward－facing step

本文计算条件为：来流马赫数为2.5，静压为1.5×104Pa，静温为170K，单位长度雷诺数为1.8×107。如图10，台阶高度为0.433in，计算网格为三块对接网格，y 的最大值为1.25 算例假设壁面条件为无滑移绝热壁面，计算区域的上边界为滑动壁面。分别采用ROE格式和本文发展的HLL－HLLC格式计算无粘通量，湍流模型采用SA一方程湍流模式。

图10 超声速后台阶流动计算网格Fig.10 Computational grid of backward－facing step

图10中，x轴和y轴均采用台阶高度做归一化处理，按照实验中的方法，压力采用x／H＝－3.39处的压力值进行归一化处理。图11显示，两种格式计算的空间压力分布云图一致。从图12马赫数云图可以看出，Prandtl－Meyer膨胀波、再附着激波和分离线非常明显。图13表明，本文两种格式计算的再附点位置均为x／H＝2.5512附近，文献［15］采用SA湍流模型计算的分离点位置为x／H＝2.53。从图14可以看出，壁面压力分布和实验值吻合较好。其中HLL－HLLC格式和不加熵修正的ROE格式计算结果一致，均表现出较高的分辨率；在分离区的压力分布好于加熵修正的ROE格式。从本算例可以看出，在中等超声速问题中本文发展的HLL－HLLC格式，和ROE格式具有相似的耗散性和精度。

图11 压力分布云图比较Fig.11 Comparison of Cp distribution

图12 马赫数云图比较Fig.12 Comparison of Ma distribution

图13 流线图比较Fig.13 Comparison of streamline

图14 压力分布与实验比较Fig.14 Comparison of aerodynamic coefficient

2.3 高超声速钝头体数值模拟

计算条件为：来流马赫数10，单位长度雷诺数1×105，壁面为等温壁面274K。模型为半径1m的半球。分别采用 HLLC，HLL－HLLC、加熵修正的ROE和AUSM＋四种格式计算，计算模式为层流。

从图15可以看出，HLLC格式计算结果出现了激波不稳定现象，捕捉到的激波为非物理解；另外三种格式HLL－HLLC、加熵修正的ROE和AUSM＋均正确反映了钝头体前的弓形激波。其中加熵修正的ROE格式捕获的激波厚度最厚，HLL－HLLC和AUSM＋格式计算得到的激波厚度基本相当；HLLHLLC相对于AUSM＋格式得到的激波后压力分布更加光滑。另，作者尝试在本算例中使用未加熵修正的ROE格式，但计算失败。由此可以看出，本文发展的HLL－HLLC近似黎曼求解器具有较低的数值耗散性和较高的精度。

图15 压力分布云图比较Fig.15 Comparison of Cp distribution

3 结 论

（1）近似黎曼求解器是否出现激波不稳定现象强烈依赖于其是否显式捕捉接触间断。HLLC格式精确的显式表达接触间断，在高超声速流场中，强激波附近容易出现激波不稳定现象；HLL格式由于抹平接触间断不会出现这种现象。本文发展的HLLHLLC格式有效的解决了HLLC格式的激波不稳定现象，同时相对于HLL格式具有更低的耗散性和更高的精度。

（2）本文设计的激波捕捉函数能够准确捕捉强激波。在强激波附近，HLL－HLLC表现为 HLL性质的格式，以克服激波不稳定现象；而在远离强激波区域过渡到HLLC性质的格式。通过三个算例证明了本文所构造格式的有效性和鲁棒性。

［1］TORO E F.Riemann solvers and numerical methods for fluid dynamics［M］.Springer：Berlin，1999.

［2］TORO E F，SPRUCE M，SPEARES W.Restoration of the contact surface in the HLL－Riemann solver［J］.Shock Waves，1994，4（1）：25－34.

［3］PANDOLFI M，D’AMBROSIO D.Numerical instabilities in upwind methods：analysis and cures for the‘Carbuncle’phenomenon［J］.Journal of Computational Physics，2001，166（2）：271－301.

［4］NIU Y Y，CHANG C H，TSENG W Y，et al.Numerical simulation of an aortic flow based on a HLLC type incompressible flow solver［J］.Communications In Computational Physics，2009，5（1）：142－162.

［5］HARTEN A，LAX P D，LEER B V.On upstream difference and Godunov－type schemes for hyperbolic conservation laws［J］.SIAM Review，1983，25（1）：35－61.

［6］KITAMURA K，ROE P，ISMAIL F.Evaluation of Euler fluxes for hypersonic heating computations［J］.AIAA Journal，2009，47（1）：44－52.

［7］REMAKI L，XIE Z Q，HASSAN O，et al.A high order finite volume－HLLC solver and anisotropic Delaunay mesh adaptation［R］.AIAA 2009－1498，2009.

［8］BATTEN P，CLARKE N，LAMBERT C，et al.On the choice of wave speeds for the HLLC Riemann solver［J］.SIAM Journal on Scientific Computing，1997，18（6）：1553－1570.

［9］NICHOLS R H，TRAMEL R W，BUNING P G.Solver and turbulence model upgrades to OVERFLOW 2 for unsteady and high－speed applications［R］.AIAA－2006－2824，2006.

［10］LIOU M S.Mass flux schemes and connection to shock instability［J］.Journal of Computational Physics，2000，160（2）：623－648.

［11］KIM S S，KIM C G，RHO O H，et al.Cure for the shock instability：development of a shock－stable Roe scheme［J］.Journal of Computational Physics，2003，185（2）：342－374.

［12］QUIRK J J.A contribution to the great Riemann solver debate［J］.International Journal for Numerical Methods in Fluids，1994，18（6）：555－574.

［13］REINARTZ B U，BALLMANN J，BOYCE R.Numerical investigation of wall temperature and entropy layer effects on double wedge shock／boundary layer interactions［R］.AIAA 2006－8137，2006.

［14］SMITH H E.The flow field and heat transfer downstream of rearward facing step in supersonic flow［R］.ARL－67－0056，1967.

［15］DIPPOLD V，MOHLER S R.Validation of the wind－US unstructured flow solver for wall－bounded flows［R］.AIAA 2006－637，2006.