三擾流片推力矢量固體火箭發動機推力調節特性研究①

張汝衡,楊 軍,姚保江,楊石林,張兵峰

(1. 中國航天科技集團有限公司四院四十一所,西安 710025; 2. 中國航天科技集團有限公司第四研究院,西安 710025)

0 引言

隨著艦艇導彈防御技術的發展,要求小型艦載導彈具有發射快、定高懸停、抗風干擾等特點。常規舵機操縱舵面偏轉的控制方式逐漸難以滿足作戰需要,因此導彈自身的控制裝置顯得尤為重要[1]。擾流片式推力矢量裝置一般安裝在導彈底部,通過舵機控制擾流片偏轉產生側向力,對導彈進行俯仰、偏航兩通道的姿態調控。其產生的側向力具有響應快、抗干擾性強等優點,能夠大幅提高舷外導彈的機動性和敏捷性[2]。研究擾流片式推力矢量發動機的推力調節特性,對導彈低速大攻角飛行過程中的姿態穩定顯得尤為重要。

目前,國內外學者針對擾流片的研究主要集中在流場特性和主推力特性兩方面。KIM等[3]通過數值模擬研究了擾流片振蕩下的流場特性,論證了擾流片在主動控制技術中的應用潛力。GONZALEZ[4]通過對擾流片不同插入深度下的流場形態進行分析,解釋了發動機主推力曲線變化的原因。以上研究都為進一步分析擾流片表面氣動力特性提供了依據。HOLLSTEIN[5]從理論層面提出了噴管壁面和擾流片表面存在壓差,壓差導致了軸向力和側向力的產生。WENTZ等[6]進一步論證了不對稱的壓力分布源自噴管壁面的激波分離,壓差分別沿噴管軸向和垂直噴管軸向投影產生軸推力和側向控制力。YEUNG和EVERITT等[7-8]通過實驗研究了不同傾斜角度的擾流片的表面壓力,認為具有適當傾斜角度的擾流片能夠改善飛行器的氣動特性。崔業兵等[9]和韓文超等[10-11]通過數值模擬研究了擾流片伸入高度、阻塞面積率和推力偏轉角對氣動力大小的影響。叢戎飛等[12]研究了在軸對稱噴管中不同幾何外形的擾流片的氣動特性,王曉輝等[13]研究了擾流片的氣動設計原則,二人分別從優化擾流片型面和減小燃氣核心射流的角度提出了減小推力損失的方法。大量關于擾流片氣動特性的研究為進一步探究擾流片矢量調節特性奠定了基礎。目前主推力調節大多采用擾流片平動遮擋發動機噴管出口的方式[14],例如李國占[15]建立的擾流片式推力調節模型中,兩對擾流片通過“一進一出”的平移運動進行推力調節。但隨著舷外導彈小型化發展過程中噴口尺寸不斷縮小,擾流片平動的運動方式逐漸顯露局限性,需要著眼于新的作動方式,如擾流片旋轉。馬武舉等[16]研究了定軸旋轉的擾流片其旋轉軸和旋轉角對主推力的影響,但未研究擾流片非對稱旋轉下的側向控制力。

本文研究一種擾流片組合差動的推矢調節裝置,用于導彈快響應系統中的推力調節和姿態調節。首先采用重疊網格技術對舵機伸縮驅動擾流片旋轉的動態過程展開數值模擬,再將數值結果通過RBF神經網絡方法進行擬合,從而獲得依據舵機位移指令計算發動機主推力和側向力的推力矢量計算方法,并通過點火試驗驗證該計算方法的正確性。

1 數值模擬

1.1 幾何模型

本文研究的推力矢量裝置的結構簡圖如圖1所示。圖1中梯形擾流片通過帶轉軸的連接件與舵機相連,舵機位移變化范圍為0～20 mm。試驗測得舵機位移與擾流片旋轉角之間近似為線性關系,舵機每行進1 mm,擾流片旋轉1.5°。擾流片在裝載位置,舵機位移為0 mm,擾流片與噴管軸線夾角為53°;舵機位移增加到最大值時,擾流片與噴管軸線夾角為23°。擾流片從裝載位置到最大插入流場過程中,旋轉角變化30°,擾流片插入流場的程度可以根據擾流片與噴管軸線的夾角與旋轉角區間的比值計算。由于旋轉角與舵機位移近似呈線性關系,也可近似按照舵機位移的占比計算。噴管擴張段型面參數和擾流片型面參數如表1所示。

(b)YOZ diagram of schematic structure圖1 推力矢量裝置結構簡圖Fig.1 Schematic diagram of thrust vector device

表1 推力矢量機構的外形尺寸Table 1 Profile parameters of thrust vector device

該裝置的矢量調節原理:舵機伸縮傳遞位移信號,伺服機構響應位移指令控制擾流片定軸旋轉,使擾流片以不同深度插入噴管尾部流場。燃氣流在擾流片插入后流動方向發生突變,噴管擴張段內出現斜激波以及邊界層分離的現象,造成擾流片周圍出現壓差。當圖1中圓周分布的擾流片對稱旋轉時,側向氣動力相互抵消,推力矢量裝置只產生軸向力。該力與發動機主推力方向相反,稱為推力損失,起到對主推力的調節作用;當擾流片非對稱旋轉時,推力矢量裝置同時產生側向力和軸向推力損失,其中Z、X方向的側向力分量可對導彈進行俯仰-偏航方向的聯合姿態調控。

由推力矢量裝置的作動原理可知,舵機位移與擾流片旋轉角度之間可由幾何關系計算得到。故在建立仿真模型時,可對推力矢量裝置簡化,只保留噴管和擾流片,將轉軸簡化為一條軸線,便于后續對擾流片組合旋轉下的矢量特性進行仿真計算。根據表1列出的型面參數建立幾何模型如圖2所示。

圖2 幾何模型Fig.2 Geometry model

1.2 網格劃分

對1.1節建立的幾何模型進行結構網格劃分,采用重疊網格技術對擾流片旋轉進行數值模擬。重疊網格處理部件旋轉的網格運動問題簡便高效,不易出現負體積報錯,且能夠提高結構網格對不規則形狀的適應能力。圖3為背景網格,圖4為部件網格之一,重疊區域是直徑為67 mm的球面。網格總數共3 400 000,沖擊壁面采用邊界層網格加密處理,通過流體最大雷諾數計算y+值,保證y+<2。

1.3 數值方法

考慮發動機內流場粘性影響顯著,且擾流片所在位置存在激波和附面層分離,故湍流模型選用Realizablek-ε方程。假設流場中的流動為定常流動,忽略燃氣與噴管壁面的摩擦傳熱,不考慮流動過程中化學反應。設定燃燒室入口為壓力入口,總壓6.75 MPa,湍流強度5%,溫度3300 K。噴管出口為壓力出口,出口壓強為標準大氣壓,溫度300 K。噴管內壁及擾流片表面均設為絕熱無滑移邊界,流體工質選用理想氣體的燃氣。選擇密度基進行三維瞬態數值求解,流動方程選擇二階迎風格式,湍流方程選擇一階迎風格式,庫朗數取為1。壓力-速度耦合采用Simple算法,壓力、速度和湍流變量的收斂標準為相對殘差小于1×10-4。

圖3 背景網格Fig.3 Background grid

圖4 部件網格Fig.4 Component grid

通過CFD軟件進行擾流片表面壓強對面積的積分得到氣動阻力,發動機主推力F、側向力FR分別按式(1)、式(2)計算:

(1)

(2)

式中mt為噴管出口質量流率;ve為噴管出口燃氣流速;pe為噴管出口截面壓強;pa為環境大氣壓強;Ae為發動機噴口面積;Fy、Fz、Fx分別為軸向、俯仰方向、偏航方向的氣動阻力。

1.4 網格無關性驗證

為使數值計算結果與網格數量無關,以舵機位移為20 mm時擾流片裝置的推力損失為指標,進行網格無關性驗證。按照1.3節敘述的仿真設置,對網格數量分別為1 740 000、3 470 000、4 920 000和5 540 000的網格模型進行數值計算,不同網格數的推力損失如圖5所示。

圖5 網格無關性驗證Fig.5 Grid independence verification

結果表明,網格數量超過3 470 000時,擾流片軸向推力損失受網格數量影響較小。考慮到計算成本和精度,在后續數值計算中選用3 470 000的網格模型。

2 結果與討論

2.1 主推力調節特性

三個舵機均按照圖6所示的位移指令控制擾流片旋轉,研究擾流片對稱旋轉下的主推力調節特性。通過數值計算得到主推力和軸向推力損失如圖7所示。

在1～14 s內,舵機位移遞減,舵機回縮,控制擾流片旋轉角減小,擾流片插入噴管尾部流場的深度增加。由圖7可見,主推力和軸向推力損失都隨舵機位移減小而增加。

圖7 主推力調節特性Fig.7 Regulating characteristic of main thrust

在15.8～17 s內,舵機對主推力進行12.5 mm以內的小位移調節,主推力隨舵機位移的變化關系如圖8所示。由圖8可見,在15.8～16.2 s、16.7～17 s內舵機絕對位移較大,主推力曲線在舵機位移突變段斜率較小,推力臺階變化不明顯;在16.2～16.7 s內舵機絕對位移較小,主推力曲線在舵機位移突變段斜率較大,推力臺階變化較為明顯。

圖8 主推力變化曲線Fig.8 Main thrust change curves

進一步研究推力臺階變化處的噴管尾部流場形態。提取圖8中主推力曲線上的4個點,編號為1～4,分別對應舵機位移為12.5、10、5.5、0.5 mm。仿真得到擾流片附近的流場馬赫數云圖如圖9所示。由圖9可見,噴管出口的燃氣流受到擾流片擾動開始減速增壓,超聲速流動貼近擾流片表面產生斜激波,斜激波的出現造成了主推力曲線的波動。舵機位移為12.5 mm時,流場中僅存在斜激波,如圖9(a)所示;舵機位移為10 mm時,斜激波交匯處形成正激波,此時流場內正激波與斜激波并存,如圖9(b)所示;舵機位移為5.5 mm時,斜激波幾乎消失,此時噴管尾部流場內形成馬赫盤,如圖9(c)所示,馬赫盤的出現造成流場劇烈變化,主推力曲線波動幅度加劇,臺階變化明顯;舵機位移為0時,插入深度達到最大,馬赫盤逐漸向噴管出口處前移,貼附于噴管口,如圖9(d)所示。

圖10為1～4采樣點處擾流片表面的壓力梯度云圖。可見,在擾流片插入流場過程中,壓力梯度較大的區域分布在沿擾流片底部邊緣。

2.2 側向力調節特性

擾流片式推力矢量發動機產生的側向力主要包括擾流片非對稱旋轉受到的側向力和噴管內流動受擾流片影響產生不對稱而產生的噴管內側向力。

對于擾流片受到的側向力,圖11為從噴管尾部看去舵機安裝示意圖,舵機1控制擾流片旋轉產生俯仰方向的負力,舵機2、3控制擾流片產生俯仰方向的正力;舵機1、2組合控制擾流片產生偏航方向的正力,舵機1、3組合控制擾流片產生偏航方向的負力。

F1、F2、F3為三個擾流片表面氣動力沿圓周法向的投影,根據圖11所示幾何關系,俯仰、偏航方向的側向力計算公式為

(3)

(4)

三個舵機按照圖12所示的位移指令控制擾流片旋轉,研究擾流片組合差動下的側向力調節特性。側向力隨舵機位移的變化關系如圖13所示。

(a)Mach diagram at point 1 (b)Mach diagram at point 2

(c)Mach diagram at point 3 (d)Mach diagram at point 4圖9 主推力曲線上特征拐點處的馬赫線圖Fig.9 Mach diagram at characteristic inflection point on main thrust curve

(a)Pressure gradient (b)Pressure gradient (c)Pressure gradient (d)Pressure gradient at point 1 at point 2 at point 3 at point 4圖10 主推力曲線上特征拐點處擾流片表面壓力梯度云圖Fig.10 Pressure gradient contours of jet tabsurface at characteristic inflection point on main thrust curve

圖11 舵機及擾流片安裝示意圖Fig.11 Installation diagram of steering gear and jet tabs

圖12 舵機位移指令Fig.12 Displacement command of steering gear

圖13 擾流片上側向力隨舵機位移變化曲線Fig.13 Lateral force on tabs varies withdisplacement of steering gear

由圖13可見,在0～3 s內,舵機2、3保持位移為8 mm,僅舵機1位移從15 mm減小到8 mm。對應的擾流片受到的俯仰方向側力Fz隨L1減小而增加,偏航方向側力Fx隨L1增加保持不變;在6～9 s內,舵機1位移保持不變,舵機2、3位移同時減小,造成Fz減小,Fx保持不變。表明舵機1控制擾流片產生俯仰方向的負力,舵機2、3控制擾流片產生俯仰方向的正力。在3～6 s內,舵機2位移減小,造成Fz和Fx均隨L2減小而減小;在9～12 s內,舵機1、2位移同步減小,Fz隨L1、L2減小而增加,Fx隨L1、L2減小而減小。表明舵機2控制擾流片產生俯仰和偏航方向的正力,舵機1、2控制擾流片產生俯仰方向的負力和偏航方向的正力。

由此可見,數值計算得到的側向力變化趨勢與式(3)、式(4)描述一致。對于流動不對稱產生的噴管內側向力,數值計算得到的噴管內側向力如圖14所示。

圖14 噴管內側向力隨舵機位移變化曲線Fig.14 Curves of force inside nozzle withdisplacement of steering gear

由圖14可見,在0～3 s和6～9 s內,分別由舵機1位移減小和舵機2、3位移同時減小控制擾流片旋轉。對應的噴管內流動不對稱性產生的噴管內側向力,俯仰方向的噴管內側向力Fz_p隨舵機1位移減小而增加,隨舵機2、3位移同時減小而減小;偏航方向的噴管內側向力Fx_p基本保持不變。表明擾流片1旋轉造成噴管受到俯仰方向的負力,擾流片2、3同步旋轉造成噴管受到俯仰方向的正力。在3～6 s和9～12 s內,分別由舵機2位移減小和舵機1、2位移同時減小控制擾流片旋轉。相應地,Fz_p隨舵機2位移減小而減小,隨著舵機1、2位移同時減小而增加;Fx_p隨舵機2位移減小和舵機1、2位移同時減小而增加。表明擾流片2旋轉造成噴管受到俯仰和偏航方向的正力,擾流片1、2同步旋轉造成噴管受到俯仰方向的負力和偏航方向的正力。

根據式(5)計算側向力調節幅度:

(5)

式中FR為側向控制力;Forg為未調節推力。

計算得發動機側向力調節幅度為42.3%,表明擾流片式推矢裝置的能夠起到良好的側向力調節作用。

綜上所述,對于擾流片非同步旋轉產生的側向力和噴管內流動不對稱造成噴管受到的側向力,數值計算結果表明,俯仰方向的側向力由舵機1、2、3的位移組合決定;偏航方向的側向力由舵機2、3的位移組合決定。由于實際控制中俯仰-偏航通道存在耦合,舵機1位移也對偏航方向的側力有所影響。

3 按位移-調節幅度擬合的推力矢量計算方法

為進一步獲得推力矢量與三個舵機位移之間的量化關系,本節采用RBF神經網絡(Radial Basis Function Neural Network)方法對數值樣本進行擬合,獲得矢量調節模型。

對圖12所示的0～12 s的位移指令,從對應的推力矢量數值結果中,每隔3 ms選取一組結果作為樣本,共得到4000組位移組合(L1,L2,L3)下的軸向推力損失Fy以及俯仰、偏航方向的側向力Fz、Fx。設δFx、δFy、δFz為推力調節量,則推力矢量Fx、Fy、Fz與未調節推力Forg之間的關系表示為

Fx=δFx·Forg,Fy=δFy·Forg,Fz=δFz·Forg

(6)

處理樣本數據可以得到對應的δFx、δFy、δFz值。將舵機位移組合(L1,L2,L3)作為輸入數據集,(δFx、δFy、δFz)作為目標數據集,利用神經網絡工具箱對數據進行訓練及擬合,獲取位移-調節幅度函數關系。

徑向基(RBF)神經網絡對于高維問題的近似具有強大的非線性函數擬合能力,其計算精度高、計算速度快[17-22],可作為推力矢量調節函數關系式的擬合方法。神經網絡包含輸入層、輸出層和隱含層三部分,其原理是利用NT個樣本點按式(7)進行建模:

(7)

(8)

采用默認設置建立兩層神經網絡模型,包含輸入層、隱含層和輸出層。設置隱含層神經元數目為15,將數據的70%作為訓練集,15%作為測試集,15%作為驗證集,算法選擇Levenberg-Marquardt得到矢量調節模型。訓練結果及相關系數如圖15、圖16所示。

(a)Training of δFx

(b)Training of δFy

(c)Training of δFz圖15 神經網絡訓練結果Fig.15 Training results of neural network

(a)Correlation coefficient of δFx

(b)Correlation coefficient of δFy

(c)Correlation coefficient of δFz圖16 擬合結果Fig.16 Fitting results

由圖16可見,經過訓練后,模型誤差達到了10-6數量級,相關系數R達到0.99以上,擬合良好。說明擬合得到的推力矢量調節函數可以用于計算主推力和側向力。

擬合得到的推力調節幅度計算公式:

δFi=a0+a1·L1+a2·L2+a3·L3+a4·L1L2+

(9)

式中δFi(i=z,x,y)分別為推力矢量發動機在俯仰、偏航和軸向的推力調節幅度;a0～a9為位移前系數,如表2所示。

表2 推力調節幅度擬合系數Table 2 Fitting coefficients of thrust amplitude

根據擬合結果,可以計算得到擾流片式推力矢量發動機的主推力和側向力:

F=(1-δFy)·Forg

(10)

(11)

4 試驗驗證

為驗證第3章擬合得到的推力矢量計算式的正確性,進行原理樣機推力開環試驗。試驗選用高精度六分力試驗臺測量推力矢量。六分力試車原理圖如圖17所示。

六分力試驗臺上布置有6個傳感器,使發動機處于靜定平衡狀態。根據剛體平衡原理測量發動機的主推力和側向力,圖17中,規定傳感器受拉為正,受壓為負,力矩正方向按右手螺旋法則確定。選擇直角坐標系Oxyz,1～6號傳感器分別測量6路分力F1～F6,由空間力系平衡條件可得

FX=-(F1+F2+F4),FZ=-(F3+F5),

FY=-(F6+mg)

(12)

試驗發動機噴管喉徑為10 mm,測得燃燒室壓強為6.75 MPa,溫度為3300 K;裝藥量為6.86 kg,比沖為2030 N·s/kg。根據試車壓強積分、裝藥量可獲得特征速度為1382 m/s,根據藥柱肉厚和燃燒時間得到燃速為6.4 mm/s。推進劑密度為1.62 g/cm3,發動機采用石墨喉襯,試車過程中無燒蝕。

按照圖18所示的舵機位移指令進行試驗。將位移指令代入第3章獲得的推力矢量調節計算式,計算得到的推力矢量與試驗結果對比如圖19所示。

由圖19可見,計算式計算結果與試驗結果基本相同。試驗表明,對于相同的舵機位移指令,采用擬合關系式計算與點火試驗得到的主推力偏差在4%以內,側向力偏差在3%以內。誤差在可接受范圍內,認為基于RBF神經網絡進行擬合的推力矢量計算式具有有效性。

造成偏差的原因,主要包含數值計算結果的偏差、試驗工況及試驗臺測量的偏差、神經網絡訓練過程的偏差三方面。對于數值結果的偏差,考慮到數值模擬不能完全貼近試驗工況,如采用理想燃氣、不考慮兩相流、忽略燃氣與壁面之間的摩擦傳熱等假設;對于試驗工況和測量的偏差,考慮到發動機地面點火試驗的復雜性,擾流片可能存在燒蝕、以及試驗中存在發動機推力偏心等;對于神經網絡訓練的偏差,考慮到樣本數據集有限、隱含層和隱含單元數目等方面。

(a)Comparison of main thrust

(b)Comparison of lateral force圖19 擬合計算結果與試驗實測值比較Fig.19 Fitting calculation compared with measured value

5 結論

(1)擾流片對稱旋轉時,推力損失與流場插入程度呈正相關。舵機進行12.5 mm以內的小位移調節時,發動機主推力曲線臺階變化明顯。三擾流片的舵機位移組合為(10.5 mm,10.5 mm,10.5 mm)時,流場內出現斜激波交會、馬赫盤前移的現象。

(2)擾流片組合差動時,推矢機構產生的側向控制力為舵機位移的函數。

(3)通過RBF神經網絡方法對數值計算結果進行訓練及擬合,其模型誤差達到了10-6數量級,相關系數R達到0.99以上,擬合度良好。由此提出了一種按舵機位移和推力調節幅度計算發動機主推力和側向力的推力矢量計算方法。

(4)試驗表明,在燃燒室總壓6.75 MPa的工況下,基于RBF的推力矢量計算方法得到的主推力預示結果偏差在4%,側向力預示結果偏差在3%以內,驗證了推力矢量計算方法的有效性。