0 引 言
激波是可压缩流体动力学中最重要的物理现象，模拟激波自然也构成计算流体力学(CFD)主要的研究内容。激波的模拟方法分为捕捉和装配两种。1953年，Lax提出自动捕捉双曲型守恒率方程间断解的一阶格式，在应用中发现数值解存在非物理波动；1959年，Godunov基于弱解理论，利用间断满足的Riemann分解的特性，构造出物理背景清晰计算稳定的一阶格式[1]。由于这两种格式都只有一阶精度，计算得到的激波需要若干网格点过渡(所谓激波较宽)。直到1983年Harten提出TVD格式才实现二阶精度，使得激波分辨率得以提高。其后发展了许多高精度格式，如NND、ENO、WENO等[2]，近几年三阶乃至更高阶的WCNS也得到实际应用[3, 4, 5]。1966年，Moretti根据激波上下游流动参数满足R-H关系式的特点提出激波装配法[1]，高精度高效率地模拟出钝体头部脱体激波，二维情况下装配法使用10×16个网格得到的压力分布的精度不逊于当代捕捉法使用385×513个网格的精度[6]，极大地推动了当时超声速流动的应用研究，直到上世纪80年代后期依然是主流算法[7]。然而，这种算法需要事先确定激波位置，应用受到极大限制，随着高精度捕捉法的发展，装配法的研究热情逐渐消退。
激波空间尺度只有几个分子的平均自由程[8]。在海平面标准状况下，马赫数M1=2.95的正激波厚度大约为6.6×10-8m，而常见的飞行器在10m量级，目前航天航空领域CFD应用中10-3m量级网格很少见，因此采用捕捉方法模拟出来的激波宽度要比实际的激波厚度大得多。激波内部处于热力学非平衡状态，而建立Euler或NS方程时采取了热力学局部平衡假设、密度不存在剧烈变化的Stokes假设以及采用了忽略松弛时间的状态方程等，从理论角度看，从Euler或NS方程出发捕捉激波，即使未来几十年后计算机发展使得激波厚度内布置若干网格成为常态，只要过度区有1个网格点，也会产生非物理波动。描述激波最完善的数学理论是双曲型守恒率方程的弱解理论。它把激波看作无厚度的间断，流动参数在此发生突跃。
由于捕捉法计算激波存在诸多理论缺陷，文献[6]中认为捕捉法经过36年(1966-2001)的发展依然无法超越装配法。近几年国外一些学者重新审视装配法，解决了高精度格式应用于高马赫数流动时在激波附近出现波动的问题[9, 10]。Paciorri等人采用在捕捉法得到的流场基础上装配浮动激波的新思路突破传统装配法的限制，极大推动了应用领域拓展[11, 12, 13]。浮动激波装配法通过局部网格重构来形成内部激波边界，由于激波运动以后需要在背景网格上重新定位和插值，流场信息传递使其难以用于非定常流动[11, 12, 13]。
本文利用文献[14]建立在ALE(Arbitrary Lagrangian-Eulerian)方程基础上的非结构动网格技术发展能够用于非定常流动的装配法。 1 ALE形式控制方程及其计算方法
采用ALE描述的三维可压缩非定常N-S方程的积分形式：
其中，Ω为控制体，∂Ω为控制边界，n为控制体边界外法向向量，dV为体积微元，dS为面积微元，xc为网格运动速度。
计算方法采用文献[14]中介绍的基于混合网格的有限体积方法，对于守恒变量存贮在单元中心的格心格式，空间离散第i个控制体：
其中，Fk、nk、Sk分别为第k个面元的通量、外法向单位矢量和面积，Nf为控制体表面个数。假设物理量在单元内线性分布得到空间二阶精度的有限体积法格式，第i个控制体的第k个表面中心点的变量值： 其中，∇Qi是单元内物理量梯度，采用Green公式计算梯度，rik是单元中心到表面中心的矢径，φik是为了抑制激波附近可能出现的非物理振荡而引入的限制器，本文采用Venkatakrishman限制器。两侧单元在表面k处的重建值形成物理量的间断，将其视为近似的一维Riemann问题，可以利用Riemann问题的各种求解方法计算单元k处的无粘通量。在本文软件中包含有Van Leer、Roe、AUSM、HLLC等多种分裂格式。时间离散采用二阶时间精度的四步Runge-Kutta方法。计算中稳定性要求时间步长满足CFL条件： 式中，U为流体速度，ai是当地声速，Vi是单元i的体积。对四步Runge-Kutta格式，CFL≤2
。全场时间步长是所有当地时间步长的最小值，即Δt=min(Δti)。
计算时还需要考虑离散几何守恒律(Discrete Geometric Conservation Law—DGCL)；很多时候网格变形算法需要结合网格重构实现大变形或大位移，这又涉及到新旧网格之间的信息传递。 1.1 离散几何守恒律
对于均匀流场，NS方程的质量方程：
在离散空间，从n推进到n+1时刻，微元体积从Vn变化为Vn+1，对于显式算法，第i个面Sin以速度xci运动形成的体积增量记为ΔVin=xci·niSin，对于全隐算法，第i个以速度xci运动形成的体积增量记为ΔVin+1=xci·niSin+1，显然有： 其中，Ns表示单元所含面元个数。这种几何计算误差必然影响流场计算，因此需要构建DGCL来消除。
根据所修正的参数把国内外文献构造的DGCL分为三类算法：第一类修正网格面积，采用平均面积：
代替Sni或Sn+1i，这种算法最早和几何守恒律概念[15]一起提出，难以应用于涉及2个以上时间层的高阶隐格式；第二类修正网格运动速度，如果时间二阶格式离散左端项：
在文献[19]中通过数学模型证明离散几何守恒律是流场参数假设为常数条件下ALE坐标形式流体力学方程的退化形式，不是独立的新约束条件，不满足几何守恒律引起非物理解的机理是在网格变形过程中体积增量与面积运动形成的体积不等而产生误差。在非定常流场的计算过程中，随着时间推进步数增加，误差不断累积，因此必须构建几何守恒律算法来消除这种误差。基于这种认识分析以上3类DGCL算法，单一修正离散方程右端项，必然导致与求解流场的空间离散格式相关。本文作者提出修正离散方程左端项中体积的新思路[20]，采用体积V″代替Vn进行计算：
网格重构以后，采用插值传递新旧网格之间流场信息，线形插值导致激波分辨率降低，ENO插值引起非物理波动，很难构造适合激波的3阶插值。为了避免新旧网格之间流场信息传递过程产生新误差，文献[22]提出一种基于动网格的信息传递方法，下面以二维问题为例对其基本原理进行介绍。
问题如图 1(a)所示，已知n时刻旧网格中O点的流场参数，需要得到n+1时刻新网格的P点和Q点流场参数。如图 1(b)所示，在旧网格的基础上时间方向推进Δt，计算得到n+1时刻流场参数的同时，采用动网格技术保持O点及周围所涉及的网格不变形但是同一速度运动，使得O点平移到P点；同样，在旧网格基础上时间方向推进Δt，将O点平移到Q点。由于P和Q所在的空间位置正好是新网格的点，如图 1(c)所示，因此通过动网格技术得到n+1时刻的流场参数，在此基础上放弃旧网格，采用新网格继续计算。可以证明，对于以上时间空间二阶精度的有限体积法格式，这种方法不引入误差，因此称为高精度信息传递算法[14]。时间显式格式的稳定性限制网格运动速度，在多维条件下少量网格需要多步时间推进才能实现统一时间Δt流场，为此文献[21]提出了一些改进措施。
 |
| 图 1 信息传递原理图Fig. 1 Schematic of information transfer |
1990年Batina首次将弹簧比拟法用于三角形网格的变形[23]。它的原理是将连接网格边视为弹簧，整个计算区域网格构成弹簧系统，在弹簧系统内部节点(网格点)列出力的平衡方程组，边界发生运动或变形后使弹簧网络受到拉伸和挤压，整个弹簧系统受力发生变化，通过求解弹簧网络的平衡方程来更新网格节点的位置。这种方法只考虑了弹簧拉压，扭转变形易出现网格折穿(snap-through)。1998年Farhat等人在二维网格顶点处施加扭转弹簧增加抗扭转能力取得实效，但在推广到三维网格遇到困难[24, 25]。2000年Blom等人根据三角形单元边对应的顶角来修正线弹簧刚度系数，提出半扭转弹簧方法[26]，2003年郭正和刘君针对四面体，采用不包含该边的二面角作为顶角，把这种方法推广到了三维情形，同时通过求解固体导热方程得到的内部节点温度作为参数来增加动边界附近网格层的弹簧刚度系数，提高了网格变形能力[27]。本文采用[27]发展的弹簧比拟法，所求解的力平衡方程组对角占优，采用Jacobia迭代能较快收敛。 2 基于非结构动网格的激波装配法
半圆柱前形成脱体激波结构的无粘流场，计算网格如图 2所示。为了加快收敛速度和提高非定常流场计算的初值精度，本文借鉴文献[11, 12, 13] 的新思路，在捕捉法流场基础上按照如下算法确定装配法的初始激波位置：(1)寻找全场密度梯度最大单元，记为E0，然后找出该单元顶点周围密度梯度最大的单元，记为E1；(2)增加排除标记过的单元、到Ei-1格心大于Ei与Ei-1格心距离约束，找出与Ei单元共点的密度梯度最大的单元，记为Ei+1，直到边界；(3)采用Bézier方法进行拟合单元格心坐标，得到的初始激波位置如图 2中曲线所示，用作装配法生成网格的边界。
 |
| 图 2 计算网格Fig. 2 Computational grids |
对于装配法重新生成的新网格，除了激波边界网格外，按照上面介绍的信息传递方法给出流动参数。对于激波边界处新网格的流场参数，先在旧网格中定位，然后取旧网格中对应单元的相邻单元中最大压力作为激波后压力p2，结合激波前压力p1可以计算激波在静止环境中的运动速度：
在均匀来流速度为v1的风轴系下，激波相对运动速度：
其中，vn1=μ·v1·n，μ=signv1·n。
本文采用基于非结构动网格的激波装配法，不需要像文献[11, 12, 13]那样处理浮动激波，网格按照以上激波速度进行运动，即xc=vs。为使激波边界新网格的流场参数相容，引入相对激波马赫数M1=vn1+xc·n，根据R-H关系式确定波后密度和速度：
其中vτ1为均匀来流在界面上切向分量。
按照以上捕捉法方法编制程序模拟马赫数大于6的钝体无粘绕流场。如果初场取自由来流参数，先用一阶迎风格式得到收敛流场，然后用二阶精度格式进一步提高分辨率。对于如图 2所示较稀疏的计算网格，随着马赫数增大激波增强，采用二阶格式后需要探索调整限制器参数、计算时间步长等才能避免发散或振荡，到了马赫数20时，采用Roe、AUSM、HLLC等空间格式很难收敛到稳定解，基于非结构动网格的激波装配法没有出现不收敛现象。图 3为马赫数等于20时采用二阶Vanleer分裂格式捕捉和装配得到的马赫数等值线图的比较。物面压力和Lyubimov的经典文献[28]、文献[7]进行比较，符合很好，由于篇幅有限，这里不再赘述。
 |
| 图 3 马赫数等值线图Fig. 3 Mach iso-contours |
本文采用基于非结构动网格的激波装配法成功地计算出超高马赫数下的定常流场。下面对如图 4所示的整流罩分离动态过程进行模拟。整流罩初始开启角度为5°，在气动力作用下整流罩绕法兰旋转打开，达到30°时整流罩分离进行六自由度自由运动。该流场存在激波、膨胀波等诸多复杂的非定常流动现象。在计算过程两种方法均需要多次网格重构，重构生成新网格在激波后区域尺度接近和数目相当，采用前面介绍的新旧网格之间信息传递方法给出新网格上参数。文献[14]提供了采用捕捉法模拟超声速多体分离干扰流场的大量验证算例，本文以捕捉法的结果作为考核标准，模拟飞行马赫数6条件下整流罩动态分离过程。对比整流罩分离过程中捕捉法和装配法得到的不同时刻的马赫云图，如图 5所示，以及质心和姿态角随时间变化曲线，如图 6、7所示。可以看出，采用装配法得到的流场结构、气动力驱动运动特性与捕捉法符合较好。对比表明本文采用的基于非结构动网格的激波装配法可以很好地模拟非定常激波。
 |
| 图 4 激波装配计算初始网格Fig. 4 Initial computational grids for shock-fitting |
 |
| 图 5 Ma=6整流罩分离过程中不同时刻马赫云图Fig. 5 Mach contours change versus time in the process of fairing separation at Ma= 6 |
 |
| 图 6 整流罩质心运动( S-C:捕捉法; S-F:装配法)Fig. 6 Fairing′s barycenter movement versus time (S-C:shock-capturing; S-F:shock-fitting) |
 |
| 图 7 整流罩姿态角变化(A1:上片角度;A2:下片角度)Fig. 7 Fairing′ angle of pitch change versus time (A1:angular of upper panel A2:angular of lower panel) |
比较了捕捉法一阶迎风和二阶格式的收敛解得到的初始激波位置，几乎没有区别，这一现象定性合理，因为不同精度格式对于已经分辨出来的激波主要差异在于宽度(网格点数)，最大梯度的位置差异较小。如果格式在激波区光滑过渡，尽管宽度不同，但是若干网点后必然恢复到激波后参数，基于这样的认识，提出采用旧网格的相邻网格单元中最大压力作为激波边界网格压力的计算方法来消除不同精度格式对初场的影响。通过图 6和图 7的对比发现，采用两种方法得到的计算结果基本一致。由于激波装配法无法消除捕捉法提供初场带来的误差，可以得出结论：采用激波装配方法的计算过程没有引入新的误差。
随着马赫数增加，二阶格式的捕捉法难以正常运行，需要不断在计算发散以后进行人工干预，基于一阶迎风格式捕捉到激波位置作为初始边界的装配法计算没有任何问题，图 8给出马赫数20时分离过程中不同时刻的压力云图。
 |
| 图 8 Ma=20整流罩分离过程中不同时刻压力云图Fig. 8 Pressure contours change versus time in the process of fairing separation at Ma=20 |
本文提出了以一种新的基于非结构动网格技术的激波装配方法，采用一阶格式捕捉法得到的流场确定装配法所需的初始激波位置，将激波装配法成功地应用到非定常流动产生的运动激波的模拟，复杂外形产生的非规则形状激波的模拟，解决了高精度格式捕捉强激波遇到的稳定性问题。
| [1] | Li S B. Theory of dissipative conservation scheme[M]. Beijing: Higher Education Press, 1997. (in Chinese) 李松波. 耗散守恒格式理论[M]. 北京: 高等教育出版社, 1997. |
| [2] | Zhang H X,Shen M Y. Computational fluid dynamics-theory and application of difference scheme[M]. Beijing: National Defense Industry Press, 2003. (in Chinese) 张涵信, 沈孟育. 计算流体力学—差分格式原理和应用[M]. 北京: 国防工业出版社, 2003. |
| [3] | DENG X G, ZHANG H X. Developing high-order weighted compact nonlinear schemes[J]. Journal of Computational Physics, 2000, 165(1):22-24. |
| [4] | Li Q, Sun D,Zheng Y K, Zhang H X. On a class of center-typed third order difference scheme orienting to engineering utilizations[J]. Acta Aerodynamics Sinica, 2013, 31(4): 466-493. (in Chinese) 李沁, 孙东, 郑永康, 张涵信. 一类中心型三阶格式及其应用[J]. 空气动力学学报, 2013, 31(4): 466-493. |
| [5] | Zhang S H, Li H, Liu X L. A direct numerical simulation of the complex multi-scale flow with shock, vortex and sound wave[J]. Computer Engineering & Science, 2012, 34(8): 99-107. (in Chinese) 张树海, 李虎, 刘旭亮. 含激波、旋涡和声波的复杂多尺度流动的直接数值模拟[J]. 计算机工程与科学, 2012, 34(8): 99-107. |
| [6] | Moretti G. Thirty-six years of shock fitting[J]. Computers & Fluids, 2002, 31:719-723. |
| [7] | Ye Y D. Numerical simulation of the modified space shuttle orbiterinviscid flowfield[D]. Changsha: National University of Defence Technology, 1991. (in Chinese) 叶友达. 航天飞机简化外形无粘流场的数值模拟[D]. 长沙: 国防科技大学, 1991. |
| [8] | Tong B G, Kong X Y, Deng G H. The gas dynamics[M]. Beijing: Higher Education Press, 2012. (in Chinese) 童秉纲, 孔祥言, 邓国华. 气体动力学[M]. 北京: 高等教育出版社, 2012. |
| [9] | Prakash A, Parsons N, Wang X, Zhong X. High-order shock-fitting methods for direct numerical simulation of hypersonic flow with chemical and thermal nonequilibrium[J]. Journal of Computational Physics, 2011, 230: 8474-8507. |
| [10] | Najafi M, Hejranfar K, Esfahanian V. Application of a shock-fitted spectral collocation method for computing transient high-speed inviscid flows over a blunt nose[J]. Journal of Computational Physics, 2014, 257: 954-980. |
| [11] | Paciorri R, Bonfiglioli A. A shock-fitting technique for 2D unstructured grids[J]. Computers & Fluids, 2009, 38: 715-726. |
| [12] | Paciorri R, Bonfiglioli A. Shock interaction computations on unstructured, two-dimensional grids using a shock-fitting technique[J]. Journal of Computational Physics, 2011, 230: 3155-3177. |
| [13] | Bonfiglioli A, Grottadaurea M, Paciorri R, Sabetta F. An unstructured, three-dimensional, shock-fitting solver for hypersonic flows[J]. Computers & Fluids, 2013, 73: 162-174. |
| [14] | Liu J,Bai X Z, Guo Z. The computational method of unstructured moving grids[M]. Changsha: National University of Defence Technology Press, 2009. (in Chinese) 刘君, 白晓征, 郭正. 非结构动网格计算方法[M]. 长沙:国防科技大学出版社, 2009. |
| [15] | Thomas P D, Lombard C K. The geometric conservation law-a link beween finite-difference and finite-volume methods of flow computation on moving grids[R]. AIAA 1978-1208. |
| [16] | Mavriplis D J, Yang Z. Construction of the discrete geometric conservationa law for high-order time-accurate simulations on dynamic meshes [J]. Journal of Computational Physics, 2006, 213: 557-573. |
| [17] | Farhat C, Geuzaine P. Design and analysis of robust ALE time-integrators for the solution of unsteady flow problems on moving grids[J]. Computer Methods in Applied Mechanics and Engineering, 2004, 193: 4073-4095. |
| [18] | Zhang L P, Wang Z J. A block LU-SGS implicit dual time-stepping algorithm for hybrid dynamic meshes[J]. Computers & Fluids, 2004, 33: 891-916. |
| [19] | Liu J,Bai X Z, Zhang H X, Guo Z. Discussion about GCL for diforming grids[J]. Aeronautical Computing Technique, 2009, 39(4): 1-5. (in Chinese) 刘君, 白晓征, 张涵信, 郭正. 关于变形网格“几何守恒律”概念的讨论[J]. 航空计算技术, 2009, 39(4): 1-5. |
| [20] | Liu J,Xu C G, Liu Y. New scheme for interface boundary of fluid/structure interaction problems[C]//2012 National conference of CFD, 2012: 247-253. (in Chinese ). 刘君, 徐春光, 刘瑜. 数值模拟流固耦合的界面新算法[C]//2012全国计算流体力学会, 2012: 247-253. |
| [21] | Xu C G. Research on numerical simulation and application of complicated jet flow field[D]. Changsha: National University of Defence Technology, 2002. (in Chinese) 徐春光. 复杂喷流流场数值模拟及应用研究[D]. 长沙: 国防科技大学, 2002. |
| [22] | Liu J, Bai X Z, Guo Z, et al. A new method for transferring flow information among meshes[J]. Computational Fluid Dynamics Journal, 2007, 15(4):509-514. |
| [23] | Batina J T. Unsteady Euler airfoil solutions using unstructured dynamic meshes[J]. AIAA Journal, 1990, 28(8): 1381-1388. |
| [24] | Farhat C, Degand C, Koobus B, et al. Torsional springs for two-dimensional dynamic unstructured fluid meshes[J]. Computer Methods in Applied Mechanics and Engineering, 1998, 163: 231-245. |
| [25] | Degand C, Farhat C. A three-dimensional torsional spring analogy method for unstructured dynamic meshes[J]. Computers & Structures, 2002, 80: 305-316. |
| [26] | Blom F J. Considerations on the spring analogy[J]. International Journal for Numerical Methods in Fluids, 2000, 32: 647-668. |
| [27] | Guo Z, Liu J, Qu Z H. Dynamic unstructured grid method with applications to 3D unsteady flows involving moving boundaries[J]. Acta Mechanics Sinica, 2003, 35(2): 140-146. (in Chinese) 郭正, 刘君, 瞿章华. 非结构动网格在三维可动边界问题中的应用[J]. 力学学报, 2003, 35(2): 140-146. |
| [28] | Lyubimov A N, Rusanov V V. Gas flows past blunt bodies, PART II: Tables of the gas dynamic functions[R]. NASA TT-F-715, 197. |