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

(1.北京航空航天大學能源與動力學院，北京100191；2.中國燃氣渦輪研究院，四川成都610500)

激波誘導矢量噴管動態(tài)數(shù)值模擬

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

(1.北京航空航天大學能源與動力學院，北京100191；2.中國燃氣渦輪研究院，四川成都610500)

對一二元收擴激波誘導矢量噴管進行了二維動態(tài)數(shù)值模擬，研究了不同落壓比、不同次流加速時間下，噴管氣動特性和流場結(jié)構(gòu)隨次流增加的變化規(guī)律。結(jié)果表明：噴管在過膨脹狀態(tài)下，次流流量達到一定值時，其推力會有較強的振蕩；在完全膨脹和欠膨脹狀態(tài)下，振蕩較小。其原因是，次流增加過程中，噴管出口附近有復雜波系存在和回流產(chǎn)生，而回流對提高推力和增加矢量角有益。

航空發(fā)動機；射流推力矢量噴管；激波誘導；動態(tài)；數(shù)值模擬

1 引言

射流推力矢量噴管由于結(jié)構(gòu)簡單、重量輕，及可降低飛機可探測性等諸多優(yōu)點，在飛行器推進技術(shù)發(fā)展中極具前景，國內(nèi)外對此都開展了大量研究。射流推力矢量噴管產(chǎn)生推力矢量，主要基于激波矢量控制法、喉道偏移法和反流法三種[1-2]。激波矢量控制法是在噴管擴張段內(nèi)一側(cè)注入次流，與主流作用產(chǎn)生斜激波，從而使得主流方向改變獲得矢量推力，用該方法能獲得較大的矢量角。國內(nèi)外主要通過試驗和數(shù)值模擬兩種手段，以矢量角、矢量效率、推力系數(shù)為噴管性能指標，研究各參數(shù)變化對其影響及相關(guān)機理[3-10]。數(shù)值模擬方面，大多是對噴管進行定常模擬，較少研究動態(tài)變化參數(shù)對噴管性能的影響[11]。吳雄等[12]對固體火箭發(fā)動機的啟動和關(guān)閉的非定常特性做了研究，指出動態(tài)與穩(wěn)態(tài)特性相差很大，瞬態(tài)側(cè)向控制力、燃燒室壓力、推力等都會顯著變化，因此對于燃氣二次噴射系統(tǒng)的非定常特性研究，將成為飛行器控制系統(tǒng)設計的重點。Saha等[13]利用FLUENT軟件，對激波管驅(qū)動的超聲速噴管的啟動過程進行了時間相關(guān)的數(shù)值模擬，并與實驗值比較；計算捕捉到了流場的主要特征，表明利用商用軟件可實現(xiàn)對瞬態(tài)激波相互作用流場的計算。

本文以一二元收擴噴管為研究對象，利用FLU?ENT軟件，對次流質(zhì)量流量從零增加到最大值這一過程進行數(shù)值模擬，研究不同落壓比、不同次流加速時間下噴管特性的變化規(guī)律，詳細分析該過程中流場結(jié)構(gòu)變化，以期為射流推力矢量噴管研制提供參考。

2 計算方法及算例驗證

2.1 計算格式

通過時間推進的有限體積法求解二維非定常雷諾平均N-S方程，湍流模型為RNGk-ε模型。近壁采用標準壁面函數(shù)，空間離散采用二階迎風格式，時間離散采用二階隱式格式。

2.2 幾何模型與網(wǎng)格劃分

采用的噴管模型為NASA蘭利研究中心的試驗模型，設計落壓比NPRD=8.78，詳細參數(shù)見文獻[14]。計算域如圖1所示，長度和寬度分別為噴管喉道高度的23、28倍。計算網(wǎng)格在壁面、次流噴縫、噴管喉道等物理量變化梯度較大位置加密處理。

圖1 模型計算域Fig.1 Model computational domain

2.3 邊界條件

噴管的主流入口設定為壓力入口條件，給定總壓、總溫；計算域上、下、后邊界都設為壓力出口條件，給定總壓(為環(huán)境大氣壓)和總溫；噴管壁面及計算域前邊界設為絕熱、無滑移壁面條件。次流入口設為質(zhì)量入口條件，使用FLUENT中的UDF功能實現(xiàn)對次流的質(zhì)量流量控制。次流質(zhì)量流量隨時間的變化函數(shù)為：

2.4 算例驗證

為驗證數(shù)值模擬的準確性，對網(wǎng)格數(shù)和時間步長進行了敏感性分析，并與文獻[14]中的試驗數(shù)據(jù)做了對比。

在NPR=4.60、次流總壓與主流總壓之比SPR= 0.7，對網(wǎng)格數(shù)為10萬、20萬、40萬的網(wǎng)格進行定常數(shù)值模擬；對NPR=4.60，ts=1.00 s，時間步長為0.01 s、0.005 s、0.002 5 s進行動態(tài)數(shù)值模擬。從圖2和圖3中結(jié)果可得出，數(shù)值模擬結(jié)果與所選取的網(wǎng)格數(shù)和時間步長無關(guān)，并與試驗符合較好。下面模擬均選擇網(wǎng)格數(shù)為20萬的網(wǎng)格，時間步長0.005 s。

圖2 噴管上、下壁面靜壓數(shù)值結(jié)果與試驗結(jié)果對比Fig.2 Computational and experimental static pressure on nozzle wall

3 計算及結(jié)果分析

對NPR=4.60、8.78、10.00，ts=0.25 s、0.50 s、1.00 s共9種工況進行計算。

3.1 噴管參數(shù)計算

無次流情況下，理想推力計算公式為：

式中：pt為噴管入口氣體總壓，Ath為噴管喉道面積，pamb為環(huán)境壓力。

有次流情況下，理想推力為主流產(chǎn)生的推力Fpri與次流產(chǎn)生的推力Finj之和。Fpri可直接用式(2)求出，F(xiàn)inj由式(3)計算。

圖3 矢量推力隨時間的變化Fig.3 Vector thrust vs.different time-steps size

式中：pt,inj,ave為次流進口質(zhì)量加權(quán)平均總壓，且分別為次流入口邊界網(wǎng)格單元面上的質(zhì)量流量和總壓，k、R為空氣的物性參數(shù)。

圖4 實際推力變化曲線Fig.4 Real thrust

圖5 推力系數(shù)變化曲線Fig.5 Thrust ratio

3.2 噴管特性分析

3.2.1 推力及推力系數(shù)

由于m˙inj是流動時間t的函數(shù)，由式(1)可知，當無量綱時間t/ts相等時，不同ts下的也相等。從圖4中可看出，NPR越大推力越大，不同NPR下推力隨t/ts的增大而增大。圖4(a)中t/ts=0.85附近，推力有明顯振蕩；圖5(a)中，隨著t/ts的增大，推力系數(shù)先增后減，在振蕩后趨于平穩(wěn)。推力系數(shù)剛開始增加是因為NPR=4.60時，噴管處于過膨脹狀態(tài)，初始流量小，流動損失小，但次流可改善過膨脹狀態(tài)使推力系數(shù)增加，而后t/ts繼續(xù)變大，流動損失增大，推力系數(shù)下降。圖5(b)、圖5(c)中，噴管都不處于過膨脹狀態(tài)，次流的注入會引起推力損失，推力系數(shù)基本上隨t/ts的增大而減小。

3.2.2 矢量角及矢量效率

圖6 矢量角變化曲線Fig.6 Vector angle

圖7 矢量效率變化曲線Fig.7 Vector efficiency

3.3 噴管流動動態(tài)過程分析

3.3.1 推力系數(shù)及推力隨t/ts的變化

如圖4、圖5所示，NPR=4.60時，隨著t/ts的增大推力系數(shù)先增后減，推力則單調(diào)增加；但在t/ts= 0.80～1.00時，推力和推力系數(shù)均有明顯振蕩。下面分析振蕩原因。

圖8示出了圖4(a)中ts=1.00 s、t/ts=0.80～1.00這段內(nèi)推力各分量的變化曲線。從圖8(d)中容易看出，t/ts=0.82、0.85、0.86、0.90四個時刻相差很大，t/ts= 0.85是一振蕩峰值。圖9為這四個時刻下的馬赫數(shù)云圖，可見噴管出口附近剛好有一復雜波系，且隨著m˙inj的增加，該波系逐漸由噴管出口外移動到噴管內(nèi)。t/ts從0.80增大至0.85的過程中，該波系部分跨越噴管出口截面，動量項減小較快，但壓力項增加得更快(圖8(a)、圖8(b))，結(jié)果是軸向推力突升，在t/ts=0.85處達到峰值(圖8(d))。t/ts=0.85～0.86這一段，該激波系中正激波部分已完全進入噴管內(nèi)(圖9(c))，動量項減小很快(圖8(a))，而此時壓力項增加有限，整體表現(xiàn)為推力又迅速減小。t/ts=0.86～0.90這一段，該波系則進一步進入噴管內(nèi)。通過分析可得出：在NPR=4.60時，噴管出口附近有一復雜波系，該波系從噴管出口截面外側(cè)進入噴管出口截面內(nèi)側(cè)，使推力產(chǎn)生振蕩。

NPR=8.78和10.00(圖5(b)、圖5(c))時，推力系數(shù)基本隨t/ts的增大而減小，但在t/ts=4.0～5.5這一段出現(xiàn)了小波動，且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時的馬赫數(shù)云圖Fig.9 Computational mach contours，with NPR=4.60，ts=1.00 s

圖10為推力各分量及理論推力的變化曲線?？梢姡簍s=1.00 s和ts=0.25 s的理論推力變化曲線幾乎重合(圖10(f))，都隨t/ts的增大而線性增大，且在t/ts= 0.40～0.50這一段內(nèi)，兩種ts下的理論推力增量都超過了實際推力增量(圖10(e)、圖10(f))。t/ts=0.42～0.44這一段內(nèi)，ts=1.00 s的噴管實際推力各分量斜率出現(xiàn)較大變化，這就是出現(xiàn)波動的原因(圖5(b))。t/ts=0.44～0.48這一段內(nèi)，ts=0.25 s的噴管實際推力各分量斜率出現(xiàn)較大變化，推力系數(shù)下降變慢。圖11、圖12分別為ts=1.00 s和ts=0.25 s下三個時刻的流場圖。圖10中，t/ts=0.42～0.44這一段內(nèi)，ts=1.00 s時實際推力增速變快。圖11(a)中t/ts=0.42時刻，噴縫后主流再次附著壁面，形成的角渦區(qū)封閉，此時無回流；而圖11(b)中顯示主流與壁面分離，角渦區(qū)不再封閉，出現(xiàn)回流；圖11(c)中的回流區(qū)更大。結(jié)合圖10、圖11可得出，圖10中實際推力增長斜率從t/ts=0.42開始突然變大，是因為噴縫后主流與壁面分離，出現(xiàn)回流，噴管出口外面壓力高的氣體進入噴管內(nèi)，改善了噴管的過膨脹狀態(tài)，增加了推力。這就使得圖5(b)中ts=1.00 s時的推力系數(shù)出現(xiàn)小的波動。ts=0.25 s時的流動也經(jīng)歷了這一過程(圖12)，只是實際推力增長斜率變化不大，且回流出現(xiàn)在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增加的變化規(guī)律與推力的變化規(guī)律基本相同，矢量角隨t/ts先增加，在t/ts=1.00后趨于穩(wěn)定。矢量效率隨t/ts的變化較復雜，NPR=4.60時，如圖7(a)所示，矢量效率先急劇變小而后急劇增大，再相對緩慢變小至穩(wěn)定。NPR= 8.78時，如圖7(b)所示，矢量效率先急劇變小而后急劇增大，再緩慢減小至平穩(wěn)一段后變大，之后變小趨于穩(wěn)定。NPR=10.00與NPR=8.78的類似，變化曲線整體呈階梯狀。

圖7(b)、圖7(c)中，矢量效率在t/ts=0.40～0.55區(qū)間內(nèi)增大，隨次流的增加出現(xiàn)了回流，而回流(圖11、圖12)改善了次流后面噴管的過膨脹狀態(tài)，使得噴管上壁面壓力升高，有助于增大矢量角。

3.3.3 ts對推力系數(shù)及推力的影響

圖4中，NPR=4.60時三個ts值下的推力系數(shù)變化趨勢完全一致；NPR=8.78和NPR=10.00時，三個ts值下的推力系數(shù)變化趨勢基本一致，但在t/ts= 0.40～0.55，不同ts值下波動出現(xiàn)的時刻和強度不同。ts=1.00 s的波動出現(xiàn)得最早也最明顯，ts=0.25 s的波動出現(xiàn)得最晚(這里早、晚是在無量綱時間下而言)也最弱。波動出現(xiàn)時刻不同的原因，是因為ts= 0.25 s情況下次流增加較快，主流來不及響應，所以分離出現(xiàn)較晚；同樣可解釋ts=1.00 s時分離出現(xiàn)較早的原因。ts=0.25 s時推力系數(shù)波動最為平坦，是因為相比ts=1.00 s，從回流產(chǎn)生至其影響最大需要更長時間，這點可從圖10中看出，所以波動不明顯。

圖11 NPR=8.78、ts=1.00 s時的馬赫數(shù)云圖和流線圖Fig.11 Computational Mach contours and streamlines,with NPR=8.78,ts=1.00 s

圖12 NPR=8.78、ts=0.25 s時的馬赫數(shù)云圖和流線圖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下定常流動時的流場圖?？梢姡嗤瑃/ts下，ts=0.25 s沒有回流出現(xiàn)，ts=1.00 s回流區(qū)最大，這時的流動最接近該工況下的定常流動。這也是t/ts=0.44下ts=1.00 s的矢量角最大、ts=0.25 s的矢量角最小的原因。不同ts下，矢量效率隨t/ts的變化與矢量角的完全一致。

圖13 NPR=8.78、t/ts=0.44時噴管上壁瞬態(tài)壓力及該次流流量下的定常靜壓分布Fig.13 Pressure distribution along upper wall at transient and steady condition,withNPR=8.78,t/ts=0.44

4 結(jié)論

(1)NPR=4.60時，噴管處于過膨脹狀態(tài)，當次流流量增加到一定程度后，噴管推力隨次流的增加而振蕩。其原因是噴管出口附近存在一復雜波系，隨著次流的增加由噴管出口外向噴管內(nèi)移動，導致推力劇烈變化。

(2)NPR=8.78與NPR=10.00下，隨著次流的增加，推力系數(shù)總體呈下降趨勢，此過程中會出現(xiàn)一定波動，原因是次流流量達一定值后，噴縫后的主流與噴管壁面分離，出現(xiàn)回流，使得次流噴縫后噴管壁面的負壓狀態(tài)得到改善，減小了推力損失。

(3)回流的產(chǎn)生有利于提高推力系數(shù)、矢量角和矢量效率。

(4)在無量綱時間尺度下，次流加速時間越短，回流產(chǎn)生就相對越滯后。

圖14 馬赫數(shù)云圖與流線圖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]張群鋒，呂志詠，金捷，等.軸對稱射流矢量噴管的試驗和數(shù)值模擬[J].推進技術(shù)2004，25(2)：139—143.

[8]王慶偉，張相毅，徐學邈，等.射流角度對雙縫射流噴管流場影響的研究[J].燃氣渦輪試驗與研究，2006，19(4)：27—29.

[9]張強，楊永，李喜樂.雙喉道射流推力矢量噴管的數(shù)值模擬研究[J].西北工業(yè)大學學報，2009，27(6)：754—759.

[10]秦亞欣，于軍力，高歌.激波誘導圓形矢量噴管數(shù)值研究[J].航空動力學報，2009，24(10)：2208—2212.

[11]宋亞飛，高峰，馬岑睿，等.激波誘導矢量噴管內(nèi)流場動態(tài)特性數(shù)值研究[J].系統(tǒng)仿真學報，2013，25(1)：36—40.

[12]吳雄，焦紹球，張為華，等.燃氣二次噴射對推力矢量控制發(fā)動機內(nèi)流場的影響[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-)，男，江西九江人，碩士研究生，研究方向為發(fā)動機內(nèi)流氣動力學。