捆綁火箭氣動載荷分布不確定性分析

董欣心，劉莉，葛佳昊，王志

（北京理工大學 宇航學院，北京 100081）

火箭飛行過程中承受靜態作用力和動態作用力，前者產生靜載荷，后者產生動載荷。靜載荷主要包括氣動載荷、發動機產生的控制力、液體推進劑產生的晃動力等；動載荷主要包括陣風載荷、抖振載荷等［1］。火箭氣動學科是總體設計的重要專業之一，氣動特性分析得到的各項氣動系數是載荷、姿態控制、結構設計等學科的先決條件，在火箭小回路設計中起到重要作用。然而在火箭飛行過程中包含許多不確定性因素，這種不確定性因素既有來自模型本身的不確定性，如在生產前無法完全準確獲知結構尺寸及機械強度［2］，也有由于對載荷條件知識缺乏造成的來自自然界的不確定性，如外界風場對于載荷的影響。準確的載荷分析是火箭結構設計的先決條件，優良的設計應能保證在任何外載荷的激勵下都能可靠完成任務，因此，需對載荷的不確定性展開分析。

傳統的火箭設計過程中，對于不確定性因素常采用計入安全系數的方式進行考慮，這種方式的問題在于安全系數的選取對于設計人員經驗的依賴性較強，由于缺乏對不確定性的量化，安全系數過小無法保證安全，而過保守主義又會提升零件生產的經濟成本，并且會帶來不必要的質量增加。

近年來，不確定性優化設計［3-4］的概念逐漸在設計領域受到重視，其前提是不確定性量化（uncertainty quantification），主要思想是在設計過程中通過對模型輸出響應的不確定性定量描述提升模型的可信度。不確定性分析可以幫助設計人員更好地量化模型響應的不確定性，獲知模型置信度。對于氣動載荷的不確定性分析主要集中在對氣動力系數的分析，包括升力系數、阻力系數及俯仰力矩系數，可以據此計算機體或箭體所受的氣動合力，為姿態控制專業提供上游數據。宋鑫等［5］對對稱翼型開展了不確定性分析，得到了翼型氣動性能定量變化區間，并以阻力系數最小為目標對其開展魯棒性優化設計。針對由馬赫數及迎角不確定性導致的氣動性能波動問題，鄔曉敬等［6］以NACA0012翼型為例，分析了跨聲速階段該現象產生的機理。

然而，僅關注氣動力系數難以為火箭結構設計提供有效的指導，在結構設計過程中，需給出箭體分站上每一截面的內力，以便對火箭各部段分別進行設計。尤其是對于捆綁火箭，芯級與助推器捆綁處存在復雜的氣動干擾問題。國內外學者對此類干擾問題開展了相關研究。Uematsu等［7］開展了高超聲速流動條件下助推器與芯級分離的風洞試驗，分別分析了不同攻角下不同截面形狀的助推器在分離過程中對芯級的干擾沖擊波形式。沈丹等［8］采用數值仿真模擬了高超聲速條件下助推器頭部對芯級的激波干擾，并對比了直頭助推與斜頭助推的干擾特性。由于捆綁火箭氣動特性存在的復雜性及不確定性，了解火箭自身參數及大氣條件等因素對氣動載荷分布不確定性的影響是十分必要的。

不確定性分析方法主要可以分為概率方法［9-11］（probabilistic methods）和可能性方法［12-13］（possibilistic methods）。傳統的概率方法多采用統計方法，如蒙特卡羅方法，其問題是所需樣本數量較大、計算量較大及計算效率低，尤其是對于氣動分析問題，單次計算花費時間較長，計算成本難以接受。廣義多項式混沌方法（general polynomial chaos，gPC）作為一種非統計的概率方法，具有所需樣本點數目少和計算精度高的優點，在氣動分析領域得到廣泛應用。Bettis等［14］依據配點非侵入式多項式混沌方法獲取隨機響應面，對運載火箭多學科分析獲得的輸出參數混合不確定性進行分析。Hosder等［15］采用稀疏的非侵入式多項式混沌方法及全局非線性靈敏度分析方法，對返回艙后體再入過程的輻射加熱現象進行分析。宋賦強等［16］采用稀疏的非侵入式多項式混沌方法量化了乘波體的氣動特性不確定性。

如果能夠獲知各項因素對響應不確定性的貢獻程度，分析不確定性的主要來源，就能通過控制這些因素提升模型的置信度，此時需要引入靈敏度分析［17］。同時，由于不確定性分析方法的計算量會隨輸入參數的增加而大幅提升，合理選定輸入參數非常重要，而靈敏度分析也為不確定分析中參數的選取提供依據。

首先，提出了依據多項式混沌理論，對捆綁火箭氣動載荷分布開展全局靈敏度分析及不確定性分析的方法。然后，采用文獻中給出氣動特性實驗數據的兩捆綁火箭模型進行CFD分析，驗證流場求解結果，并以此模型為基礎建立外形參數化模型，對提出的方法進行驗證。在來流參數及外形參數存在擾動的情況下，采用拉丁超立方采樣獲取樣本點，通過CFD仿真分析獲取樣本點壓力系數及黏性應力系數，分段積分得到氣動力系數沿箭體分布，并由此計算氣動載荷分布。最后，依據計算結果分析了流場的流動情況及氣動載荷波動的產生原因。

1 氣動載荷分布不確定性分析方法

1.1 靈敏度分析方法

常用的全局靈敏度分析方法有基本效應法［18］、基于導數的靈敏度分析法［19］、基于方差的靈敏度分析法［20-21］、矩獨立靈敏度分析法［22］等，其中方差靈敏度分析法又稱為Sobol方法，通過輸入變量的方差在系統響應方差中的貢獻度來衡量該變量的重要程度，能充分表征變量對響應的影響程度，因此，本文采用Sobol方法分析各因素對氣動載荷的靈敏度。

將模型的不確定性輸入定義為X，不確定性響應定義為Y，則該不確定模型函數可表示為Y=f（X），其中X=（X1，X2，…，Xn），為不失一般性，對各不確定性輸入進行歸一化處理，使n維輸入變量在[ 0，1]內服從均勻分布，變量間相互獨立，通過高維模型展開可以得到

可以進行展開的條件為

說明展開式中的所有項均兩兩正交，其表示為

對式（1）兩側同時取方差，得到

Sobol指數為［21］

Sobol指數的取值在［0，1］之間，且取值越大說明該輸入變量對響應的影響越大。

1.2 不確定性量化方法

采用多項式混沌展開方法對系統進行不確定性分析，其思想在于將一個隨機過程利用屬于某特定分布類型的正交多項式混沌之和替代。假設有隨機過程如下：

式中：ξ=［ξ1，ξ2，…，ξn］為服從一定分布的隨機變量，并且相互獨立。

則可將式（6）展開，并在有限階截斷得到

其中：n為隨機變量數量；s為混沌多項式階數。

依據Askey準則，正交多項式基底的選擇取決于標準隨機變量ξ的分布，自然界中大部分隨機現象可以用高斯分布來描述，本文問題中討論的不確定性變量也服從高斯分布，其對應的多項式為Hermite多項式。n維Hermite多項式表達式為

根據多項式混沌理論，得到原不確定性系統的前二階統計矩為

1.3 載荷計算方法及仿真流程

表妹和表妹夫這次徹底把心結打開，兩人有所反思，他們約定為了共同的未來一起努力付出。經過調解他們的婚姻危機，喬振宇似有所悟，他不再跟我陰陽怪氣地對話，份內的家務活也重新上手了。我們倆回到了以前的生活軌道上，可我覺得僅僅這樣是不夠的，這樣只能讓我們暫時回到貌似風平浪靜的婚姻狀態，可是，他的大男子主義不連根拔出，一旦出現新問題，我們的關系勢必回到抱怨爭執的狀態，繼續彼此傷害繼續惡化，到那時候任憑多厚的愛情基礎也扛不住這“強拆”啊！

在實際計算過程中，對式（12）沿軸向按氣動分站進行分段積分，得到每一站點的氣動力系數，近似代替該段的氣動力系數。由此結果可進一步計算得到沿箭體x軸分布的軸力及法向力為

仿真計算流程如圖1所示，主要步驟如下：

圖1 仿真分析流程Fig.1 Flowchart of simulation and analysis

步驟1 依據多項式混沌展開階數確定采樣數量N，依據變量分布，采用拉丁超立方采樣獲得自變量X的N組樣本點。

步驟2 基于VBA語言對SolidWorks二次開發，實現參數化模型自動更新，完成對結構外形參數的更新，導出STEP格式文件。

步驟3 采用CFD前處理軟件Gambit對模型進行氣動網格劃分及邊界條件設置，依據預先錄制的Journal文件讀取步驟2生成的STEP格式模型文件，進行自動化網格劃分，并生成網格Mesh文件。

步驟4 依據預先錄制的Journal文件實現CFD仿真條件設置，讀取步驟3生成的Mesh文件，采用MATLAB腳本更改Journal文件，實現來流參數的自動更新。

步驟6 對求得的氣動特性數據開展Sobol靈敏度分析及不確定性分析。

2 模型建立及流場求解

2.1 捆綁火箭模型參數化

為驗證仿真分析結果的可靠性，以文獻［23］中提出的兩捆綁構型火箭模型為對象開展氣動載荷不確定性分析，其外形如圖2所示，具體參數數值在表1中給出。該模型平均了現有主要帶捆綁運載火箭的外形尺寸參數，且文獻中給出了該模型芯級表面壓力系數分布的風洞試驗數據。采用SolidWorks建模，并通過Excel中的VBA程序進行二次開發，可以實現外形自動化更新。

圖2 捆綁火箭參數化模型Fig.2 Strap-on launch vehicle parametric model

表1 外形參數數值Table 1 Values of shape parameters

2.2 網格劃分及流場條件設置

由于捆綁式運載火箭的外形較復雜，芯級和助推器之間存在小的間隙，難以采用結構網格來生成計算域的網格，模型氣動網格生成采用非結構網格、圓柱形的計算域。首先，生成火箭壁面和圓柱面邊界上的網格點，連接生成三角形網格。然后，采用陣面推進法生成四面體非結構網格。

為了獲取更準確的計算結果，在火箭頭部整流罩處和芯級與助推器間易產生氣流干擾處進行局部氣動網格加密。計算域的選擇應保證流動能夠充分發展，計算域選定為拋物面圍成區域，長為90 m，底面半徑為25 m，箭體頭部距拋物面頭部10 m。為降低計算量，用對稱面將模型分為2部分，僅對半模型進行分析計算。最終生成四面體網格數量為793 678，箭體、計算域表面網格及邊界條件設置如圖3所示。

圖3 網格劃分示意圖Fig.3 Schematic diagram of grid generation

模型驗證工況及分析工況的流動條件均為高度可壓縮流動，對于此類三維定常可壓縮流動問題，選用密度基求解器能夠獲得更好的求解效果。采用基于節點的Green-Gauss函數求梯度，該方法精度較高，適用于非結構化網格。假設空氣為理想氣體，服從理想氣體狀態方程。湍流模型采用Spalart-Allmaras模型，適用于具有壁面限制的流動問題，對有逆壓梯度的邊界層問題能夠求得較好的結果。

2.3 氣動特性模型驗證

為驗證網格劃分合理性與計算模型準確性，將仿真計算結果與文獻［23］中提供的風洞試驗數據進行對比。選用2種工況對計算結果進行準確性驗證。其中，工況1為Ma=0.6、α=0°，工況2為Ma=2.0、α=0°。兩助推中心軸線位于XOY平面內，假設迎角僅位于XOZ平面內（見圖4），則法向力沿z軸方向分布。

就火箭芯級靠近助推器一側，即圖4中紅色虛線所示位置的表面壓力系數Cp沿箭體軸向分布。圖5給出了本文及文獻［23］的CFD仿真結果與風洞試驗結果對比。從圖5可以看出，本文的仿真結果與風洞試驗數據吻合較好，略優于文獻［23］的仿真結果，可以滿足后續分析對模型精度的需求。

圖4 迎角方向及表面壓力系數提取位置Fig.4 Attack angle direction and surface pressure coefficient extraction location

圖5 表面壓力系數分布對比Fig.5 Comparison of surface pressure coefficient distribution

3 仿真驗證及結果分析

3.1 仿真參數設置

分析計算發現，氣動載荷變化較為劇烈的主要為火箭頭部整流罩處及捆綁與芯級發生氣流干擾處，故將不確定性輸入選為X=［θ1，d，Ma，α］，包括2個幾何外形尺寸參數，分別為整流罩頭錐半頂角θ1和芯級助推間距d，以及2個氣動參數，分別為飛行馬赫數Ma和迎角α。分析上述參數在小范圍內擾動對捆綁火箭氣動載荷的影響情況。

自然界中的變量廣泛服從正態分布，依據中心極限定理，誤差可以看作是許多量的疊加，理應服從正態分布，因此，認為選取的變量均服從正態分布。選取火箭氣動載荷分析過程中較為關注的最大動壓狀態進行分析，通過查閱相關型號彈道數據確定該工況對應的飛行狀態。假設Ma均值為1.3，95%概率位于［1.2，1.4］，即Ma～N（1.3，0.052）；α均值為5°，95%概率位于［4°，6°］，即α～N（5，0.52）；θ1均值為17°，d均值為1 000 mm，均取變異系數cv為0.06，即θ1～N（17，1.022），d～N（1 000，602）。計算采用三階多項式混沌方法，依據式（8）可得三階四維問題的展開式項數為34，采樣點數通常取項數的1.5～3倍［6］，選取采樣點為60。依據變量分布，采用拉丁超立方采樣獲得Ma、α、θ1、d的60組樣本點。

3.2 靈敏度及不確定性結果分析

在結構設計過程中，通常采用分段設計方法，選取各部段載荷最大的截面，將最大載荷作為該部段的設計載荷，因此重點關注載荷峰值。氣動載荷按沿箭體x軸方向及垂直于箭體x軸方向可以分為軸向力和法向力，依據式（13）計算得到各組樣本對應軸向力及法向力分布。

依據載荷分布結果，分別選定軸向力分布在頭錐處的峰值（峰值1）、法向力在頭錐處的峰值（峰值2）及芯級和助推發生氣流干擾處的峰值（峰值3）大小作為響應值，峰值所在位置如圖6（a）、圖6（b）中所示，分析頭錐半頂角、芯級助推間距、迎角及馬赫數對響應的靈敏度。分析得到Sobol指數如圖7所示。

圖6 軸向x及法向力不確定性分布Fig.6 Uncertainty distribution of axial force and normal force

圖7 不確定性因素Sobol靈敏度指數Fig.7 Sobol sensitivity indices of uncertainty factors

從靈敏度分析結果可以看出，飛行馬赫數對軸向力峰值的影響較大，而飛行迎角對法向力峰值有明顯影響，對軸向力峰值影響很小。頭錐半頂角對頭錐處軸向力影響較大，而對法向力影響偏小；芯級助推間距對二者發生氣動干擾處的法向力有一定影響，但其影響小于迎角及馬赫數。馬赫數對軸向力及法向力的作用包括力系數和動壓2個部分，因此，對各響應值均有較大影響。

依據多項式混沌展開方法對軸向力和法向力分布的期望和方差進行分析，得到其不確定性分布如圖7所示。圖中：實線表示軸向力和法向力分布的期望值μ，即平均水平；虛線表示偏離均值一倍標準差位置μ+σ，即對應發生概率為65.26%；藍色邊界表示偏離均值二倍標準差位置μ+2σ，即對應發生概率為95.44%，顏色越淺表示載荷到達該值的概率越低。所關注載荷峰值出現位置、對應的均值、標準差和不確定度結果如表2所示。

表2 不確定性分析結果Table 2 Uncertainty analysis results

從不確定性分析結果可以看出，軸向力的不確定性主要出現在頭錐處，峰值1處不確定度為10.90%，其后在火箭芯級截面無變化處幾乎為零，變化十分平緩，僅在出現助推器位置有微弱波動；相比之下法向力的不確定性較為顯著，在頭錐處、助推干擾處及尾段均有明顯波動，峰值2處不確定度為6.327%，峰值3處不確定度為10.79%，在芯級截面不變且未與助推器產生氣流干擾處同樣存在一定的不確定性，主要是由于來流參數波動產生的。

3.3 流場結果分析

為了更直觀地分析火箭幾何外形對氣動特性產生影響的原因，將迎角及馬赫數固定為Ma=1.3，α=5°，分別分析d=1 000 mm時θ1=16°、θ1=17°、θ1=18°三 種 工 況 及θ1=17°時d=800 mm、d=1 000 mm、d=1 200 mm三種工況下的流動情況。選取受工況改變影響較為顯著的火箭芯級靠近助推器一側的表面壓力系數分布情況進行對比，結果如圖8所示。

圖8 不同工況下表面壓力系數分布情況對比Fig.8 Comparison of pressure coefficient distribution at different conditions

從圖8（a）中可以看出，隨著頭錐半頂角增大，來流受到的壓縮程度變大，頭錐處表面壓力系數峰值出現位置逐漸提前，且峰值逐漸增大；從圖8（b）中明顯可以看到，激波與膨脹波來回反射壓縮引起的壓力系數震蕩，隨著芯級助推間距增大，干擾處表面壓力系數峰值出現位置略微后移，對應的波動峰值也逐漸減小。圖9給出標稱工況下的馬赫數云圖及不同工況下的表面壓力系數云圖。

從馬赫數及表面壓力系數云圖中可以看出，頭錐處激波和膨脹波的干擾會導致沖擊波后部產生傾斜的高流速區域，即低壓區。而隨著頭錐半頂角增大，物面前端內折角度及后部外折角度均增大，導致激波波面傾角增大，頭錐表面氣流速度降低更快，氣壓波動幅度較大。芯級和助推間存在膨脹波（由于壓縮來回反射的現象），當間距增大時，反射次數減少，表面壓力系數及載荷變化都更加平衡，能量的損失也更小。單獨從氣動設計角度來講，適當減小頭錐半頂角、增大芯級助推間距更有利于載荷平衡變化，為結構設計降低不確定性。

4 結 論

1）針對帶捆綁火箭氣動載荷分布不確定的問題，提出了基于多項式混沌理論的捆綁火箭氣動載荷靈敏度分析及不確定性量化方法，為結構設計提供更準確的需用載荷。

2）通過算例分析了來流馬赫數、迎角及2個外形參數對氣動軸向力和法向力峰值的靈敏度。結果表明，來流速度是影響軸向力的主要因素，而迎角是影響法向力的主要因素。頭錐半頂角對頭錐處的軸向力影響較大，芯級助推間距對二者產生氣流干擾處的法向力有一定影響，但相對于馬赫數和迎角影響不大。

3）通過非侵入式多項式混沌方法得到軸向力及法向力不確定性分布情況。結果表明，軸向力不確定性主要出現在頭錐處，法向力不確定性主要出現在頭錐及與助推器產生干擾處。頭錐處的載荷不確定性主要原因是來流速度變化和頭錐截面變化導致激波與膨脹波波角及強度不確定；助推干擾處的法向力不確定性主要由于來流條件及間距變化加劇了膨脹波反射波動。

4）頭錐半頂角增大、芯級助推間距減小均會加劇芯級表面壓力系數波動，增大氣動載荷變化幅度。因此，為降低結構設計過程的載荷不確定性，可對氣動外形參數設計提出一定建議，如適當減小頭錐半頂角和增大芯級助推間距。

本文采用火箭半模型進行仿真計算，僅考慮了橫向氣動載荷影響，采用全模型可以綜合分析橫側向氣動載荷不確定性及來流方向與助推所在平面夾角變化對氣動載荷的影響，后續將在本文基礎上對此問題進行進一步研究。