張德權,賈軍凱,武澤平,王東輝,牛 祿,周偉華

(1.河北工業(yè)大學 機械工程學院,天津 300401;2.國防科技大學 空天科學學院,長沙 410073;3.上海航天動力技術研究所,上海 201109)

0 引言

固體姿軌控發(fā)動機(Solid Divert and Attitude Control Motor,SDACM)可以為飛行器軌道和姿態(tài)調(diào)整提供動力[1]。其工作原理是通過改變喉栓位移以控制閥門燃氣流量,從而調(diào)節(jié)單閥推力大小,進一步通過多個單閥推力合成不同大小和方向的合推力矢量,實現(xiàn)飛行器軌道和姿態(tài)控制[2]。固體姿軌控發(fā)動機中各種不確定性因素相互耦合,引起實際推力與名義推力之間的偏差,導致飛行器軌道和姿態(tài)難以達到預期控制效果[3]。為此,對固體姿軌控發(fā)動機推力偏差進行不確定性分析,研究推力調(diào)控過程中不確定性因素傳遞規(guī)律,為固體姿軌控發(fā)動機實現(xiàn)快速準確調(diào)控提供參考。

目前,研究發(fā)動機推力性能主要手段是計算流體力學(Computational Fluid Dynamics,CFD)和有限元分析等數(shù)值模擬方法[4],但是這些方法求解成本高,計算效率低。為此,代理模型方法被引入到固體發(fā)動機性能研究中以提高計算效率。代理模型僅需少量的仿真樣本就可以將原有復雜耗時仿真模型替換為簡單數(shù)學模型,顯著降低計算成本[5]。常用的代理模型有徑向基函數(shù)[6]、響應面模型[7]和Kriging模型等[8]。多位學者將代理模型方法用于構建固體發(fā)動機近似模型以降低計算成本。文謙等[9]使用徑向基函數(shù)模型對運載火箭進行氣動力和結(jié)構承載能力近似建模,大幅提高計算效率。TOLA等[10]通過響應面模型研究幾何參數(shù)與火箭發(fā)動機結(jié)構完整性和內(nèi)彈道性能之間的關系。王鵬宇等[11]采用基于Kriging模型的優(yōu)化框架,實現(xiàn)多工況下喉栓式固體發(fā)動機噴管型面優(yōu)化設計。

上述研究均從確定性角度出發(fā)分析固體發(fā)動機性能,未考慮固體發(fā)動機中的不確定性因素[12],為此,多位學者[13-16]考慮環(huán)境載荷、材料屬性和結(jié)構尺寸等多源不確定性因素對發(fā)動機藥柱貯存、藥柱結(jié)構和密封性能的不確定性進行分析。然而,固體姿軌控發(fā)動機推力性能不確定性相關研究還較為少見。固體姿軌控發(fā)動機推力性能同樣受到多源不確定性因素影響,主要包括伺服機構不確定性、加工精度不確定性以及推進劑材料屬性不確定性等。這些不確定性因素導致發(fā)動機推力實際值與名義值存在偏差,如果推力偏差超出控制系統(tǒng)容許偏差,會導致固體姿軌控發(fā)動機難以達到預期軌道控制效果。因此,考慮固體姿軌控發(fā)動機中推力影響因素的不確定性,對推力偏差進行不確定性分析具有重要意義。

本文采用Fluent進行固體軌控發(fā)動機單閥推力仿真計算獲取名義推力。然后,考慮參數(shù)不確定性進行仿真計算,得到軌控發(fā)動機單閥仿真推力樣本。根據(jù)仿真推力樣本和名義推力建立軌控發(fā)動機單閥推力偏差Kriging模型,實現(xiàn)推力快速預測。在推力偏差Kriging模型基礎上,對其進行不確定性分析,獲取軌控發(fā)動機單閥推力偏差概率曲線,并計算不同容許偏差下推力偏差可靠度?；诖砟Ｐ蛯腆w姿軌控發(fā)動機推力偏差進行不確定性分析,有望為快速量化評估實時調(diào)節(jié)過程的推力輸出不確定性提供參考,同時對固體姿軌控發(fā)動機性能可靠性設計具有一定理論指導意義。

1 軌控發(fā)動機單閥推力建模

固體姿軌控發(fā)動機一般由姿控發(fā)動機和軌控發(fā)動機組成,其中姿控發(fā)動機為飛行器姿態(tài)控制提供動力,軌控發(fā)動機實現(xiàn)飛行器軌道變換。在軌控發(fā)動機工作過程中,通過多閥協(xié)同調(diào)節(jié),控制各閥門推力,可以合成不同大小和方向的合推力矢量,實現(xiàn)飛行器姿態(tài)和軌道控制。因此,軌控發(fā)動機推力快速精確調(diào)節(jié)是實現(xiàn)飛行器高精度目標控制的基礎[17],若實際推力與名義推力間的偏差超過控制系統(tǒng)容許偏差,會導致固體姿軌控發(fā)動機脫離預定軌道,無法完成預期任務。

1.1 軌控發(fā)動機單閥推力計算模型

軌控發(fā)動機單閥推力計算公式為[18]

(1)

(2)

式中ρ為推進劑密度;C*為特征速度;a為燃速系數(shù);Ab為藥柱燃面面積;At為噴管等效喉部面積;n為壓強指數(shù)。

由式(2)可以看出,燃燒室壓強主要受到推進劑材料屬性、裝藥尺寸和等效喉部面積三方面因素影響。因此,分別從三方面選取共計六種影響因素對喉栓式固體軌控發(fā)動機單閥推力進行分析計算。推進劑材料屬性相關因素為隨機性較強的推進劑特征速度和推進劑密度;由于采用端面燃燒裝藥,裝藥尺寸方面選擇藥柱直徑;等效喉部面積影響因素有喉栓位移、喉栓直徑和噴管喉徑。六種影響因素信息列于表1。

表1 軌控發(fā)動機單閥推力影響參數(shù)Table 1 Influence parameters of single valve thrust in the divert motor

計算得到燃燒室壓強后,需要對固體軌控發(fā)動機進行內(nèi)流場仿真,進而根據(jù)仿真結(jié)果計算發(fā)動機推力。圖1為固體姿軌控發(fā)動機噴管構型及裝藥示意圖。由于喉栓和噴管均為回轉(zhuǎn)體,為此采用二維軸對稱原理簡化噴管構型以提高計算效率。簡化后噴管構型如圖2所示。

圖1 軌控發(fā)動機噴管構型及裝藥示意圖Fig.1 Nozzle configuration and propellant charging schematic of divert motor

圖2 軌控發(fā)動機簡化噴管構型Fig.2 Simplified nozzle configuration of divert motor

1.2 軌控發(fā)動機單閥推力內(nèi)流場仿真方法

固體姿軌控發(fā)動機工作過程中流體滿足連續(xù)介質(zhì)假設,遵循流體力學基本方程:質(zhì)量守恒方程、動量守恒方程和能量守恒方程[19]。因此,建立流體控制方程組[20]:

(1)連續(xù)方程

(3)

(2)動量方程

(4)

(3)能量方程

·uiτji+pujfj+SE

(5)

采用Realizablek-ε兩方程湍流模型對湍動能k和耗散速度ε進行求解。仿真計算完成后,將結(jié)果代入式(1)即可計算獲取仿真推力。

為驗證本文仿真方法的準確性和有效性,建立文獻[21]中噴管構型CFD仿真模型,分別與低壓(4 MPa)和高壓(12 MPa)兩種不同工況的試驗結(jié)果進行對比,結(jié)果如表2所示。高壓狀態(tài)噴管流場云圖如圖3所示。由表2可見,低壓和高壓狀態(tài)下仿真推力與試驗結(jié)果相對誤差均在5%以內(nèi),表明本文所用數(shù)值仿真分析方法精度滿足實際工程分析需求。

表2 仿真推力與試驗推力結(jié)果對比[21-22]Table 2 Comparisons of simulation thrust and test thrust[21-22]

(a)Pressure contour

1.3 軌控發(fā)動機單閥推力仿真計算

喉栓式軌控發(fā)動機的燃燒室壓強隨喉栓位移變化,因此采用壓力入口邊界條件,總溫為3200 K。噴管和喉栓等壁面采用絕熱無滑移壁面邊界條件,其型面根據(jù)輸入條件可變化。將出口邊界條件定義為壓力出口,環(huán)境壓強及溫度設置為海平面標準大氣參數(shù)。

計算域網(wǎng)格劃分如圖4所示。軌控發(fā)動機噴管及喉栓型面變化復雜,因此采用三角形非結(jié)構網(wǎng)格,并對物理量變化較劇烈的區(qū)域,如邊界層、噴管喉部及喉栓頭部等進行網(wǎng)格加密處理。

圖4 噴管流場計算域網(wǎng)格Fig.4 Grid of calculation region of nozzle flow field

網(wǎng)格無關性驗證結(jié)果列于表3。當網(wǎng)格數(shù)量達到55 223后,仿真推力計算結(jié)果受網(wǎng)格密度影響較小。因此本文計算域網(wǎng)格數(shù)量選擇55 223,推力仿真結(jié)果云圖如圖5所示,計算得到仿真推力F0=7.863 5 kN。

表3 CFD仿真模型的網(wǎng)格無關性驗證Table 3 Verification of grid-independence of the CFD simulation model

本文建立軌控發(fā)動機單閥推力計算模型,選擇推力影響參數(shù)進行仿真計算。將計算推力作為名義推力,為軌控發(fā)動機單閥推力偏差代理模型構建提供基礎仿真數(shù)據(jù)。

2 軌控發(fā)動機單閥推力偏差代理模型構建

不確定性分析需要足夠多的仿真數(shù)據(jù)支撐以獲得精確結(jié)果,而軌控發(fā)動機推力仿真單次耗時約為0.5 h,直接使用仿真模型進行不確定性分析會造成難以承受的計算負擔。為此,通過少數(shù)仿真計算構建單閥推力偏差代理模型代替耗時仿真模型,以減少計算量。由于Kriging模型具有較為良好綜合性能[23],本文選擇Kriging模型近似軌控發(fā)動機單閥推力偏差。

2.1 Kriging模型原理

Kriging模型是基于隨機過程的半?yún)?shù)化模型,可以實現(xiàn)空間分布參數(shù)的無偏估計,僅需估計點附近的局部信息就可以給出全局近似[23]。假設有n維樣本x=[x1,x2,…,xn],其對應的響應值為Y=[Y1(x),Y2(x),…,Yn(x)],則Kriging模型一般形式可表示為[23]

g(x)=fT(x)β+z(x)

(6)

式中f(x)=[f1(x),f2(x),…,fm(x)]T為回歸模型基函數(shù);β=[β1,β2,…,βm]T為回歸系數(shù)向量,m為基函數(shù)的個數(shù);fT(x)β為線性回歸部分;z(x)為均值為0、方差為σ的正態(tài)分布隨機過程,表示模型對線性回歸部分的局部隨機偏差。

不同樣本點之間的協(xié)方差由式(7)計算:

cov(z(xi),z(xj))=σ2Rθ(xi,xj)

(7)

式中i,j=1,2,…,n;Rθ(xi,xj)為樣本點xi與xj之間的相關函數(shù);θ為相關性參數(shù)。

根據(jù)最小二乘法,未知系數(shù)β和隨機過程方差σ2可分別估計為[24]

(8)

(9)

式中F為由Fik=fk(xi)組成的矩陣(i=1,2,…,n;k=1,2,…,m);R為Ri,j=Rθ(xi,xj)組成的矩陣。

Kriging模型預測函數(shù)為

(10)

式中rT(x)=[R(x,x1),R(x,x2),…,R(x,xn)]。

2.2 軌控發(fā)動機單閥推力偏差Kriging模型構建

軌控發(fā)動機單閥推力受到噴管構型等幾何尺寸以及推進劑屬性影響,而這些因素不可避免地存在不確定性,將推力影響參數(shù)考慮為隨機變量更加符合工程實際。因此,將表1中參數(shù)設定為服從正態(tài)分布的隨機變量[20,25],則其概率分布信息列于表4。

表4 單閥推力不確定參數(shù)統(tǒng)計信息Table 4 Statistical data of uncertain parameters for the single valve thrust

考慮表4中隨機變量,實際推力與名義推力之間的偏差ε可表示為

ε(x)=F(x)-F0

(11)

式中F(x)為軌控發(fā)動機單閥實際推力;x為隨機輸入?yún)?shù)組成的六維向量;F0為名義推力,根據(jù)軌控發(fā)動機單閥推力偏差仿真計算結(jié)果,取F0=7.863 5 kN。

推力Kriging預測值為

(12)

針對文中的六維問題,一次實驗設計法采用60組初始樣本可以構建較為精確的代理模型[26]。為保證軌控發(fā)動機單閥推力偏差代理模型的精度,根據(jù)表4所示推力影響因素統(tǒng)計信息生成110組隨機樣本,通過CFD獲得仿真推力值,由輸入?yún)?shù)和對應仿真推力組成110組樣本,從中選取90組構建樣本用以建立單閥推力偏差Kriging模型,20組測試樣本來驗證代理模型精度。

四種常用代理模型精度評價指標[23]分別為復相關系數(shù)(R2)、均方根誤差(eRMSE)、相對平均絕對誤差(eRAAE)和相對最大絕對誤差(eRMAE),其表達式如下[27]:

(13)

(14)

(15)

(16)

σSTD表達式為

(17)

四種評價指標中,eRMSE和eRAAE用于評估模型全局精度,其值越接近0表明模型全局精度越高;eRMAE用來評估模型局部精度,其值越接近0,模型局部精度越高;R2用于評估模型全局精度,其值越接近于1表明模型精度越高。20組測試樣本上的評價指標值列于表5。

表5 單閥推力偏差Kriging模型精度評價指標Table 5 Evaluation index of Kriging model accuracy for single valve thrust deviation

由表5可知,所建立單閥推力偏差Kriging模型具有較高的全局精度,而其局部精度評價指標eRMAE值為0.935 2。為此,根據(jù)式(12)計算測試樣本處推力Kriging預測值,與對應樣本處CFD模擬值進行對比并求出相對誤差,推力預測相對誤差如圖6所示。由圖6可知,20組測試樣本中僅2組樣本推力預測值和模擬值的相對誤差大于0.4%,說明所建立軌控發(fā)動機單閥推力偏差Kriging模型具有足夠高的局部精度。

圖6 Kriging預測推力和CFD模擬推力相對誤差圖Fig.6 Relative error of Kriging predicted thrust and CFD simulated thrust

本文考慮參數(shù)不確定性對軌控發(fā)動機單閥推力進行仿真計算,獲得其偏差響應樣本。使用隨機輸入?yún)?shù)和推力偏差響應樣本構建推力偏差Kriging模型并驗證其精度,數(shù)據(jù)表明,所建代理模型具有較高的全局精度和局部精度。

3 軌控發(fā)動機單閥推力性能分析

3.1 軌控發(fā)動機推力影響因素敏感性排序

為研究軌控發(fā)動機推力對各影響因素的敏感性,將各參數(shù)分別設為服從正態(tài)分布的隨機變量,其他5個變量為固定值,通過單閥推力偏差的標準差度量推力對各影響參數(shù)變化的敏感程度。推力偏差的標準差越大,表示推力對相應隨機參數(shù)的變化越敏感。當?shù)趉個因素為隨機變量時,相應的軌控發(fā)動機單閥推力可表示為

(18)

式中χk=(μ1,…,μk-1,xk,μk+1,…,μ6),μk為第k個因素的均值,χk為第k個隨機變量。

在建立的軌控發(fā)動機單閥推力偏差Kriging模型基礎上,采用蒙特卡洛模擬法(MCS)獲取推力偏差的標準差,計算結(jié)果列于表6。由表6可知,當取喉栓直徑為隨機變量時,推力偏差的標準差最大,達到0.016 3 kN。取藥柱直徑為隨機變量時的推力偏差標準差最小,為0.000 8 kN。這表明,喉栓直徑對軌控發(fā)動機單閥推力偏差不確定性水平的影響最大,而藥柱直徑則導致推力偏差的最低不確定性水平。

表6 不同不確定變量下推力偏差的標準差計算結(jié)果Table 6 Results of standard deviation for thrust deviation with different uncertain parameters

為更加直觀地比較各參數(shù)對推力偏差的影響,圖7給出推力偏差敏感性分析結(jié)果柱狀圖?？芍?各參數(shù)對推力影響程度排序:喉栓直徑>喉栓位移>噴管喉徑>推進劑密度>推進劑特征速度>藥柱直徑。

圖7 推力偏差敏感性結(jié)果Fig.7 Sensitivity result of thrust deviation

由于等效喉部面積Ar對燃燒室壓強pc和固體軌控發(fā)動機推力都存在影響,表明等效喉部面積不確定性對推力不確定性影響效果更顯著。六種推力影響因素中,喉栓直徑、喉栓位移和噴管喉徑三因素與等效喉部面積相關,推進劑密度、推進劑特征速度和藥柱直徑三參數(shù)與燃燒室壓強相關。因此,喉栓直徑、喉栓位移和噴管喉徑對推力偏差影響程度大于推進劑密度,推進劑特征速度和藥柱直徑對推力偏差影響程度。

3.2 軌控發(fā)動機單閥推力偏差不確定性分析

表7給出了考慮參數(shù)不確定性時110組CFD模擬推力和105組Kriging預測推力最大值、最小值和平均值。

表7 考慮參數(shù)不確定性的軌控發(fā)動機單閥推力值Table 7 Single valve thrust values of divert motor considering parameter uncertainties

由表7可知,考慮不確定性的推力值和推力名義值相比有較大波動范圍,表明固體姿軌控發(fā)動機推力影響因素的不確定性對推力性能有不容忽視的影響。因此,需要開展軌控發(fā)動機單閥推力偏差不確定性分析,進一步探明參數(shù)不確定性對固體姿軌控發(fā)動機推力性能的影響規(guī)律。

在軌控發(fā)動機單閥推力偏差Kriging模型的基礎上,采用MCS方法對軌控發(fā)動機單閥推力偏差進行不確定性分析。110組CFD模擬推力偏差和105組Kriging預測推力偏差的前兩階矩列于表8。

表8 軌控發(fā)動機單閥推力偏差前兩階矩Table 8 The first-two order moments for single valve thrust deviation of divert motor

由于CFD數(shù)值模擬的樣本量限制,其結(jié)果僅具有一定的參考意義。兩種方法得到的推力偏差均值和標準差都相差不大,說明Kriging預測推力偏差不確定性分析結(jié)果精度較高。

圖8為軌控發(fā)動機單閥推力偏差概率密度函數(shù)(Probability Density Function,PDF)和累積分布函數(shù)(Cumulative Distribution Function,CDF)曲線。由圖8可知,在整體趨勢上,Kriging預測值的概率曲線和CFD仿真值概率曲線擬合較好。出于計算成本考慮,本文僅獲取110組CFD模擬數(shù)據(jù),難以獲得準確推力偏差概率信息,導致其PDF和CDF曲線有較大波動。因此,以105組樣本所得Kriging預測值的概率曲線分析固體軌控發(fā)動機單閥推力偏差不確定性。

(a)PDF curves of the thrust deviation

(b)CDF curves of the thrust deviation圖8 軌控發(fā)動機單閥推力偏差概率曲線Fig.8 Probability curves of single valve thrust deviation for divert motor

當軌控發(fā)動機實際推力與名義推力之間偏差小于一定范圍時,可由控制系統(tǒng)對其推力進行補償,以保證發(fā)動機實現(xiàn)預期軌道控制效果,此范圍稱為控制系統(tǒng)容許偏差。當推力偏差小于容許偏差時,認為軌控發(fā)動機推力偏差性能可靠[25]軌控發(fā)動機推力偏差性能可靠的概率稱為軌控發(fā)動機推力偏差可靠度。

表9為容許偏差分別取0.01～0.1 kN時的推力偏差可靠度。當容許偏差取0.1 kN時,推力偏差CFD仿真值與Kriging預測值的可靠度分別為99.12%和99.56%,表明此時固體姿軌控發(fā)動機控制系統(tǒng)可以補償幾乎全部推力偏差,使其達到理想控制效果。當軌控發(fā)動機單閥推力偏差可靠度設計要求為99%時,控制系統(tǒng)對發(fā)動機推力輸出的容許偏差為0.09 kN即可實現(xiàn)該設計要求。

表9 不同容許偏差下單閥推力偏差可靠度Table 9 Reliability of single valve thrust deviation subject to different tolerable deviations

為更直觀地比較不同容許偏差下軌控發(fā)動機單閥推力偏差可靠度,圖9給出對應到單閥推力偏差PDF曲線上的推力偏差可靠域,其中灰色區(qū)域面積表示單閥推力偏差可靠度。

(a)Tolerable deviations is 0.02 kN

(b)Tolerable deviations is 0.04 kN

(c)Tolerable deviations is 0.06 kN圖9 不同容許偏差下單閥推力偏差可靠域Fig.9 Reliability region of single valve thrust deviation subject to different tolerable deviations

對固體軌控發(fā)動機推力偏差進行不確定性分析,探明發(fā)動機推力調(diào)控過程中參數(shù)不確定性的影響,可為固體姿軌控發(fā)動機推力快速準確調(diào)節(jié)提供參考。同時,可以更加高效并直觀地了解該發(fā)動機推力性能是否滿足可靠性設計要求,為喉栓式固體軌控發(fā)動機設計提供理論指導。

4 結(jié)論

考慮影響固體姿軌控發(fā)動機推力性能的不確定性因素,基于CFD數(shù)值模擬構建軌控發(fā)動機單閥推力偏差Kriging模型,探究不確定性對固體姿軌控發(fā)動機推力偏差的影響。主要結(jié)論如下:

(1)根據(jù)軌控發(fā)動機推力調(diào)控原理和計算模型,將喉栓位移、喉栓直徑、噴管喉徑、藥柱直徑、推進劑密度和特征速度等參數(shù)的不確定性視為影響軌控發(fā)動機推力性能的主要因素。

(2)采用復相關系數(shù)、均方根誤差、相對平均絕對誤差和相對最大絕對誤差四種評價指標對本文構建的代理模型進行精度驗證,數(shù)據(jù)表明所建立的推力偏差Kriging模型具有較高的全局精度和局部精度。

(3)依據(jù)所建立的Kriging模型開展軌控發(fā)動機推力影響因素敏感性、不確定性和可靠性分析。結(jié)果表明,為提高固體軌控發(fā)動機單閥推力偏差可靠度,在控制系統(tǒng)容差確定后,需著重控制喉栓直徑、喉栓位移和噴管喉徑等噴管構型相關參數(shù)不確定性水平,同時有選擇地控制推進劑密度、特征速度和藥柱直徑等裝藥相關參數(shù)不確定性水平以平衡成本。