激波誘導矢量噴管動態數值模擬

馬偉，杜剛，金捷，廖華琳

(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-)，男，江西九江人，碩士研究生，研究方向為發動機內流氣動力學。