Numerical study of shock wave drag reduction mechanism by nanosecond-pulse laser energy deposition

超声速和高超声速飞行器在飞行时，前方会产生一个弓形激波，波阻大，飞行性能受到极大影响。纳秒脉冲激光具有峰值功率密度高、易于击穿空气形成等离子体的特点，并且具有方向性好、可控性强、可机载和不影响气动外形等特点，在主动流动控制领域具有广泛的应用前景。

激光减阻的基本原理是在飞行器前方沉积激光能量，激光等离子体与飞行器周围流场相互作用改变弓形激波形态，从而减小波阻。目前国内外在纳秒脉冲激光能量沉积减小波阻的机理上存在不同的观点，因此有待进一步详细研究。Myrabo等人^[1]认为激光能量沉积形成的空气锥是波阻减小的原因。Tret'yakov^[2]认为激光能量沉积导致的低密度热尾迹是减阻的原因。Sasoh^[3]和Azarova^[4]等人认为，激光能量沉积诱导的等离子体热核受弓形激波影响后形成的涡结构，是波阻降低的原因。Markhotok^[5]通过建立激波-等离子体团相互作用模型，数值计算了激波结构的变化，认为激波的折射是波阻减小的单一机理，他同时也关注到了激波-等离子体团相互作用时会有涡的产生。文明等人^[6-7]在激波管中利用纹影技术实验研究了激光等离子体与正激波相互作用的特性，给出了流速、激光能量、作用距离等几个关键参数的影响规律。

现有的激光能量沉积减阻数值模拟方法主要存在两方面不足。一方面，采用的能量沉积模型多数为球状形或椭球形的高斯分布模型或均匀分布模型。这些模型带来了适当的简化，使参数化数值研究更容易开展，但这些模型均不能模拟出等离子体热核合理的空间分布。另一方面，大多数的研究均采用量热完全气体模型，且不考虑空气在高温下的离解和电离，与真实情况相差较大，Desai等人^[8]的工作表明考虑真实气体效应对能量沉积减阻的数值模拟有重要意义。

文中针对纳秒脉冲激光能量沉积减小波阻研究中的上述问题，提出了一种采用泪滴形能量沉积模型并考虑空气离解电离的数值模拟方法，分别研究了单脉冲激光能量沉积与正激波相互作用、单脉冲激光能量沉积与钝头体超声速流场相互作用，详细揭示了单纳秒脉冲激光能量沉积减小波阻的机理。

大多数学者都认为在激光能量沉积结束的时刻，等离子体热核已经处于局部热力学平衡状态。相关的实验研究表明，热核在脉冲结束的1 ns内就已经到达局部热力学平衡状态。而激光能量沉积减阻研究均在激光脉冲结束后的微秒量级。在激光能量沉积后，空气增加的热力学内能包括平动能、转动能、振动能和电子势能，需要用平动温度、转动温度、振动温度和电子温度来描述。而采用局部热力学假设之后，可以只用一个温度来描述等离子体热核的状态，从而大大简化了能量方程的数目。

在局部热力学平衡的基础上，对激光能量沉积后热核演化过程进行如下简化：

（1）轴对称、非定常、可压缩的层流流动；

（2）激光脉冲结束时空气各组分的速度为0 m/s；

（3）不考虑空气的热力学非平衡效应、电磁效应和热辐射的影响；

（4）热核中各组分的浓度和物性参数只与温度有关；

（5）满足理想气体状态方程。

Dors在研究激光能量在静止空气中沉积时采用了非对称泪滴形能量沉积模型^[9]。该模型认为热核的初始形状并不是球形或椭球形，而是泪滴形。模型的核心假设是：温度分布在激光入射轴上时呈指数衰减，而垂直于激光入射轴上为高斯分布。

借鉴Dors的能量沉积模型，文中提出的波长1064 nm、脉冲10 ns、单脉冲激光入射能量39.9 mJ（能量沉积约10.6 mJ）的能量沉积模型，可描述如下：

以激光脉冲结束时刻为t = 0时刻。以激光入射方向为x轴方向，垂直于激光入射方向为y轴方向，以热核中点为原点。初始热核在x轴方向上的长度为1 mm。热核初始温度分布为^[10]：

式中：T(x, y)的单位为K；x、y的单位为m。

在一定的环境温度和压力下，激光能量沉积后形成的热核初始温度存在一个饱和值^[11]。当热核达到饱和温度后，继续增大激光能量并不会使热核温度继续升高。文中将饱和温度取为20000 K^[12]。热核初始温度分布如图1所示。

Figure 1.  Initial temperature of hot core

在考虑能量沉积模型和后续数值计算之前，需要对Fluent软件里的空气等离子体参数进行修正。根据前人的研究结果^[13]，300~20000 K范围内空气等离子体的物性参数，包括比热容、粘度、导热系数、相对分子质量和声速等，均可认为是温度的函数。在用Fluent软件进行计算之前，需要把这些数据拟合为多项式函数的形式，然后通过UDF导入Fluent计算过程。图2中给出了空气等离子体的部分物性参数。

Figure 2.  Specific heat capacity and viscosity of air plasma

采用上述能量沉积模型时，热核吸收的激光能量计算公式如下：

式中：x₂，x₁为热核在x轴上的左右边界，m；y₂，y₁为热核在y轴上的上下边界，m；T₀为环境温度，K；h为单位质量气体的静焓，J·kg⁻¹；$ {h}_{\alpha } $为组分α的静焓，J·kg⁻¹；$ {Y}_{\alpha } $为组分α的质量分数。以公式（1）确定的能量沉积模型为例，积分可得沉积能量约为10.1 mJ，与实验测得的10.6 mJ接近。

采用Fluent软件中层流有限速率模型。文中认为常温常压下的空气是由体积分数为79%的${{\rm{N}}_2}$和21%的${{\rm{O}}_{\rm{2}}}$组成，不计其他成分；认为激光能量沉积后形成的等离子体热核由N₂、O₂、N₂⁺、O₂⁺、NO、O、N、NO⁺、N⁺、O⁺和${{\rm{e}}^ - }$这11种组分组成。文中选取2001年Park提出的11组分模型，其中包含51个化学反应及计算化学反应速率所需的常数，具体的数值可参见参考文献[14]。

激光能量沉积后流场的演化过程满足质量守恒方程、动量守恒方程和能量守恒方程这三个基本方程，同时热核中还包含有各种组分，因此还需要满足描述传质过程的组分方程。完整的控制方程组如下所示。

连续方程

动量方程

能量方程

组分方程

其中

式中：S_m为质量源项，kg·m⁻³·s⁻¹；${\tau _{ij}}$为应力张量，Pa；H为单位质量气体的总焓，J·kg⁻¹；${J_{\alpha j}}$为组分α在j方向上的质扩散通量，kg·m⁻²·s⁻¹；S_h为化学反应热，W·m⁻³；${M_\alpha }$为组分α的相对分子质量，kg·mol⁻¹；${\omega _\alpha }$组分α的质量生成率，kg·m⁻³·s⁻¹；${h_{\alpha 0}}$为组分α的标准状态焓，J·kg⁻¹；${T_{\rm ref}}$为参考温度，取298.15 K；${c_{p_\alpha }}$为组分α的比热容，J·kg⁻¹·K⁻¹；${D_\alpha }$为组分α的质扩散系数，m²·s⁻¹。

上述方程组尚未达到封闭条件，还需要补充温度、压力和密度的约束方程。对于文中的高温高压情况，可以采用理想气体假设。因此：

式中：${M_{\rm w}}$为等离子体的相对分子质量，kg·mol⁻¹；R为气体常数，8.3145 J·mol⁻¹·K⁻¹。

激波和等离体热核是激光能量沉积的两个关键因素，需要在这两方面将数值计算结果与实验结果进行对比，以验证数值计算的可靠性。

数值计算中沉积的激光能量为10.1 mJ，与之相应，并考虑到等离子体屏蔽效应造成的能量沉积减小，在验证实验中采用脉冲Nd:YAG激光器在大气中沉积10.6 mJ的激光能量，聚焦透镜焦距15 cm，激光脉宽8 ns，波长1064 nm，利用美国Coherent公司的FieldMaxII-TOP型激光能量计监测入射激光能量。采用纹影测量技术针对激波半径和热核演化过程特性进行观测，包括脉冲氙灯闪光源、150 mm口径的纹影镜组、刀口和德国PCO公司的HSFC PRO超高速增强型相机。时序同步控制器选用美国Stanford Research公司的DG645数字脉冲信号发生器。

图3给出了数值计算得到的激波半径与实验值的对比。由图可知，实验测得的激波半径始终略大于计算值，这可能是因为实验沉积的激光能量略大于数值模拟。在激光能量沉积后的40 μs内，激波半径的计算值与实验值误差在2.5%以内，表明文中提出的数值计算方法可以很好地模拟激波的传播。

Figure 3.  Comparison between simulation and experiment of shock wave radius

图4是t=60 μs、t=90 μs和t=120 μs时热核的纹影图和计算得到的密度梯度云图对比，六幅图均为同等比例尺。需要注意的是，纹影拍摄是密度变化在光路上累加的效果，而计算得到的密度梯度云图仅为一个剖面上的结果。在t = 60 μs时，二者稍有差别，纹影图中尖刺开始萌芽，而计算结果中尖刺已经发展得较完备，这是因为初始温度剖面的不够完善，使得热核过膨胀程度过大或过膨胀过程发展过早。而t=90 μs和t=120 μs时的计算结果与实验差别很小，说明这种模拟方法在热核演化的后期十分可靠。

Figure 4.  Schlieren images and density gradient contours of hot core at t = 60 μs, t = 90 μs and t = 120 μs

利用激波管产生入射正激波，然后采用Fluent软件的UDF功能在正激波前方实现单脉冲激光能量沉积。

选择圆形激波管，高压段（0.40 m < x <0.60 m）和低压段（0 m < x < 0.40 m）长度分别为0.2 m和0.4 m。为了提高计算效率，网格的x轴方向在0.07 m < x <0.10 m区域加密（达到30网格/mm），而在其他段采用较稀疏的网格。各段在y轴方向的网格划分是完全一样的，网格沿着y轴负方向渐密，在x轴附近，y方向网格同样达到30网格/mm。对称轴x轴为axis边界（采用二维网格计算三维问题），其余边界均为无滑移壁面边界条件。

采用上述激波管产生马赫数1.92的入射正激波。低压段为13800 Pa、163 K的空气，根据激波关系，可求出合适的高压段条件之一为：2.136×10⁵ Pa、300 K的空气。验证计算得到的入射正激波马赫数1.92，速度491.08 m/s。接触界面的速度为298.22 m/s。

能量沉积区域选在x=0.085 m，y=0 m处，落在网格加密区中央。所模拟的激光入射方向为x轴负方向，能量大小约10.1 mJ。

正激波从x=0.4 m处产生后沿x轴负方向传播，与等离子体热核相互影响需要经历0.315 m的路程。根据正激波和接触界面的速度差，可求得正激波扫过热核后，热核不受接触界面干扰的自由演化时间达到270 μs，完全满足此节研究的要求。

文中计算了激光能量沉积8 μs （t=8 μs）后，热核与正激波相遇的情形。

由图5可知，在t=7 μs时，入射正激波（I）已经与激光能量沉积诱导的激波（Ⅱ）相交，但尚未与热核（Ⅲ）接触。正激波波后的一部分气体受到激波的再次压缩，压力由5.7×10⁴ Pa上升至8.3×10⁴ Pa。后续类似云图的横纵坐标单位均为m。

Figure 5.  Contours of density gradient and pressure at t = 7 μs

t = 22 μs时，正激波（I）刚刚完全穿透热核（Ⅲ），热核正在逐渐演变为涡环，如图6所示。密度梯度云图中仍能观察到热核受压缩而形成的微弱的稀疏波（Ⅳ）在传播。激波（Ⅱ）已经不再是连续曲面，在正激波左右已经割裂，这是由于激波在正激波前后的传播环境不一致所导致。正激波尚在热核内部时，称为透射激波（Transmitted shock wave），而当正激波完全穿透热核后，称为再透射激波（Ⅴ，Re-transmitted shock wave），再透射激波并不是平面，而是凸向正激波运动方向。透射激波转变为再透射激波之前，会与热核内表面相碰撞，从而产生再透射-反射波系（Ⅵ，Re-transmitted reflected shock waves）。由于热核内表面是近似对称的曲面，再透射-反射波系会形成交叉结构。再透射激波与热核之间会产生马赫杆（VII）。正激波、再透射激波和马赫杆相交于三波点（VIII），从三波点延伸出滑移层（Ⅸ）。结合压力云图可知，高压区主要有两个，一个是激波波后，压力约为正激波波后压力p₂的114%；另一个则为再透射-反射波波后，压力约为p₂的123%。低压区同样也有两个，一个是稀疏波波后，压力约为P₂的80%；另一个为热核演化形成的涡环位置，压力约为P₂的70%。

Figure 6.  Contours of density gradient and pressure at t = 22 μs

后续的计算结果表明，只有涡环所在的低压区是能够比较长时间地存在的。在文中算例中，涡环演化100 μs后涡环中心压力仍保持在p₂的70%左右。随着稀疏波的衰减，对应的低压区压力逐渐回升；随着激波、再透射-反射波系的衰减，对应的高压区压力逐渐降低。

此节选取的算例为Sasoh开展的实验^[15]。Sasoh利用吸气式风洞产生1.92 Ma的实验段气流，冲击直径20 mm的圆柱形钝头体。在圆柱形钝头体上游21.2 mm处沉积单脉冲激光能量，入射激光能量为12.0 mJ。Sasoh采用纹影系统观察流场演化，并设计力天平测量波阻。此节数值模拟中拟采用的激光能量沉积率为0.167，即沉积的激光能量为2 mJ。计算中所需的参数如表1所示。

所采用的计算网格和边界条件如图7所示。来流的方向为由左至右，设置为压力远场边界。右侧设置为压力出口，x轴上为axis边界（采用二维网格计算三维问题），圆柱形钝头体表面为无滑移壁面边界，采用C型网格。在axis边界均匀布置900个网格；压力出口边界作为与其相对应的边，同样设置900网格，并且在靠近钝头体的一端进行了加密处理。钝头体圆柱的高取为0.01 m。圆柱底面的半径和高分别均匀布置300网格；压力远场边界作为与二者相对应的边，设置为600网格，且在两端加密。这样设置可以保证热核与钝头体流场相互影响区域的网格密度达到30 网格/mm，经计算验证，满足网格无关性条件。

Figure 7.  Boundary conditions and grid division

在无激光能量沉积时，文中计算得到的基准波阻D_B为21.1 N，Sasoh实验测得基准波阻为（22.4±0.3） N，误差5.8%。相比于Sasoh的无粘、不考虑真实气体效应的数值结果20.6 N （误差8.0%），文中钝头体基准波阻的计算更准确。文中计算得到激波脱体距离为9.2 mm，实验值为9.0 mm，误差2.2%。同时，钝头体左端面边缘处的膨胀波以及由边缘延伸出的滑移层也与实验相当吻合，如图8所示。基准波阻、激波脱体距离、膨胀波以及滑移层的计算结果验证了上述数值方法的可靠性。

Figure 8.  Comparison of density gradient contour and schlieren image without laser energy deposition

钝头体波阻变化曲线如图9所示。图中波阻由钝头体表面压力积分而得，并进行了归一化处理。t = 0 μs时，激光能量沉积到流场中。t = 33 μs，波阻第一次开始上升；t = 36.5 μs时，波阻达到最大值，约为D_B （基准波阻，即无激光能量沉积时钝头体表面压力积分）的108.5%，此后波阻开始下降；t = 47.3 μs时波阻达到极小值，约为D_B的90.4%。t = 53.2 μs时波阻达到极大值，此后波阻逐渐降低至D_B的68.9%。t = 200.5 μs以后，波阻逐渐回升，并在t = 252.0 μs时达到极大值，约为D_B的107.6%。此后波阻逐渐下降，最终在D_B附近振荡，收敛于D_B。

Figure 9.  Curves of normalized wave drag with single pules laser energy deposition

为分析波阻变化的原因，下面将抽取波阻曲线上各个峰和谷对应时刻的流场图像进行研究。见图10~图14。

Figure 10.  Contours of density gradient and pressure at t = 36.5 μs

Figure 11.  Contours of density gradient and pressure at t = 47.3 μs

Figure 12.  Contours of density gradient and pressure at t = 53.2 μs

Figure 13.  Contours of density gradient and pressure at t = 150 μs

Figure 14.  Contours of density gradient and pressure at t = 220 μs

激光能量沉积诱导的激波会使钝头体波阻增大，导致波阻曲线出现第一个上升沿。图10给出了这种情况下的流场密度梯度云图和压力云图。图中可见激波（Ⅱ）和激波的反射波（XI）。事实上，激波到达钝头体左端面的瞬间正是钝头体波阻开始上升的时刻，为t = 33 μs。由于激波是球形，靠近轴线部分的最先反射。钝头体左端面附近的气体经历了激波和激波的反射波两次压缩，从而使压力明显升高。随着激波往下游传播，激波的反射波覆盖钝头体左端面的范围越来越大。在t = 36.5 μs时，激波的反射波已经全部覆盖钝头体左端面，此刻钝头体波阻达到第一个极大值。

再透射-反射波系（见图11中的Ⅵ）使钝头体波阻再次增大，波阻曲线出现第二个上升沿。图11给出了波阻曲线第二个上升沿开始时刻的流场密度梯度和压力云图。再透射-反射波系正好传播到钝头体左端面，并逐渐反射（类似于图10中激波的反射），产生再透射-反射波系的反射波（见图12中的Ⅻ）。当再透射-反射波系的反射波（Ⅻ）完全覆盖钝头体左端面时（图12所示情形），波阻出现第二个极大值。

在钝头体波阻达到第二个极大值后，由于再透射-反射波系的反射波离开钝头体，原本被包围在钝头体左端面和反射波之间的高压气体开始膨胀，从而导致钝头体波阻小幅降低。随后，具有低压力特性的涡环（图13中的Ⅲ）不断接近钝头体左端面，使得钝头体波阻进一步降低，表现为t = 53.2~200.5 μs之间波阻的大幅减小。如图13所示，t=150 μs，波阻降低至D_B的93%；而在t=200.5 μs时，钝头体波阻达到最小值，仅为D_B的68.9%，此刻低压涡环恰好到达钝头体左端面的边缘。随后涡环将离开钝头体左端面，而波阻也逐渐上升。经上述分析可知，热核经过弓形激波后，由于Richtmyer-Meshkov不稳定性而形成的低压涡环，是减阻的最主要机理。

当低压涡环离开钝头体左端面后，钝头体左端面会经历一个再压缩过程，导致波阻回升。这个过程会在钝头体端面的边缘首先形成压缩波（图14中的XⅢ），压缩波不断向弓形激波传播，钝头体左端面压力逐渐回升。在t=252 μs时，压缩波达到弓形激波，波阻达到最大值，约为D_B的107.6%。此后流场逐渐恢复到未受扰动的状态，波阻恢复到D_B。再压缩过程会导致波阻超调，如果抑制这个再压缩过程，那么减阻效果将会持续，这是进一步研究高重频激光能量沉积的突破点。

下面分析减阻百分比η和能量效率S。Sasoh在实验中不仅研究了单脉冲激光能量沉积减阻，还研究了高重频激光能量沉积减阻。高重频激光频率为4 kHz，每个脉冲的入射激光能量和单脉冲情形一样，均为12 mJ/pulse。高重频情形下实验测得的波阻降低为3%。数值模拟发现，单脉冲激光能量沉积的阻力变化曲线周期正好为250 μs，即文中图9的t=33~283 μs。因此，在4 kHz的高重频激光能量沉积作用下，钝头体的阻力变化曲线即是图9的t=33~283 μs时间段曲线的不断重复，这个结论在Sasoh的论文中也有所体现。将该段曲线积分可得波阻冲量I：

因此图9的t=33~283 μs时间段波阻的时间平均值等于

可得减阻百分比η （由于能量沉积而减小的波阻与基准波阻之比）：

文中数值模拟得到的减阻百分比（3.28%）稍大于实验测得的3%，表明数值模拟中沉积的单脉冲能量（2 mJ）大于实际值。这说明在Sasoh的实验工况下，当入射激光能量为12 mJ时，激光能量沉积率小于16.7%。而Sasoh的基于完全气体模型、求解欧拉方程的数值模拟方法，需要将能量沉积率提升至45%，才能得到与实验相符的减阻效果。这表明基于球形能量沉积模型和完全气体模型的数值模拟方法对于减阻效果的计算有较大误差，而文中的基于泪滴形能量沉积模型和真实气体效应的数值模拟方法具有较高的可靠性。

（1）单脉冲激光能量沉积与正激波的相互影响时，在正激波波后流场中，激波会产生一个高压区，稀疏波会产生一个低压区，再透射-反射波系会产生一个高压区，热核受正激波影响而形成的涡环会形成一个低压区，并伴随一个高压区。只有涡环所在的低压区以及它伴随的高压区压力是稳定的，并且低压区占据主导地位。

（2）单脉冲激光能量沉积与钝头体超声速流场相互作用时，激波和再透射-反射波系与钝头体碰撞后产生的反射波会使钝头体阻力增大，但由于作用时间短，对减阻效果的削弱作用较小；热核离开钝头体后流场再压缩过程是削弱减阻效果的最主要因素。热核穿过弓形激波后形成的低压涡环导致钝头体表面压力急剧减小，是波阻减小的主要机理。