喉栓式可調噴管推力模型優化研究

陳曉麗，牛 祿，潘科瑋，王志新，陳浩然

(上海航天動力技術研究所，上海 201109)

0 引言

固體姿軌控發動機利用燃氣直接橫向噴射,獲得作用于飛行器的橫向推力,從而產生機動控制力或力矩,可應用于導彈[1]、逃逸系統[2]等。喉栓式可調噴管是固體姿軌控發動機實現推力調節的重要部件,通過調節喉栓位置,可控制工質流量,實現推力的調節。

20世紀60年代,美國開始研究喉栓式推力可調固體發動機,并將喉栓式可調噴管應用于標準-3導彈的推力連續調節固體姿軌控動力系統[3], 2013年完成了鑒定試驗,2014年標準-3 Block IB導彈進入海軍正式部署階段。國內在喉栓式可調噴管的理論分析、數值模擬和試驗研究等方面也開展了較多研究[4-6],但工程應用尚有差距。噴管推力是喉栓式可調噴管的重要性能,也是固體姿軌控發動機推力調控系統設計的基礎。喉栓式可調噴管的推力常采用數值模擬的方法獲得[7-9],該方法在獲得噴管推力等性能的同時,還能對噴管內流場進行分析,但固體姿軌控動力系統存在諸多非穩態,噴管將面臨多種工況,增加了數值模擬的工作量及工作難度,數值模擬的大量處理及計算工作也不利于噴管的設計優化。此外,現有喉栓式可調噴管的推力測量手段較少,測量難度較大,也給喉栓式可調噴管的推力性能研究帶來困難。

針對喉栓式可調噴管,建立了幾何喉部面積計算模型和理論推力模型,研究幾何喉部面積和噴管理論推力隨喉栓位置的變化情況,并建立了三維仿真模型,對喉栓在不同位置下的噴管內流場進行了數值模擬,分析喉栓在不同位置時的噴管內流場分布,研究噴管質量流率、推力隨喉栓位置的變化規律,與理論推力進行了對比。建立了喉栓式可調噴管氣動等效喉部面積的計算方法,對等效喉部面積計算模型和推力模型進行了優化,并進行了冷氣試驗驗證。研究結果可為后續喉栓式可調噴管的設計及優化、固體姿軌控發動機推力調控系統的設計提供參考。

1 喉栓式可調噴管理論推力模型

1.1 喉栓式可調噴管工作原理

喉栓式可調噴管的結構如圖1所示,進氣段與噴管軸線成一定角度,喉栓與噴管內壁形成環形氣體通道,即等效喉部,喉栓式可調噴管通過改變喉栓的位置改變等效喉部面積,從而控制噴管輸出推力。喉栓處于初始位置時,噴管開度最大,此時噴管的等效喉部面積最大,噴管的推力達到最大值,喉栓向外移動,噴管開度減小,通過的工質流量減小,噴管推力減小,當喉栓達到最大位移時,噴管完全關閉,噴管推力下降為0。

圖1 非同軸喉栓式可調噴管示意圖Fig.1 Non-coaxial pintle flow regulation nozzle schematic diagram

在固體姿軌控發動機推力調控系統的設計及仿真研究中[10-13],等效喉部面積計算模型和推力模型是整個系統的重點,等效喉部面積取喉栓和噴管型面構成的幾何喉部面積,即噴管流道面積最小處的面積,推力取基于等熵流動假設的理論推力。

1.2 幾何喉部面積計算模型

當噴管收斂段和喉栓頭部型面均為圓弧時,設噴管軸線為x軸,方向如圖2所示,O1,O2分別為噴管收斂段和喉栓頭部型面圓弧中心,其半徑分別為R1,R2,圓心初始坐標為(xo1,yo1),(xo2,yo2),設喉栓向外運動了x(0 ≤x≤xmax,xmax為喉栓伸入喉部的最大位移,此時喉栓型面與噴管型面相切,噴管處于關閉狀態),則O2坐標變為(xo2+x,yo2),A,B分別圓O1和圓O2上的一點,坐標為(xA,yA),(xB,yB),RA,RB分別為A,B到噴管軸線的距離,則其環道面積S為[10]:

(1)

S的最小值即噴管的幾何喉部面積At,At的求解是一個有界的非線性規劃問題,可以采用內點法進行求解。At是x的函數,且隨著x近似線性變化[9]。

圖2 幾何喉部示意圖Fig.2 Geometric throat schematic diagram

1.3 理論推力模型

噴管的推力是噴管內、外表面所受氣體壓力的合力,根據動量守恒定律,可得推力公式[14]:

(2)

基于一維等熵流動假設,噴管質量流率和噴氣速度分別為[14]:

(3)

(4)

(5)

式中:R,k,Tf,pc分別為氣體常數、比熱比、推進劑絕熱燃燒溫度和噴管入口截面處的壓強。代入推力公式得噴管的理論推力F為[14]:

F=CFpcAt

(6)

(7)

入口壓強pc、喉部面積At、出口面積Ae已知,求得噴管出口壓強即可得到推力大小。膨脹壓強比與噴管面積比關系為[14]:

(8)

根據壓強比和面積比的關系,可由噴管面積比迭代求得對應的膨脹壓強比,進而得到喉栓在不同位置時的噴管推力輸出。

設噴管入口壓強為9 MPa,出口壓強為大氣壓強,工質比熱比為1.25,喉栓從位置0運動到xmax處,噴管幾何喉部面積和理論推力等性能的變化如圖3所示。其中,F0和At0分別為喉栓處于初始位置時噴管的理論推力和幾何喉部面積。噴管理論推力和幾何喉部面積在喉栓處于初始位置時達到最大值,兩者隨喉栓向外運動而近似線性減小,幾何喉部面積在喉栓達到最大位移處減小為0,而理論推力在喉栓達到最大位移前就已減小為0,此后出現負值,這是由于上述推力模型僅適用于噴管的完全膨脹、欠膨脹和輕微過膨脹狀態,在喉栓運動后期,噴管擴張比劇增,噴管處于過膨脹狀態,公式已不再適用。

圖3 噴管性能隨喉栓位置的變化Fig.3 The variation of nozzle performance with the position of pintle

2 喉栓式可調噴管內流場數值模擬

為了研究喉栓位置對噴管性能的影響,建立更為準確的推力模型,對喉栓在不同位置下的噴管內流場進行了數值模擬,分析喉栓在不同位置時的噴管內流場分布,研究噴管質量流率、推力等隨喉栓位置的變化規律。

2.1 數值計算方法

采用FLUENT進行內流場數值模擬,RANS方程表達式為[15]:

1)連續方程

(9)

2)動量方程

(10)

3)能量方程

(11)

采用RNGk-ε模型,可以更好地處理高應變率及流線彎曲程度較大的流動,其湍動能k和耗散率ε方程分別為[15]:

(12)

(13)

式中:Gk是由平均速度梯度引起的湍動能;Gb是由浮力引起的湍動能;YM代表可壓縮湍流中脈動膨脹的貢獻;C1ε,C2ε,C3ε為經驗常數,αk,αε分別為湍動能和耗散率對應的普朗特數;Sk和Sε為源項。

2.2 網格劃分及無關性驗證

幾何模型由噴管內型面、喉栓頭部型面等組成,根據結構的對稱性,計算采用二分之一三維模型。對計算區域分塊進行網格劃分,在進氣段等部位建立結構網格,對幾何形狀較為復雜的部位采用非結構網格,計算區域通過交界面進行參數傳遞。流體流經喉栓頭部和噴管收斂段時流速會發生較大變化,因此需要對網格進行局部加密,局部網格示意圖見圖4。

圖4 噴管內流場局部網格示意圖Fig.4 Local grid diagram of flow field in nozzle

將邊界條件設置為:

1)入口邊界條件,采用壓力入口邊界條件,設置入口壓強為9 000 000 Pa,總溫為1 800 K;

2)出口邊界條件,采用壓力出口邊界條件,設置出口壓強為101 325 Pa,總溫為300 K;

3)固體壁面條件,采用標準壁面函數,絕熱無滑移條件。

4)假設流體介質為單一理想氣體,不考慮多相流動,不考慮重力影響。

為排除各個工況因網格數量不同而造成的差異,需對已建立的網格進行無關性驗證。將網格數分別為758 046和1 846 890的網格模型進行比較,驗證噴管質量流率和推力大小,結果見表1。

表1 兩種網格模型的計算結果Table 1 Calculation results of two grid models

兩種網格模型計算得出的質量流率差異小于0.1%,推力值差異小于0.01%,可認為模型1已達到網格無關。

2.3 計算結果及討論

通過計算得到了噴管內流場、噴管推力等隨著喉栓伸入喉部的變化,圖5給出了喉栓向外運動時噴管對稱面上的壓強和馬赫數分布變化。

喉栓距離噴管喉部較遠時對流場產生的擾動很小,如圖5(a)所示,在噴管喉部及擴張段前半部分,壓強逐漸減小,由于存在輕微的過膨脹,在擴張段接近出口處壓強基本不再變化,馬赫數沿著噴管軸線方向均勻增加,在噴管喉部達到Ma1,此時噴管每秒質量流率和推力最大。

喉栓向外運動到一定位置后,燃氣在喉栓表面與噴管收斂段內壁之間加速膨脹,在喉栓頭部附近形成等效喉部,由于非同軸進氣,等效喉部先在進氣一側形成,如圖5(b)所示,此時等效喉部面積與幾何喉部面積相差較大。

喉栓繼續伸入喉部,在喉栓與噴管內壁圍成的環形通道處形成等效喉部,在喉栓頭部形成了渦流區,且喉栓頭部附近流場分布非對稱性明顯,如圖5(c)所示。隨著噴管開度減小,形成等效喉部的位置向上游移動且等效喉部面積減小,氣流在等效喉部后的膨脹加劇,喉栓頭部的渦流區縮小。喉栓處于全關位置附近時,噴管內流場分布愈加紊亂,噴管處于過膨脹狀態,出現氣流分離現象,如圖5(d)所示。

圖5 喉栓在不同位置時噴管對稱面上的壓強和馬赫數分布圖Fig.5 The pressure and Mach number distribution on the symmetry plane of the nozzle with different positions

噴管質量流率和推力隨喉栓位置的變化如圖6、圖7所示,隨著喉栓向外運動,噴管每秒質量流率和幾何喉部面積近似線性減小。噴管推力在喉栓運動初期近似線性減小,而隨著喉栓位移增大,對流場的擾動增加,且喉部面積減小,噴管擴張比增大,噴管過膨脹,推力隨喉栓位置的變化率呈減小的趨勢。

圖6 每秒質量流率隨喉栓位置的變化Fig.6 The change of mass flow rate per second with the position of pintle

理論推力隨喉栓位置的變化率大于CFD仿真結果,在喉栓運動初期,理論推力大于CFD結果,且隨著喉栓向外運動,兩者差距逐漸減小,當喉栓運動至某一位置時,理論推力與CFD計算結果相吻合,而喉栓繼續伸入喉部時,理論推力小于CFD計算結果且兩者差距隨喉栓伸入喉部而增大。

圖7 理論推力和CFD仿真結果對比Fig.7 Comparison of theoretical thrust and CFD simulation result

經分析,等效喉部在接近喉栓頭部的位置形成,燃氣經過喉栓頭部后流動方向發生較大變化,氣體的慣性作用較為明顯,存在流動損失,且由于非同軸的影響,喉栓上下型面形成等效喉部的位置不對稱,使得幾何喉部面積大于實際等效喉部面積,理論推力計算結果大于CFD仿真結果。隨著喉栓伸入喉部,喉部出現狹縫式喉徑,導致噴管的擴張比劇增,噴管處于過膨脹狀態,氣流分離加劇,不符合等熵流動假設,1.3節中的推力公式已不適用。因此,需要對推力模型進行優化。

3 喉栓式可調噴管推力模型優化研究

3.1 氣動等效喉部面積

等效喉部面積計算模型直接影響推力模型,是推力模型優化的基礎。定義噴管的氣動等效喉部面積計算公式為:

(14)

代入CFD仿真結果,可得噴管的氣動等效喉部面積。改變入口邊界條件和工質參數,出口邊界條件不變,得到噴管在不同位置時的等效擴張比Ae/At,q如表2所示。

表2 不同喉栓位置和工況下的等效擴張比Table 2 The equivalent expansion ratio under different pintle positions and working conditions

表3 各工況條件Table 3 Working conditions

由結果可知,喉栓在初始位置時,當入口壓強由9 MPa降為3 MPa,氣動喉部面積略減小,噴管等效擴張比增加了0.046 2(0.36%,工質為燃氣)和0.067 2(0.53%,工質為冷氣);當入口壓強不變,工質由燃氣變為冷氣時,氣動喉部面積略有增加,噴管等效擴張比減小了0.093 8(0.74%,入口壓強為9 MPa)和0.072 8(0.57%,入口壓強3 MPa),各工況下的噴管等效擴張比標準差為0.056。

喉栓向外運動1個單位后,入口壓強和工質發生變化時,對應的噴管等效擴張比變化率分別為0.28%,0.31%,0.47%,0.44%,各工況下的噴管等效擴張比標準差為0.058,較喉栓在初始位置時各工況下等效擴張比的標準差略有增加。

3.2 等效喉部面積計算模型

由于噴管內存在流動損失,噴管幾何喉部面積大于實際等效喉部面積,給噴管推力計算帶來誤差,為優化噴管推力計算模型,采用噴管氣動等效喉部代替推力模型中的幾何喉部。

由CFD仿真結果可知,喉栓在同一位置而工況不同時,噴管等效擴張比的差異很小,因此,忽略工況對氣動等效喉部面積的影響,近似認為氣動等效喉部面積隨喉栓位置的變化規律只與喉栓和噴管型面相關。

由圖3可知,幾何喉部面積與喉栓位置x近似成一次方函數關系,而根據CFD仿真結果,氣動等效喉部面積與喉栓位置x近似成三次方函數關系,因此建立以下修正公式:

(15)

應用非線性最小二乘法,求得系數a1=0.92,a2=-0.75/(yo1-yo2),a3=-1.25(R1+R2)/(yo1-yo2)。

幾何喉部面積、修正后的等效喉部面積及工況1下氣動等效喉部面積對比如圖8所示,氣動等效喉部面積小于幾何喉部面積,結合圖5,隨著喉栓向外運動,等效喉部形成的位置向上游移動,與流速方向接近垂直,與幾何喉部面積差距減小。修正后的等效喉部面積計算公式能較好擬合氣動喉部面積隨喉栓位置的變化規律。

圖8 等效喉部面積計算結果對比Fig.8 Comparison of equivalent throat area calculation results

3.3 推力模型優化

將修正后的等效喉部面積代入理論推力計算公式,計算得到工況1下噴管推力隨喉栓位置變化如圖9所示。

圖9 初步修正后的推力計算結果Fig.9 The thrust calculation result after preliminary correction

由圖可知,在喉栓運動初期,推力得到較好修正,與CFD仿真結果相吻合,當喉栓運動至某一位置后,推力計算結果與CFD仿真結果的偏差逐漸增大,此時噴管內流場不符合等熵流動假設,需要對推力計算模型再次進行修正,最終推力計算模型為:

(16)

(17)

(18)

(19)

式中:b1=46;b2=0.015;b3=-22.5/k;b4=-0.09。

上述推力模型對動推力和靜推力兩部分皆進行了修正,利用優化后的推力計算模型,計算得到工況1下噴管推力隨喉栓位置變化如圖10所示。

圖10 優化后的推力模型計算結果Fig.10 The optimized thrust model calculation results

表4 不同喉栓位置下的推力計算結果對比Table 4 Comparison of thrust calculation results under different pintle positions

將理論推力、優化模型計算結果與CFD仿真結果進行對比,結果如表4所示,均方根誤差由0.063減小至0.002,優化后推力模型計算結果與CFD仿真結果吻合較好。

4 試驗驗證及分析

4.1 試驗系統

試驗系統由對稱螺旋進氣管、壓強傳感器、電磁閥、靜態標定組件、噴管、板簧、螺桿和力傳感器等組成,如圖11所示。進氣管路在結構上采用了螺旋方式進氣且在兩側對稱分布,高壓氣瓶由減壓閥、高壓軟管與螺旋管路連接,通過減壓閥可初步調整氣體流出的壓強,進氣端通過連接板支架固定在底板上,由于氣流流經對稱螺旋管路后同時流出,能對氣流起緩沖和穩定的作用,最后經由三通管將穩定的氣流送入主管道。

圖11 試驗系統示意圖Fig.11 Diagram of test system

高壓氣體通過兩端對稱螺旋管路匯入主管路,經過壓強傳感器,待氣流穩定后由壓強傳感器記錄腔體壓強值,再由電磁閥控制主管路氣流的開閉,保障主管路氣流的穩定性,同時可以檢驗管路的密閉性是否良好。喉栓作動時,噴管朝上噴氣,通過板簧和螺桿將力傳遞給力傳感器。對整個試驗系統進行標定后,可進行噴管冷流試驗,測得不同腔體壓強下的噴管推力。

4.2 試驗結果及分析

測量喉栓在不同位置時噴管的穩定輸出推力,測試結果見表5。

腔體壓強分別在9 MPa和3 MPa左右時,推力模型計算結果和試驗結果的最大絕對百分比誤差分別為10.72%和 11.07%,均發生在x/xmax=0.9處,即噴管即將關閉時。經分析,推力模型計算結果與試驗結果存在偏差是由于推力模型的建立過程中存在簡化,比如忽略了工況對氣動喉部面積的影響、噴管外流場的影響等,且隨著喉栓伸入喉部,對流場的擾動增加,噴管內部流動復雜,給推力計算帶來了困難。此外,試驗系統存在進氣擾動、管路和板簧等帶來的推力傳遞損失等,也給推力測量帶來了誤差。

腔體壓強分別在9 MPa和3 MPa左右時,不同喉栓位置的噴管推力計算結果與試驗結果的平均絕對百分比誤差分別為3.10%和4.99%,修正后的推力模型計算結果與試驗結果基本吻合。

表5 試驗結果與推力模型計算結果Table 5 Test results and thrust model calculation results

4 結論

1)建立了三維仿真模型,可準確地模擬喉栓對噴管內流場的影響,隨著喉栓深入喉部,噴管推力和出口質量流率減小,喉栓對流場的擾動加劇,噴管擴張段出現氣流分離現象。

2)建立了氣動等效喉部面積的計算方法,修正了喉部面積計算模型,氣動等效喉部面積與入口邊界條件、工質參數、喉栓和噴管型面相關,入口邊界條件與工質參數對等效喉部面積的影響在一定范圍內可以忽略。

3)優化了推力模型,該模型可預示噴管推力隨喉栓位置的變化規律,喉栓向外運動時,噴管推力先近似線性減小,后喉栓伸入喉部,噴管推力不滿足線性變化,隨著喉栓向外運動逐漸降為0。

4)進行了噴管冷氣推力測試,獲得了不同喉栓位置和腔體壓強下的噴管推力,試驗結果與推力模型計算結果基本吻合。