等离子体合成射流对钝头激波的控制与减阻

2021-08-03陈加政胡国暾樊国超陈伟芳

航空学报 2021年7期

陈加政，胡国暾，樊国超，陈伟芳,*

1.浙江大学航空航天学院，杭州 310027 2.中国运载火箭技术研究院研究发展中心，北京 100076

等离子体合成射流作为一种主动流动控制方法[1]，由于其激励器具有体积小、响应快和射流速度高等工作特点，可在不改变飞行器整体外形的情况下对其流场进行流动控制，从而改善飞行器的气动特性和可操控性，为飞行器的设计和优化提供冗余空间。其近年来在边界层控制[2-4]和高速流动控制领域[5-7]逐渐展现了良好的工程应用前景。

为了探究等离子体合成射流对流场的作用机理和控制效果，国内外大量学者对其进行了数值模拟研究。Grossman等[8]发展了一个瞬时加热模型，其认为激励器腔体内能量沉积过程瞬间完成，将激励器腔内气体设定较高初始温度和压力，对激励器腔体及其出口流场进行了数值模拟。Cybyk等[9-10]也采取了类似的思路，研究了激励器出口直径、腔体大小及能量沉积大小对“SparkJet”激励器性能的影响。国内的张宇[11]也利用了该方法，研究了两电极等离子体合成射流激励器的直接力工作特性。刘朋冲等[12]则根据分子平转动迅速平衡的特点，直接利用激励器内放电等离子体的转动温度为800 K的光谱结果，将其设为仿真的温度边界条件，通过控制持续时间来调整能量注入大小，研究了单次放电激励器流场的演化规律以及不同放电时间尺度对激励器出口速度的影响。

由于Grossman的瞬时加热模型不能描述能量沉积过程，仿真的初始条件严重依赖于简化的理论分析模型的精确性，难以真实反映激励器的射流速度、流场结构、工作周期等特性，因此学界又进一步发展了相应的能量源项模型。单勇等[13]认为，能量沉积的平均效果可以通过在一定时间内注入不变的能量来表示，从而能在能量控制方程中添加源项来模拟放电后在腔体内产生的高能量密度环境，并在考虑固体外壁面与外界环境的辐射换热的条件下，数值模拟了等离子体合成射流激励器的流场，获得了流场的典型结构。Lv等[14]则在此基础上，进一步研究了激励器频率、壁面换热系数以及壁面温度对激励器性能的影响。顾仁勇等[15]也利用该方法对运输机后体流动分离进行了主动控制的数值研究。与此同时，国防科技大学的罗振兵团队[16-18]也通过能量源项模型研究了等离子体合成射流的一些特性。他们在腔体内能量注入的基础上，还考虑了等离子体气体的部分热力学特性。

在能量源项模型的基础上，本文进一步完善了等离子体在热完全气体效应下的热力学属性和输运参数，从而提高了数值模拟中激励器放电功率的上限与计算的稳定性，改善了能量源项模型的计算精度与适用范围。同时，通过对超声速平板流场的动态模拟，对比分析不同热力学参数下合成射流对主流的干扰作用以及流场结构的演化过程，探究等离子体合成射流对边界层流动发展的作用效果，同时验证等离子体物性参数拟合的必要性。通过对马赫数为3超声速球头流场的动态模拟，探究等离子体合成射流对球头流场的调控作用，并分析其对球头激波及其阻力的影响规律，为钝头体飞行器的减阻布局设计提供数值依据。

1 等离子体合成射流数值建模

1.1 能量源项模型

针对等离子体合成射流，考虑其物理机制主要为气体放电的电加热作用，因此可以将激励下复杂的物理化学流动简化为唯象的能量源项模型，把放电的能量沉积添加到控制方程的源项中，从而模拟等离子体射流的形成。其数值控制方程为

(1)

剪切应力τij的表达式为

(2)

式中：μ为动力黏性系数；V为速度矢量；δij为应力张量分量。

(3)

(4)

式中：Ein为激励器消耗总能量大小；ηe为能量转化效率；t为时间；n为周期数；T为周期时间；V为放电区域体积；τ为每个周期内输入能量的持续时间；Qav为放电总能量转化成气体热能的能量密度，W/m3。

1.2 等离子体的热力学属性与输运系数

当等离子体合成射流激励器发生放电作用后，腔体内的温度将迅速上升，空气发生电离形成具有较高温度的等离子体。此时，腔内气体将不再是量热完全气体，其热力学性质将发生改变。因此，为求解式(1)，必须先给定定压比热cp、热传导系数κ与黏性系数μ等参数，其随温度和压力的变化而变化。但腔体内等离子体成分方程复杂，而射流的非定常特性又使得直接采用化学非平衡动力学方程求解异常低效。因此本文将等离子体视为热完全气体(单一组分)，采用工程拟合公式来计算等离子体的物性参数，从而快速地求解等离子体控制方程式(1)，同时满足激励器不同放电功率下的计算精度要求。而对于等离子体物性参数的工程拟合公式，本文选取Capitelli[19]与Gupta[20]等的研究成果，其温度范围在50～30 000 K 左右，具体参数可参考文献[19-20]。

等离子体的定压比热为

cp=exp(aζ4+bζ3+cζ2+dζ+e)

(5)

等离子体的黏性系数为

(6)

等离子体的热传导系数为

(7)

式中：

ζ=ln(T/10 000)

(8)

(9)

γ(T,b,c)=e-[(T-b)/c]2

(10)

式(5)～式(7)中唯一自变量为热力学温度T，其他参数均为拟合参数。当温度小于 30 000 K 时，定压比热cp、焓值h、黏性系数μ和热传导系数κ随温度变化的曲线如图1所示。其中，cp_Constant表示定压比热为常数1 006.43 J/(kg·K)，适用于量热完全气体；cp_Gupta91适用于处于局部热力学平衡态的空气，为式(5)的拟合结果。h_Linear和h_Gupta91分别表示cp为常数和拟合积分结果对应的焓值。μ_Sutherland表示采用Sutherland公式计算的黏性系数；μ_Capitelli 表示式(6)拟合结果。κ_Sutherland表示黏性采用Sutherland公式计算对应的热传导系数；κ_Capitelli为式(7)的拟合结果；κ_Experiment为实验结果[19]。

图1 定压比热、焓值、黏性系数和热传导系数随温度的变化曲线

由图1可见，在2 000 K以下等离子体气体的定压比热、黏性系数和热传导系数与完全气体相差不大。但2 000 K之后随着温度的逐渐升高，等离子体气体的热力学属性与输运系数相比量热完全气体的误差逐渐增大，甚至达到1～2个数量级的差异。而激励器中在剧烈放电作用下温度可达约10 000 K，将其当作量热完全气体处理存在较大的误差。

2 合成射流流动控制

2.1 超声速平板边界层控制

算例采用三维结构网格数值仿真，对称面计算网格如图2所示，左端压力入口距喷口130 mm，右端压力出口距喷口100 mm，上端压力远场距喷口100 mm，展向宽度为100 mm。外部流场网格数目为221×101×101，总网格量(含激励器)为320万。喷管直径为3 mm，长5.5 mm，激励器为直径7.5 mm、高10 mm的圆柱体，放电区域为腔体中间1/4区域，放电时间为20 μs，注入能量为1.92 J。来流条件与文献[16]中实验一致，静温为166.7 K，静压为12 950 Pa，马赫数为2。壁温取来流静温。计算时间步长为1×10-8s，内迭代步数为20。湍流模型采用RNGk-ε模型[21]。

图2 平板流场计算域及其网格

图3展示了激励器放电能量效率为10%时采用式(5)～式(7)计算的结果。为了方便与文献[16]的实验纹影图结果进行对比，采用密度梯度灰度图显示，其中A系列代表计算结果，B系列代表实验纹影图，H为等离子体合成射流诱导的弓形激波最大高度。由图可见，对于超声速平板流场的基本结构，数值模拟结果在流场大尺度结构上与实验较为符合，但在时间上存在约10 μs的滞后。这主要是由于计算中能量是在放电阶段均匀注入，而真实激励器工作电压在放电前半阶段较大并逐渐变小，大部分能量在前半段注入，从而造成了约为放电时间一半即10 μs的时滞。对于数值模拟结果，在35 μs左右，已可以观察到明显的由合成射流引起的弓形激波，并开始向下游发展。约从60 μs开始，大尺度涡开始形成，并由于主流作用向下游传播。150 μs左右，扩大的弓形激波转变为斜激波，并在200 μs时强度达到最大。此后斜激波逐渐减弱，且涡结构也逐渐耗散消失。

图3 合成射流作用下平板流场结构演化实验[16]与数值模拟结果比较

为了更好地分析能量效率和定压比热的选取对合成射流计算结果的影响，图4给出了不同能量效率下等离子体合成射流诱导的弓形激波最大高度H随时间的变化。

图4中百分数表示能量效率，虚线表示定压比热采用常数形式，取值为1 006.43 J/(kg·K)，而实线表示定压比热采用拟合公式的计算结果，Exp为文献[16]中的实验数据。由图4可见，对于常数定压比热的情况，其弓形激波最大高度的时间演化规律较为一致，能量效率越高，其值越大，曲线的斜率也有所上升，5%能量效率的结果与实验吻合较好，但曲线斜率相较实验结果更大。而采用拟合公式定压比热的结果，10%能量效率的结果与实验吻合很好，尽管能量效率选取较大，但弓形激波最大高度却上升不多，其斜率小于常数定压比热的曲线。这是因为拟合公式获得的定压比热在等离子体温度上升到较高值时比常数值要大得多，气体的焓值高,吸热能力上升了。对于相同的能量密度，拟合公式的结果使气体内能增加而动能减少，流场结构的发展更加缓慢，表现为弓形激波高度和增长斜率都会减小，从而得到与实验更为符合的结果。

图4 弓形激波最大高度时间演化曲线

2.2 钝头激波控制与减阻

算例采用轴对称仿真，其计算区域如图5所示，壁面网格加密，总网格量为26万。来流静温107.1 K，静压2 758 Pa，马赫数为3，单位雷诺数为7.5×106。钝头、腔体和喷管壁设为无滑移等温壁，壁温取来流静温。钝头直径50 mm，喷管直径9 mm，长4.5 mm，激励器为直径12.4 mm、高12.4 mm的圆柱体，放电区域依旧为腔体中间1/4区域，放电时间为15 μs。稳态条件下(射流未启动)数值模拟结果显示激波脱体距离为5.5 mm，总阻力(包括压差阻力和黏性阻力)大小为35.0 N。

图5 球头流场计算域及其网格

图6给出了实际注入热能Qin=Ein·ηe=270 mJ的合成射流作用下流场结构密度梯度灰度图。由图6可知，在24 μs左右，射流开始在钝头头部出现。40 μs时刻在激励作用下，头部激波的顶端被推向上游并形成一个凸起。60 μs和80 μs，头部激波进一步被推向上游，高温低密度射流在头部激波和钝头体之间扩散，并有部分开始沿钝头体表面向下游流动。在100 μs左右，射流的强度达到峰值，头部激波的前移距离达到最大。此时由于射流的速度较大，密度较低，为了平衡球头弓形激波后的低速高密度特性，在射流内部产生了一道正激波。之后，射流变得越来越弱，空气开始重新填充激励器腔体。在280 μs左右，射流基本喷射完成，将开始下一个周期的运动。

图6 合成射流作用下球头流场结构随时间的演化

图7给出了在前2个周期内，钝头流场弓形激波的脱体距离随时间演化的情况。其中Steady表示无射流定常状态下的激波脱体距离，Transient表示激励器工作时的结果。由图可见激波脱体距离随时间呈明显周期性变化。对于第1个周期，激波脱体距离在125 μs时达到最大,为14.86 mm;而在290 μs达到最小,为3.77 mm。相较于未施加射流时5.5 mm的情况，激波脱体距离最大增加了2倍，由此可见在射流周期的喷阶段，等离子体合成射流对头部激波有较强的控制作用。

图7 激波脱体距离随时间的变化

由图8可知，等离子体合成射流的喷口(Exit)平均速度(V)和钝头阻力(D)同样呈明显的周期性变化，且阻力变化相较速度有迟滞现象。在主射流结束后，激励器会出现二次射流，但其强度较主射流有大幅下降，阻力和射流速度变化幅度都明显减小。在第1个周期内即主射流作用下，相较于无射流定常状态，平均阻力降低了6.3%，而阻力最大则降低了32.0%。可见合成射流有较强的减阻效应。

图8 喷口平均速度与阻力随时间的变化

为了更好地解释球头阻力特性的演化规律，图9给出了0～42 μs激励器腔体的压力变化云图。由图9可知，对于图8中阻力演化初期阶段的2个跳跃，是由激励器内气体运动造成的。在阻力的极大值处，由10～20 μs的压力云图可见，等离子体放电产生的高温高压区形成正激波向两侧传播，传播到激励器两端使激励器壁面压力骤然上升，而左侧肩部壁面面积较底部小，使得合力为阻力，从而使阻力有剧烈上升但持续时间很短；在极小值处，由20～42 μs的压力云图可见，在喷口的拐角处由于膨胀而产生了2个低压区，同时入口肩部气流由于低压抽吸以及喷射作用下的气体流失也产生了低压区，并与出口处相持。但与此同时激波在底部反射离开而使底部压力下降，阻力减小。之后反射波到达肩部使其压力迅速上升，从而使阻力急剧减小。50 μs后激励器内部激波消散，阻力变化主要由激励器外部射流和主流作用引起。

图9 射流早期激励器腔体与出口压力云图

图10则说明了阻力在射流各阶段变化的原因。对于射流的演化而言，其喷阶段弓形激波外推，波后的压力总体降低约10 kPa，并且射流形成了涡对低压区，2种共同作用使阻力减少。而吸阶段激波脱体距离逐渐减小，波后压力随之升高，腔体内形成高压区，从而使阻力增大，且阻力相比速度周期具有滞后性。滞后性则是因为在刚开始的转换阶段，比如喷转吸，低压区不会马上消失，而是慢慢减小，故阻力不会立刻转换。

图10 典型时刻流场压力云图

图11说明尽管能量会被耗散或传递至下游，但还是有相当一部分高温高压气体被二次吸入致使腔体增压，在腔体内形成高压区，进而产生强度较小的二次射流，又使脱体距离增加，形成周期性振荡。

图11 典型时刻流场焓值云图

为使计算结果更具普遍性和参考，下面给出了实际放电气体热能Qin注入为108、162、216、270、324 mJ的数值模拟结果对比。由图12、图13 可以看出，整体上不同热能注入的计算结果阻力和射流速度变化规律基本相同，分为喷阶段和吸阶段。热能的大小主要影响喷阶段，随热能的增大，射流速度的最大速度和能量增大导致激波外推距离变大，波后压力减小使减阻效应更显著。而吸阶段阻力和回填速度变化不大，其与注入热能的大小几乎无关。并且喷吸各阶段持续时间不随热能大小变化而变化，整个作用周期为固定值。

图12 不同热能下出口速度变化

图13 不同热能下阻力变化

等离子体合成射流控制钝头激波的最终目的是减小钝头飞行器的能量消耗，减阻所节约的能量需大于激励器所消耗的能量。

激励器消耗的总能量为

(11)

式中：Ein为电容能量；Cd为电容大小；Vb为击穿电压。击穿电压需实验测量，其与两极距离、介质等有关。

激励器所消耗的总能量存在电路的耗散和其他的能量损失，因此放电实际注入热能为

(12)

式中：Td为放电时间；U(t)为瞬态电压；I(t)为瞬态电流。

一个放电周期减阻收益的能量为

(13)

式中：D0为定常状态阻力；D(t)为瞬态阻力；U∞为自由来流速度。表1给出了不同能量注入的减阻收益对比,表中ηb为减阻收益大于消耗总能量的最小能量效率。

由表1可以看出，减小的阻力所节约的能量是远大于放电实际注入的热能。但由于放电电路的耗散和其他的能量损失，激励器消耗的总能量不能全部转化为所注入的热能。钝头激励器由于弓形激波的压缩作用,其激励器腔体内压力为33 kPa,较高的环境压强能明显提升激励器的工作性能[22]。并且文献[16, 23]中提出的三电极激励器，其通过对电极间隙、电容、电压等激励器参数选取优化，放电效率最高可达80%以上，即能实现减小整个系统能量消耗的目的。

表1 减阻收益对比

由于等离子体合成射流单次放电的控制作用时间短，在工程应用中一般以一定频率多次脉冲放电的形式为主。而等离子体激励器单个工作周期为能量沉积、射流喷出和吸气复原3个阶段，下个放电时刻应在激励器吸气复原之后且避开二次射流射出阶段，否则激励器工作效果不佳。本文选取复原刚好完成时刻进行下次放电，周期为365 μs，该放电频率为饱和频率2 740 Hz。由图14 可以看出，各周期的射流速度和阻力变化规律与单周期基本相同，第1个周期与后续周期略有不同，后续周期射流峰值较大。图示4个放电周期的阻力平均值相比定常减小6.1%，与单周期6.3%相差不大。这说明当放电频率小于饱和频率，多放电周期可以看作单周期的叠加。