包含激波的复杂流动精细化模拟是计算流体力学（CFD）领域的前沿和难点问题。激波捕捉格式是指能够计算包含间断的一系列数值方法。在激波捕捉方法中，针对间断本身没有采用任何特殊处理方法，流体使用守恒形式控制方程进行描述，任何激波与间断都在计算场中一起求解。激波捕捉方法一直是CFD研究的重点和难点，在有限差分和有限体积框架下，研究者们发展了大量经典的激波捕捉方法[1-3]，虽然现代间断捕捉方法已经能够解决激波捕捉问题，但是针对一些复杂流动现象，激波捕捉方法仍然需要进一步发展和完善。作为复杂流动现象的一个典型代表，以激波/湍流干扰为典型特征的可压缩湍流是自然界和工业设计中常见的流动问题。由于传统的基于雷诺平均N-S方程方法无法满足精细模拟的需求，因此开展了直接数值模拟（DNS）和大涡模拟（LES）方法的研究。以DNS和LES为代表的精细化湍流模拟技术对数值离散误差有严格要求，此时需要采用高阶精度格式。可压缩湍流中同时存在激波和湍流并且两者发生动态干扰，这对数值方法提出两类相互矛盾的数值特性要求：既要低耗散模拟湍流小尺度结构，又要有足够耗散抑制激波附近虚假振荡[4]。此外，由于需要降低数值耗散实现复杂流动结构精细化模拟，高阶高分辨率激波捕捉方法更容易出现非物理负值解，从而导致计算不稳定，甚至失败。即便是激波捕捉性能优异的WENO类方法，也无法完全避免强激波模拟中负值解的出现。已有方法更多通过局部降阶、后处理修正等方式尝试消除负值解对计算稳定性的影响，人为干预性强且效果不一而足。因此，可以考虑发展严格保正计算方法，从格式构造的数学原理解决负值解问题，从而有效加强激波捕捉方法的稳定性。本文构建了后验保持解特征[5-6]激波捕捉格式，通过结合高精度低耗散格式和高分辨率激波捕捉格式各自优势，实现光滑区域高精度求解和间断区域抑制数值振荡的目的。因此，该方法本质上属于混合格式（hybrid scheme）[7-9]。本文在有限体积方法和间断Galerkin方法框架下分别发展了基于该格式的高分辨率激波捕捉方法，并且验证了发展的方法在可压缩湍流问题、低压低密度问题中的初步应用。1　控制方程和离散方法1.1　Euler方程和间断Galerkin（DG）离散方法Euler方程可写成如下形式∂ρ∂t+∇⋅ρu=0∂ρu∂t+∇⋅ρuu+pδ_=0∂ρe∂t+∇⋅(ρe+p)u=0 (1)式中, ρ为密度；u为速度矢量；p为静压；δ_为Kronecker张量； e 为单位质量能量。将上述Euler方程写成如下形式∂Q∂t+∇⋅fcQ=0 (2)式中，Q为守恒型变量；fc为无黏通量项。引入有限元解Qh和多项式测试函数 vh ，通过分部积分的方式可以获得DG离散的控制方程∫Kνh∂Qh∂tdΩ+∫∂Kνhfc⋅n⃗dσ-∫K∇νh⋅fcQhdΩ=0 (3)式中，Qh=∑l=1N(p)νh,lQl(t)；Ql为DG解自由度；N(p)为自由度数目。1.2　Navier-Stokes（N-S）方程和有限体积（FV）离散方法不考虑源项的三维可压缩N-S方程可写成如下形式∂Q∂t+∇⋅fcQ=∇⋅fv(Q,∇Q) (4)式中，fv为黏性通量项；∇Q代表变量梯度。基于结构网格FV方法的离散控制方程形式可写成dQ¯i,j,k(t)dt+f^c,ijk+1/2- f^c,ijk-1/2=f^v,ijk+1/2-f^v,ijk-1/2 (5)式中，Q¯i,j,k为单元平均值；f^c和f^v分别为单元界面数值无黏通量和数值黏性通量；ijk±1/2代表了单元三个方向的界面。2　后验保持解特征重构方法2.1　基于后验测试的混合激波捕捉格式参考文献[10]提出了MOOD（multi-dimensional optimal order detection）方法，通过多步的后验测试和降阶重算，MOOD方法能够实现很多优良的数学特性，如DMP（discrete maximum principle）、保正性等。基于MOOD后验测试理念，研究[11-13]发展了多种类型混合格式，包括迎风重构类型[12]、中心重构类型[11]等。这里，首先给出基于MOOD后验测试构建混合格式的思路。首先，采用某种低耗散格式求解方程（5）可以得到ΔQ*Δt+Rhscentral(Qn)=0 (6)式中，Q*称为候选解。这里采用添加人工黏性的改进型中心格式（central scheme）[11]获得候选解。然后判断候选解是否满足某种测试准则：满足测试准则的单元，候选解就作为下一时刻解；不满足测试准则的单元，采用更加稳定的激波捕捉格式进行重算（recompute），同时重算单元相邻单元进行残差重组（reassemble）。对于测试准则,本文选取包含小量的DMP准则[14]，如式（7）、式（8）所示min(Q¯mn(tn))- δQ¯ijk*(tn+1)max(Q¯mn(tn))+δ (7)δ=max10- 5,10- 4·max(Q¯mn)-min(Q¯mn) (8)除了DMP准则，测试准则中也包括密度和压力为正值的物理准则。可以看到，上述方法通过后验测试准则确定单元界面的数值通量最终将由计算候选解的中心格式还是重算过程的激波捕捉格式来计算，因此，方法属于混合格式。2.2　保持解特征重构格式BVD原理[15]要求尽量减小通过重构得到的网格边界两侧的物理量之前的差，从而能够有效地控制黎曼求解器中的数值黏性。通过选用适当候补函数，建立多步BVD选择过程，可以构建一类高分辨率激波捕捉格式[16-19]，光滑区域利用高阶重构降低耗散，间断区域有效抑制数值振荡并提高间断分辨率。参考文献[5]针对双曲守恒律方程系统实现了该类激波捕捉格式的解保正性，他们将方法总结为基于BVD的解特征保持方法。参照上述工作，这里简要给出本文格式构造中涉及的候补函数。候补重构函数（1）：五阶线性迎风重构qi+1/2L=130q¯i-2- 1360q¯i-1+4760q¯i+920q¯i+1- 120q¯i+2 (9)候补重构函数（2）—二阶THINC重构qi(x)=q¯min+q¯m21+θtanhβx-xi-1/2xi+1/2-xi-1/2-x˜i (10)式中，多步BVD选择算法以及参数选择等均与参考文献[16]保持一致，这里就不再赘述。2.3　后验保持解特征重构方法利用2.1节后验测试的概念，结合2.2节保持解特征重构格式，可以构建后验保持解特征重构方法，具体实现过程如下。（1）采用中心型格式求解式（6），获得候选解Q*。根据后验测试准则，标识重算单元和重组单元。（2）对于重算单元，采用保持解特征重构方法进行重算。这里以一维重构为例，介绍基于BVD原理对重算单元实现保持解特征的重构过程。假设目标单元i为重算单元，那么首先获得单元j（j =i-4, i-3, …, i+4）上采用五阶线性迎风重构（见式（9））和二阶THINC重构（见式（10））得到的单元界面变量，记作qj±1/2L,R。按照式（11）计算单元总变差（TBV）TBVj=qj-1/2L-qj-1/2R+qj+1/2L-qj+1/2R (11)式中：j =i-3, i-2,…, i+3。通过多步BVD选择算法（详细过程可参见文献[16]），获得单元j（j =i-1,i和i+1）的重构函数，进而给出单元界面i±1/2两侧重构变量，应用近似黎曼求解器，即可获得单元i的通量变化。对于重组单元，针对和重算单元相接的界面，用上述保持解特征重构计算的界面通量替换中心型格式计算的界面通量，然后和其他界面通量重新组合获得单元新的通量变化。最后，利用新的通量可以得到重算和重组单元上的更新解。为进一步增强数值方法稳定性，可以对重算和重组的单元进行二次后验测试。测试准则仍然采用修改的DMP准则和密度压力为正的物理准则，但DMP准则中的小量见式（8），调整为10-2和10-1，并且不满足二次后验测试的单元将直接使用一阶格式计算，从而确保解的保正性。数值经验表明，多步BVD选择算法能够有效抑制振荡，因而可以有效避免非物理负值解出现。对于本文考虑的算例，几乎所有单元在大部分计算时间都不需要启用一阶格式。依据作者经验，如果能有效控制计算过程启用一阶计算的单元数量，能够在实现保正计算的同时避免数值方法精度受到过多影响。在下文的3.1~3.4节，算例验证部分采用了多种激波捕捉格式。为便于区别，将后验保持解特征方法用MOOD+BVD或MOOD+BVCD代表，原始方法用BVD或BVCD代表，其中BVCD代表的是BVD算法引入低耗散中心型（Central）重构方法，BVCD代表的方法可以参见文献[20]，五阶Mapped WENO用WM5代表。此外，采用上述二次后验测试方法的算例会特别说明。2.4　保持解特征DG/FV方法混合DG/FV方法[14,21-22]的基本思想是，依据激波指示器（shock indicator），判断流场单元处于光滑区域还是间断区域，将光滑区域的单元标记为DG单元，将间断区域的单元标记为FV单元。设计激波指示器可以根据当前时刻DG解的信息，也可以利用后验计算得到的下一时刻解的信息[14]。本文采用的激波指示方法[23]利用了BVD选择重构函数的信息。具体来说是根据BVD选择线性多项式还是THINC重构，选择前者的单元倾向处于光滑区，选择后者的单元倾向处于间断区，具体方法请参见文献[23]，这里就不再进一步展开。为充分利用DG解多项式多自由度信息，FV单元通常划分为多个亚格子单元，而在亚格子单元上，可以采用成熟的激波捕捉格式来求解。本节重点介绍针对混合DG/FV实现保持解特征重构以及密度压力保正性计算的方法。首先，本文采用参考文献[23]的方法构建亚格子单元，该方法中亚格子的数目和DG解自由度的数目相同，并且亚格子解定义为DG多项式解在亚格子上的单元平均值，这不仅符合守恒律的要求，对于准确求解间断也是非常重要的。然后，在构建的亚格子单元上，应用保持解特征重构方法。需要注意的是，此处并不采用中心型格式去计算候选解，而是利用修改的BVD重构过程[23]，获得变量重构信息，直接标识DG单元和FV单元。在FV的亚格子单元上，采用保持解特征重构格式，计算亚格子解。然后，研究给出DG/FV方法实现激波亚格子分辨的原理。如图1所示，在FV单元上划分了4个亚格子单元。假定亚格子单元i包含间断，根据保持解特征重构原理，在单元i上应选择间断分辨高，即压缩性强的非线性激波捕捉方法。以THINC格式为例，如图1所示，在单元i上BVD方法会选择更大β值（βl）的THINC重构，而在非间断单元上，BVD方法会选择更小β值（βs）的THINC重构。这样就保证了数值方法在间断附近具有足够低的耗散（即激波捕捉更加锐利）。类似地，在BVD选择过程中还可以加入2.2节中的五阶线性重构，进一步降低数值耗散。因此，通过在亚格子单元上实现间断高分辨率捕捉，DG/FV方法能够实现激波亚格子分辨。10.12050/are20220401.F001图1保持解特征DG/FV方法实现激波亚格子分辨示意图[24]Fig.1Illustration of subcell resolution by DG/FV method with solution property preserving property[24]最后，为了获得严格的保正数学特性，在标识DG单元和FV单元过程中，加入了对DG单元积分点位置的密度和压力检测。如果发现任一积分点出现负值，该单元标记为FV单元。同时，对FV单元的亚格子界面处的密度和压力进行检测，如果重构得到的密度或压力为负，那么该亚格子单元直接采用THINC重构（β=1.1）计算。最后，为了保证下一时刻解严格为正，采用后验测试方法，测试准则只应用密度压力为正的物理准则。如果亚格子单元候选解出现负值，则重算过程中直接采用一阶格式重算，从而确保解的保正性。3　数值算例数值算例分为两大部分：3.1~3.4节主要考察MOOD+BVD方法在FV离散框架下求解可压缩湍流问题的性能；3.5~3.6节主要考察保持解特征DG/FV方法求解低压低密度问题的性能。3.1　可压缩湍流算例：Taylor-Green涡该算例本质上是不可压流动，实际模拟时给定马赫数Ma=0.08，雷诺数Re=1600。图2和图3给出了不同网格上使用不同格式得到的体平均动能和动能耗散率随时间演化的结果，参考解为文献[25]中密网格上得到的DNS解。可以看到，在同等网格下，本文方法要好于原始格式（BVD）的结果，在963的网格上，本文方法已经能获得接近DNS的结果，并且要好于WM5格式在1283网格上的结果。10.12050/are20220401.F002图2Taylor-Green涡，动能随时间变化曲线Fig.2Evolutions of kinetic energy for the Taylor-Green vortex10.12050/are20220401.F003图3Taylor-Green涡，动能耗散率随时间变化曲线Fig.3Evolutions of dissipation rate for the Taylor-Green vortex3.2　可压缩湍流算例：各向同性衰减湍流本小节选取可压缩各向同性衰减湍流算例，考察方法捕捉激波和分辨小尺度结构的能力，算例的详细设置可参考文献[4]和[26]。首先选取湍流马赫数Mat=0.6，泰勒微尺度雷诺数Reλ=100的算例[4]，此时流场中已出现局部超声速区域，具有较强的压缩性。图4、图5给出了不同方法计算的动能和温度随时间的演化结果，后验保持解特征方法仍然要优于原始格式在同等网格以及WM5格式在密网格上的结果，图6则给出了目前方法在不同网格上得到的动能能谱。随着网格加密，能谱能很好地逼近DNS解，并且963网格的结果已经基本和1283网格的结果重合。然后选取了压缩性更强的Mat=1.2，Reλ=72的算例[26]，该工况下需要采用第2.3节二次后验方法避免压力或密度出现负值。图7~图10给出方法在1923的网格得到的体平均量随时间演化的结果，参考解为文献[26]在3843网格得到的DNS解。图11、图12则给出2563网格上t=τ （τ代表大涡翻转时间尺度）时刻流场Q图和速度散度（dilation）的结果，从Q图可以看到，局部流场出现较高的当地湍流Ma数，速度散度的结果也显示了狭长的激波列（shocklet）结构和强膨胀区域无序地分布在流场中，类似结果可以参见文献[26]中的图9。这些强压缩和膨胀结构对激波捕捉格式的鲁棒性提出了更高的要求，目前的方法能够模拟强可压缩湍流，并且在相对较粗的网格获得了和DNS基本吻合的结果。10.12050/are20220401.F004图4Mat=0.6各向同性衰减湍流，动能随时间变化曲线Fig.4Evolutions of kinetic energy for isotropic turbulence with Mat=0.610.12050/are20220401.F005图5Mat=0.6各向同性衰减湍流，温度随时间变化曲线Fig.5Evolutions of temperature for isotropic turbulence with Mat=0.610.12050/are20220401.F006图6Mat=0.6各向同性衰减湍流，t=4τ时刻动能能谱Fig.6Velocity spectra at t=4τ for isotropic turbulence withMat=0.610.12050/are20220401.F007图7Mat=1.2各向同性衰减湍流，密度脉动随时间变化曲线Fig.7Evolutions of density fluctuations for isotropicturbulence with Mat=1.210.12050/are20220401.F008图8Mat=1.2各向同性衰减湍流，solenoidal耗散率随时间变化曲线Fig.8Evolutions of solenoidal dissipation rate for isotropicturbulence with Mat=1.210.12050/are20220401.F009图9Mat=1.2各向同性衰减湍流，dilatational耗散率随时间变化曲线Fig.9Evolutions of dilatational dissipation rate for isotropicturbulence with Mat=1.210.12050/are20220401.F010图10Mat=1.2各向同性衰减湍流，pressure-dilatation随时间变化曲线Fig.10Evolutions of pressure-dilatation for isotropicturbulence with Mat=1.210.12050/are20220401.F011图11Mat=1.2各向同性衰减湍流，t=τ时刻Q图（Mat数染色）Fig.11Iso-surface of Q results colored with Mat at t=τ for isotropic turbulence with Mat=1.210.12050/are20220401.F012图12Mat=1.2各向同性衰减湍流，t=τ时刻无量纲化的dilation结果Fig.12Counters of non-dimensional dilatation at t=τ forisotropic turbulence with Mat=1.23.3　可压缩湍流算例：超声速槽道湍流选取超声速槽道湍流考察方法模拟壁湍流的能力。主要计算参数设置为Maτ=0.082（Ma=1.5）, Reτ=222，下标τ代表摩擦速度。算例的详细设置可参考文献[27]。本算例采用网格数目为130×100×80，图13、图14给出了统计平均量的计算结果，方法能够给出较为准确的预测。图13中平均温度的结果和谱方法的结果[28]存在偏差，根据可压缩规范壁湍流温度和速度对应关系[29]，速度场的结果也会存在偏离。这里指出，参考文献[27]中隐式大涡模拟和直接数值模拟得到的温度计算结果也和谱方法的结果有类似偏差，参见文献[27]的图3。同样，高阶雷诺应力在具体量值上仍然和参考文献[28]谱方法的结果存在一定偏差，这一点也和参考文献[27]中部分方法的计算结果类似。表1给出了图14中的雷诺应力极值，将不同方法和参考解的偏差作对比可以发现，目前方法虽然采用较少的网格，但仍然可以获得和参考文献[27]采用高阶有限差分（FD）相近的偏差结果。最后，图15瞬时流场Q图也表明了方法能够给出较为丰富的流场结构，这些对比验证了后验保持解特征方法能够适用超声速壁湍流高精度模拟。10.12050/are20220401.F013图13超声速槽道流，平均流向速度、密度和温度分布Fig.13Mean streamwise velocity, density and temperature forsupersonic channel flow10.12050/are20220401.F014图14超声速槽道流，平均雷诺应力分布Fig.14Mean Reynolds stresses for supersonic channel flow10.12050/are20220401.T001表1超声速槽道流，平均雷诺应力极值对比Table 1Comparisons on mean Reynolds stresses forsupersonic channel flowu'u'v'v'w'w'谱方法DNS[28]，参考解7.4380.6451.102六阶FD TENO[27]，1293网格7.8950.5870.963FV MOOD+BVD7.0780.5900.95910.12050/are20220401.F015图15超声速槽道流，MOOD+BVD计算瞬时流场Q图（Ma数染色）Fig.15Instantaneous Q results colored with Ma numbeobtained by MOOD+BVD3.4　可压缩湍流算例：高亚声速空腔流动本算例涉及旋涡、强剪切、激波膨胀波干扰等复杂流动现象，主要参数设置包括来流Ma∞=0.8,基于空腔长度ReL=8.6×105，空腔长深比L/D≈0.42。算例的详细设置可参考文献[30]。采用隐式大涡模拟，计算网格约1.2百万，属于相对较粗的网格。计算过程中应用了第2.3节的二次后验测试。图16给出流场瞬时Q准则等值图和z=0平面密度梯度云图。可以看到流场中存在压力波的来回运动和多次反射，Q图显示了混合层中出现三个组织较好的大尺度旋涡结构，这三个旋涡结构在展向具有一定的二维特性，旋涡S1的尺寸最大，S2尺寸稍小并与S1的结构相似，旋涡S3的尺寸很小，几乎由单根涡管构成。图17给出了唇口下游x=0.4站位处的流向速度、横向速度及雷诺应力分布，参考结果为Larchevêque等[31]以及Forestier等[32]试验结果。目前方法在相对较粗的网格上仍然能合理地预测混合层内的速度和雷诺应力分布。10.12050/are20220401.F016图16高亚声速空腔流动，瞬时Q准则等值图和z=0平面密度梯度云图Fig.16Instantaneous Q results and density gradient at z=0 for cavity flow10.12050/are20220401.F021图17高亚声速空腔流动，x=0.4站位速度和雷诺应力分布Fig.17Mean velocity and Reynolds stresses at x=0.4 for cavity flow3.5　低压低密度算例：激波绕后向台阶流动激波绕后向台阶是测试方法保正性的典型算例。主要计算参数设置为Ma=5.09，后向台阶为90°，算例的详细设置可参考文献[24]~[33]。本算例采用边长为1/32的四边形网格，此处说明该算例计算过程中并没有启动2.4节提到的后验一阶计算过程。图18给出了t=2.3时刻密度解，采用亚格子解显示，图19给出解的局部放大图并显示了主网格。计算结果表明，方法能够在亚格子尺度上获得基本无震荡的解，并且能够在一个单元内部实现激波捕捉。10.12050/are20220401.F017图18激波绕后向台阶，t=2.3时刻密度解（以亚格子解显示）Fig.18Density solutions at the subgrid level at t=2.3 for shock diffraction over a step10.12050/are20220401.F018图19激波绕后向台阶，图18的局部放大图（网格线显示主网格）Fig.19Enlarged view of Fig. 18 with gridlines indicating main element3.6　低压低密度算例：Mach 800 喷流该算例对算法的保正性提出了更高的要求，有研究指出[34]，在马赫数非常高的情况下，不做任何修改的二阶TVD格式也无法保证获得正值的密度和压力解。主要计算参数设置为喷流速度为800，其他详细设置可参考文献[24]和[34]。本算例采用网格数目为400×600，图20、图21给出了t=0.002时刻密度解和识别的FV单元，密度解和文献[34]结果保持一致，FV单元基本围绕间断结构并有较好的对称性，证明了方法具有良好的数值特性。10.12050/are20220401.F019图20Mach 800 喷流，密度解（区间[0.0159, 10.8]，对数率显示）Fig.20Density solution between [0.0159,10.8] with color bars showing logarithmic values for Mach 800 jet flow10.12050/are20220401.F020图21Mach 800 喷流（红色单元为FV单元）Fig.21Mach 800 jet flow (FV cells are showed with red color)4　结论本文发展了后验保持解特征重构方法。该方法本质上属于混合型计算方法，结合了高精度低耗散格式和高分辨率激波捕捉格式各自优势，实现光滑区域高精度求解和间断区域抑制数值振荡的目的。方法在FV和DG框架通过不同的途径进行了实现：在FV框架下，通过后验测试和重算，耦合中心型格式和保持解特征重构实现高精度计算；在DG框架下，通过将流场中的“光滑”和“间断”区域标记为DG和FV单元，在FV单元实现亚格子保持解特征重构，并利用后验计算实现密度压力保正性，从而获得保持解特征DG/FV方法。本文采用可压缩湍流问题、低压低密度问题分别考察方法在复杂流动模拟、保正性计算方面的性能，得到主要结论如下：（1）后验保持解特征方法（MOOD+BVD）相比原始保持解特征方法（BVD）具有更低的数值耗散，在粗网格上求解可压缩湍流问题性能更加优异，部分算例验证了方法要优于其他高阶格式在更密网格上的计算结果，并且方法兼具低耗散和强鲁棒性的特点，能够适用于强可压缩湍流、壁湍流问题的隐式大涡模拟。（2）针对低压低密度问题，保持解特征DG/FV方法具备压力密度保正性，能够在亚格子尺度上获得基本无震荡的解，并且在单元内部实现激波捕捉。