激波诱导矢量喷管动态数值模拟

马伟，杜刚，金捷，廖华琳

(1.北京航空航天大学能源与动力学院，北京100191；2.中国燃气涡轮研究院，四川成都610500)

激波诱导矢量喷管动态数值模拟

马伟1，杜刚1，金捷1，廖华琳2

(1.北京航空航天大学能源与动力学院，北京100191；2.中国燃气涡轮研究院，四川成都610500)

对一二元收扩激波诱导矢量喷管进行了二维动态数值模拟，研究了不同落压比、不同次流加速时间下，喷管气动特性和流场结构随次流增加的变化规律。结果表明：喷管在过膨胀状态下，次流流量达到一定值时，其推力会有较强的振荡；在完全膨胀和欠膨胀状态下，振荡较小。其原因是，次流增加过程中，喷管出口附近有复杂波系存在和回流产生，而回流对提高推力和增加矢量角有益。

航空发动机；射流推力矢量喷管；激波诱导；动态；数值模拟

1 引言

射流推力矢量喷管由于结构简单、重量轻，及可降低飞机可探测性等诸多优点，在飞行器推进技术发展中极具前景，国内外对此都开展了大量研究。射流推力矢量喷管产生推力矢量，主要基于激波矢量控制法、喉道偏移法和反流法三种[1-2]。激波矢量控制法是在喷管扩张段内一侧注入次流，与主流作用产生斜激波，从而使得主流方向改变获得矢量推力，用该方法能获得较大的矢量角。国内外主要通过试验和数值模拟两种手段，以矢量角、矢量效率、推力系数为喷管性能指标，研究各参数变化对其影响及相关机理[3-10]。数值模拟方面，大多是对喷管进行定常模拟，较少研究动态变化参数对喷管性能的影响[11]。吴雄等[12]对固体火箭发动机的启动和关闭的非定常特性做了研究，指出动态与稳态特性相差很大，瞬态侧向控制力、燃烧室压力、推力等都会显著变化，因此对于燃气二次喷射系统的非定常特性研究，将成为飞行器控制系统设计的重点。Saha等[13]利用FLUENT软件，对激波管驱动的超声速喷管的启动过程进行了时间相关的数值模拟，并与实验值比较；计算捕捉到了流场的主要特征，表明利用商用软件可实现对瞬态激波相互作用流场的计算。

本文以一二元收扩喷管为研究对象，利用FLU⁃ENT软件，对次流质量流量从零增加到最大值这一过程进行数值模拟，研究不同落压比、不同次流加速时间下喷管特性的变化规律，详细分析该过程中流场结构变化，以期为射流推力矢量喷管研制提供参考。

2 计算方法及算例验证

2.1 计算格式

通过时间推进的有限体积法求解二维非定常雷诺平均N-S方程，湍流模型为RNGk-ε模型。近壁采用标准壁面函数，空间离散采用二阶迎风格式，时间离散采用二阶隐式格式。

2.2 几何模型与网格划分

采用的喷管模型为NASA兰利研究中心的试验模型，设计落压比NPRD=8.78，详细参数见文献[14]。计算域如图1所示，长度和宽度分别为喷管喉道高度的23、28倍。计算网格在壁面、次流喷缝、喷管喉道等物理量变化梯度较大位置加密处理。

图1 模型计算域Fig.1 Model computational domain

2.3 边界条件

喷管的主流入口设定为压力入口条件，给定总压、总温；计算域上、下、后边界都设为压力出口条件，给定总压(为环境大气压)和总温；喷管壁面及计算域前边界设为绝热、无滑移壁面条件。次流入口设为质量入口条件，使用FLUENT中的UDF功能实现对次流的质量流量控制。次流质量流量随时间的变化函数为：

2.4 算例验证

为验证数值模拟的准确性，对网格数和时间步长进行了敏感性分析，并与文献[14]中的试验数据做了对比。

在NPR=4.60、次流总压与主流总压之比SPR= 0.7，对网格数为10万、20万、40万的网格进行定常数值模拟；对NPR=4.60，ts=1.00 s，时间步长为0.01 s、0.005 s、0.002 5 s进行动态数值模拟。从图2和图3中结果可得出，数值模拟结果与所选取的网格数和时间步长无关，并与试验符合较好。下面模拟均选择网格数为20万的网格，时间步长0.005 s。

图2 喷管上、下壁面静压数值结果与试验结果对比Fig.2 Computational and experimental static pressure on nozzle wall

3 计算及结果分析

对NPR=4.60、8.78、10.00，ts=0.25 s、0.50 s、1.00 s共9种工况进行计算。

3.1 喷管参数计算

无次流情况下，理想推力计算公式为：

式中：pt为喷管入口气体总压，Ath为喷管喉道面积，pamb为环境压力。

有次流情况下，理想推力为主流产生的推力Fpri与次流产生的推力Finj之和。Fpri可直接用式(2)求出，Finj由式(3)计算。

图3 矢量推力随时间的变化Fig.3 Vector thrust vs.different time-steps size

式中：pt,inj,ave为次流进口质量加权平均总压，且分别为次流入口边界网格单元面上的质量流量和总压，k、R为空气的物性参数。

图4 实际推力变化曲线Fig.4 Real thrust

图5 推力系数变化曲线Fig.5 Thrust ratio

3.2 喷管特性分析

3.2.1 推力及推力系数

由于m˙inj是流动时间t的函数，由式(1)可知，当无量纲时间t/ts相等时，不同ts下的也相等。从图4中可看出，NPR越大推力越大，不同NPR下推力随t/ts的增大而增大。图4(a)中t/ts=0.85附近，推力有明显振荡；图5(a)中，随着t/ts的增大，推力系数先增后减，在振荡后趋于平稳。推力系数刚开始增加是因为NPR=4.60时，喷管处于过膨胀状态，初始流量小，流动损失小，但次流可改善过膨胀状态使推力系数增加，而后t/ts继续变大，流动损失增大，推力系数下降。图5(b)、图5(c)中，喷管都不处于过膨胀状态，次流的注入会引起推力损失，推力系数基本上随t/ts的增大而减小。

3.2.2 矢量角及矢量效率

图6 矢量角变化曲线Fig.6 Vector angle

图7 矢量效率变化曲线Fig.7 Vector efficiency

3.3 喷管流动动态过程分析

3.3.1 推力系数及推力随t/ts的变化

如图4、图5所示，NPR=4.60时，随着t/ts的增大推力系数先增后减，推力则单调增加；但在t/ts= 0.80～1.00时，推力和推力系数均有明显振荡。下面分析振荡原因。

图8示出了图4(a)中ts=1.00 s、t/ts=0.80～1.00这段内推力各分量的变化曲线。从图8(d)中容易看出，t/ts=0.82、0.85、0.86、0.90四个时刻相差很大，t/ts= 0.85是一振荡峰值。图9为这四个时刻下的马赫数云图，可见喷管出口附近刚好有一复杂波系，且随着m˙inj的增加，该波系逐渐由喷管出口外移动到喷管内。t/ts从0.80增大至0.85的过程中，该波系部分跨越喷管出口截面，动量项减小较快，但压力项增加得更快(图8(a)、图8(b))，结果是轴向推力突升，在t/ts=0.85处达到峰值(图8(d))。t/ts=0.85～0.86这一段，该激波系中正激波部分已完全进入喷管内(图9(c))，动量项减小很快(图8(a))，而此时压力项增加有限，整体表现为推力又迅速减小。t/ts=0.86～0.90这一段，该波系则进一步进入喷管内。通过分析可得出：在NPR=4.60时，喷管出口附近有一复杂波系，该波系从喷管出口截面外侧进入喷管出口截面内侧，使推力产生振荡。

NPR=8.78和10.00(图5(b)、图5(c))时，推力系数基本随t/ts的增大而减小，但在t/ts=4.0～5.5这一段出现了小波动，且ts=1.00 s波动最大，ts=0.50 s的次之，ts= 0.25 s的波动则很弱。为探究其原因，下面对NPR=8.78中ts=1.00 s和ts=0.25 s的流动过程进行分析。

图8 NPR=4.60、ts=1.00 s时实际推力各分量的变化曲线Fig.8 Force components of real thrust,withNPR=4.6，ts=1.00 s

图9 NPR=4.60、ts=1.00 s时的马赫数云图Fig.9 Computational mach contours，with NPR=4.60，ts=1.00 s

图10为推力各分量及理论推力的变化曲线。可见：ts=1.00 s和ts=0.25 s的理论推力变化曲线几乎重合(图10(f))，都随t/ts的增大而线性增大，且在t/ts= 0.40～0.50这一段内，两种ts下的理论推力增量都超过了实际推力增量(图10(e)、图10(f))。t/ts=0.42～0.44这一段内，ts=1.00 s的喷管实际推力各分量斜率出现较大变化，这就是出现波动的原因(图5(b))。t/ts=0.44～0.48这一段内，ts=0.25 s的喷管实际推力各分量斜率出现较大变化，推力系数下降变慢。图11、图12分别为ts=1.00 s和ts=0.25 s下三个时刻的流场图。图10中，t/ts=0.42～0.44这一段内，ts=1.00 s时实际推力增速变快。图11(a)中t/ts=0.42时刻，喷缝后主流再次附着壁面，形成的角涡区封闭，此时无回流；而图11(b)中显示主流与壁面分离，角涡区不再封闭，出现回流；图11(c)中的回流区更大。结合图10、图11可得出，图10中实际推力增长斜率从t/ts=0.42开始突然变大，是因为喷缝后主流与壁面分离，出现回流，喷管出口外面压力高的气体进入喷管内，改善了喷管的过膨胀状态，增加了推力。这就使得图5(b)中ts=1.00 s时的推力系数出现小的波动。ts=0.25 s时的流动也经历了这一过程(图12)，只是实际推力增长斜率变化不大，且回流出现在t/ts=0.44时刻，相比ts=1.00 s稍晚。NPR=10.00时的流动过程与NPR=8.78的类似。

图10 NPR=8.78时不同ts下的实际推力分量变化曲线Fig.10 Force components of real thrust,withNPR=8.78

3.3.2 矢量角和矢量效率随t/ts的变化

从图6中可知，矢量角随t/ts增加的变化规律与推力的变化规律基本相同，矢量角随t/ts先增加，在t/ts=1.00后趋于稳定。矢量效率随t/ts的变化较复杂，NPR=4.60时，如图7(a)所示，矢量效率先急剧变小而后急剧增大，再相对缓慢变小至稳定。NPR= 8.78时，如图7(b)所示，矢量效率先急剧变小而后急剧增大，再缓慢减小至平稳一段后变大，之后变小趋于稳定。NPR=10.00与NPR=8.78的类似，变化曲线整体呈阶梯状。

图7(b)、图7(c)中，矢量效率在t/ts=0.40～0.55区间内增大，随次流的增加出现了回流，而回流(图11、图12)改善了次流后面喷管的过膨胀状态，使得喷管上壁面压力升高，有助于增大矢量角。

3.3.3 ts对推力系数及推力的影响

图4中，NPR=4.60时三个ts值下的推力系数变化趋势完全一致；NPR=8.78和NPR=10.00时，三个ts值下的推力系数变化趋势基本一致，但在t/ts= 0.40～0.55，不同ts值下波动出现的时刻和强度不同。ts=1.00 s的波动出现得最早也最明显，ts=0.25 s的波动出现得最晚(这里早、晚是在无量纲时间下而言)也最弱。波动出现时刻不同的原因，是因为ts= 0.25 s情况下次流增加较快，主流来不及响应，所以分离出现较晚；同样可解释ts=1.00 s时分离出现较早的原因。ts=0.25 s时推力系数波动最为平坦，是因为相比ts=1.00 s，从回流产生至其影响最大需要更长时间，这点可从图10中看出，所以波动不明显。

图11 NPR=8.78、ts=1.00 s时的马赫数云图和流线图Fig.11 Computational Mach contours and streamlines,with NPR=8.78,ts=1.00 s

图12 NPR=8.78、ts=0.25 s时的马赫数云图和流线图Fig.12 Computational Mach contours and streamlines，withNPR=8.78，ts=0.25 s

3.3.4 ts对矢量角及矢量效率的影响

图6中，不同ts下矢量角随t/ts的变化趋势完全一致。以图6(b)为例，在t/ts=0.44时，不同ts下的矢量角相差最大。图13为t/ts=0.44时不同ts值下喷管上壁面的压力分布，及该m˙inj下定常流动时喷管上壁面的压力分布。可见，不同ts下，喷缝后面喷管上壁面的压力分布差异较大。图14为t/ts=0.44时不同ts下的流场图，及该m˙inj下定常流动时的流场图。可见，相同t/ts下，ts=0.25 s没有回流出现，ts=1.00 s回流区最大，这时的流动最接近该工况下的定常流动。这也是t/ts=0.44下ts=1.00 s的矢量角最大、ts=0.25 s的矢量角最小的原因。不同ts下，矢量效率随t/ts的变化与矢量角的完全一致。

图13 NPR=8.78、t/ts=0.44时喷管上壁瞬态压力及该次流流量下的定常静压分布Fig.13 Pressure distribution along upper wall at transient and steady condition,withNPR=8.78,t/ts=0.44

4 结论

(1)NPR=4.60时，喷管处于过膨胀状态，当次流流量增加到一定程度后，喷管推力随次流的增加而振荡。其原因是喷管出口附近存在一复杂波系，随着次流的增加由喷管出口外向喷管内移动，导致推力剧烈变化。

(2)NPR=8.78与NPR=10.00下，随着次流的增加，推力系数总体呈下降趋势，此过程中会出现一定波动，原因是次流流量达一定值后，喷缝后的主流与喷管壁面分离，出现回流，使得次流喷缝后喷管壁面的负压状态得到改善，减小了推力损失。

(3)回流的产生有利于提高推力系数、矢量角和矢量效率。

(4)在无量纲时间尺度下，次流加速时间越短，回流产生就相对越滞后。

图14 马赫数云图与流线图Fig.14 Computational Mach contours and streamlines

[1]Deere K A.Summary of Fluidic Thrust Vectoring Re⁃search Conducted at NASA Langley Research Center[R]. AIAA 2003-3800，2003.

[2]Mason M S，Crowther W J.Fluidic Thrust Vectoring of Low Observable Aircraft[C]//.CEAS Aerospace Aerody⁃namic Research Conference.Cambridge，2002.

[3]Deere K A.Computational Investigation of the Aero⁃dynamic Effects on Fluidic Thrust Vectoring[R].AIAA 2000-3598，2000.

[4]Hamed A，Laskowski G.A Parametric Study of Slot In⁃jection Thrust Vectoring in a 2DCD Nozzle[R].AIAA 1997-3154，1997.

[5]Neely A J，Gesto F N，Young J.Performance Studies of Shock Vector Control Fluidic Thrust Vectoring[R].AIAA 2007-5086，2007.

[6]Williams R G，Vittal B R.Fluidic Thrust Vectoring and Throat Control Exhaust Nozzle[R].AIAA 2002-4060，2002.

[7]张群锋，吕志咏，金捷，等.轴对称射流矢量喷管的试验和数值模拟[J].推进技术2004，25(2)：139—143.

[8]王庆伟，张相毅，徐学邈，等.射流角度对双缝射流喷管流场影响的研究[J].燃气涡轮试验与研究，2006，19(4)：27—29.

[9]张强，杨永，李喜乐.双喉道射流推力矢量喷管的数值模拟研究[J].西北工业大学学报，2009，27(6)：754—759.

[10]秦亚欣，于军力，高歌.激波诱导圆形矢量喷管数值研究[J].航空动力学报，2009，24(10)：2208—2212.

[11]宋亚飞，高峰，马岑睿，等.激波诱导矢量喷管内流场动态特性数值研究[J].系统仿真学报，2013，25(1)：36—40.

[12]吴雄，焦绍球，张为华，等.燃气二次喷射对推力矢量控制发动机内流场的影响[J].南京航空航天大学学报，2007，39(6)：775—780.

[13]Saha S，Chakraborty D.Numerical Exploration of Starting Process in Supersonic Nozzle[J].Aeronautical Journal，2007，111(1115)：51—58.

[14]Waithe K A，Deere K A.Experimental and Computational Investigation for Multiple Injection Ports in a Conver⁃gent-Divergent nozzle for Fluidic Thrust Vectoring[R]. AIAA 2003-3802，2003.

[15]Anderson C J，Giuliano V J，Wing D J.Investigation of Hy⁃brid Fluidic/Mechanical Thrust Vectoring for Fixed-Exit Exhaust Nozzles[R].AIAA 1997-3148，1997.

Dynamic Numerical Simulation of Shock Vector Control Exhaust Nozzle

MA Wei1，DU Gang1，JIN Jie1，LIAO Hua-lin2
(1.School of Jet Propulsion，Beijing University of Aeronautics and Astronautics，Beijing 100191，China；2.China Gas Turbine Establishment，Chengdu 610500，China)

2D dynamic numerical simulation of a shock vector controlled 2D-CD exhaust nozzle was con⁃ducted to investigate the effect of increasing secondary mass flow on the internal nozzle performance and flow field structure under different nozzle pressure ratios.Results indicate that when the secondary mass reaches a certain value,the thrust of nozzle will undergo a strong oscillation when nozzle is under over-ex⁃panded condition.The oscillation will be much weaker when nozzle is under either fully expanded or un⁃der-expanded conditions.This phenomenon is attributed to a series of complex waves near the nozzle outlet and the emergence of back flow.It is also concluded that back flow can enhance the performance of nozzle and enlarge the vector angle.

aero-engine；fluidic thrust vectoring control；shock vector control；dynamic；numerical simulation

V231.3

：A

：1672-2620(2014)05-0030-08

2013-12-18；

：2014-06-30

马伟(1990-)，男，江西九江人，硕士研究生，研究方向为发动机内流气动力学。