喉栓式固體軌控發(fā)動機單閥推力偏差不確定性分析①

張德權(quán),賈軍凱,武澤平,王東輝,牛 祿,周偉華

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

0 引言

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

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

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

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

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)飛行器高精度目標(biāo)控制的基礎(chǔ)[17],若實際推力與名義推力間的偏差超過控制系統(tǒng)容許偏差,會導(dǎo)致固體姿軌控發(fā)動機脫離預(yù)定軌道,無法完成預(yù)期任務(wù)。

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

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

(1)

(2)

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

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

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

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

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

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

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

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

(1)連續(xù)方程

(3)

(2)動量方程

(4)

(3)能量方程

·uiτji+pujfj+SE

(5)

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

為驗證本文仿真方法的準(zhǔn)確性和有效性,建立文獻[21]中噴管構(gòu)型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è)置為海平面標(biāo)準(zhǔn)大氣參數(shù)。

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

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

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

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

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

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

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

2.1 Kriging模型原理

Kriging模型是基于隨機過程的半?yún)?shù)化模型,可以實現(xiàn)空間分布參數(shù)的無偏估計,僅需估計點附近的局部信息就可以給出全局近似[23]。假設(shè)有n維樣本x=[x1,x2,…,xn],其對應(yīng)的響應(yīng)值為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之間的相關(guān)函數(shù);θ為相關(guān)性參數(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模型預(yù)測函數(shù)為

(10)

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

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

軌控發(fā)動機單閥推力受到噴管構(gòu)型等幾何尺寸以及推進劑屬性影響,而這些因素不可避免地存在不確定性,將推力影響參數(shù)考慮為隨機變量更加符合工程實際。因此,將表1中參數(shù)設(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預(yù)測值為

(12)

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

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

(13)

(14)

(15)

(16)

σSTD表達式為

(17)

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

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

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

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

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

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

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

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

(18)

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

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

表6 不同不確定變量下推力偏差的標(biāo)準(zhǔn)差計算結(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ā)動機推力都存在影響,表明等效喉部面積不確定性對推力不確定性影響效果更顯著。六種推力影響因素中,喉栓直徑、喉栓位移和噴管喉徑三因素與等效喉部面積相關(guān),推進劑密度、推進劑特征速度和藥柱直徑三參數(shù)與燃燒室壓強相關(guān)。因此,喉栓直徑、喉栓位移和噴管喉徑對推力偏差影響程度大于推進劑密度,推進劑特征速度和藥柱直徑對推力偏差影響程度。

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

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

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

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

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

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

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

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

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

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

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

為更直觀地比較不同容許偏差下軌控發(fā)動機單閥推力偏差可靠度,圖9給出對應(yīng)到單閥推力偏差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ā)動機推力快速準(zhǔn)確調(diào)節(jié)提供參考。同時,可以更加高效并直觀地了解該發(fā)動機推力性能是否滿足可靠性設(shè)計要求,為喉栓式固體軌控發(fā)動機設(shè)計提供理論指導(dǎo)。

4 結(jié)論

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

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

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

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