多工況下變推力固體發動機喉栓噴管型面一體化優化設計①

王鵬宇，王政濤，武澤平，張 錫，張為華

(國防科技大學 空天科學學院，長沙 410073)

0 引言

固體火箭發動機具有結構簡單、維護方便、可靠性高、可長期處于戰備狀態等優點，目前已在各類戰術、戰略導彈中得到廣泛應用。傳統固體火箭發動機缺乏節流措施，需要通過對裝藥構型進行詳細設計以滿足特定推力方案。喉栓式變推力發動機不僅保留了傳統固體發動機的諸多優點，而且推力調節具有主動性、實時性、可大范圍無極調節等優點，可提高固體火箭發動機性能，增強導彈機動性和突防能力，大幅提升導彈射程，已發展成為固體火箭發動機推力調節的重要研究方向。但是，喉栓運動改變噴管喉部面積實現推力調控的同時，不可避免地會影響噴管內燃氣的流動情況，造成流動損失，甚至出現流動分離現象。BURROUGHS等對不同形狀喉栓的變推力發動機進行試驗和數值分析，發現在喉栓位置下游存在微弱激波。OSTRANDER等通過仿真分析喉栓介入對發動機流場的影響，結果表明，當喉栓沿軸線運動時，燃氣流動方向偏離發動機軸線，造成一定的推力損失，且在喉栓頭部及噴管擴張段內會出現流動分離現象，加重該區域內的燒蝕現象。因此，需要對噴管型面及喉栓型面進行整體優化，以減小喉栓介入所引起的流動損失。成沉等基于響應面法對喉栓式變推力發動機進行一體化優化設計，使得喉栓完全介入狀態下比沖損失最小。由于變推力發動機實際工作過程中喉栓需在一定范圍內移動，僅針對喉栓某特定位置進行型面優化無法綜合考量發動機的性能。因此，變推力發動機優化過程需喉栓不同位置的進行多工況仿真計算。

本文基于Fluent進行喉栓式噴管多工況流場仿真計算，通過對噴管型面及喉栓型面控制參數進行優化設計，以期實現喉栓式噴管性能的提升。由于CFD模型求解耗時，直接進行數值優化設計計算成本難以接受。因此，本文采用基于代理模型的優化框架，實現型面控制參數快速優化，以滿足工程優化設計效率需求。

1 優化設計模型

1.1 型面參數化

本文優化的噴管收斂段、喉部半徑保持不變，而擴張段采用圓弧及三次曲線描述，設計變量包括擴張段過渡圓弧半徑、擴張段過渡圓弧角、噴管出口傾角、噴管擴張段長度以及噴管擴張比。其中，噴管擴張比由式(1)給出：

(1)

式中為噴管出口半徑；為噴管喉部半徑。

為保證喉栓運動過程中所形成的等效喉部面積最大、最小值相同，喉栓半徑保持不變，而喉栓頭部型面則采用雙圓弧切線構型，由喉栓頭部長度、喉栓頭體過渡圓弧半徑以及喉栓頭部圓弧半徑等設計變量唯一確定，噴管和喉栓構型示意如圖1所示。

圖1 噴管及喉栓型面示意圖Fig.1 The contour of the pintle nozzle

1.2 數值計算模型

本文采用Favre平均N-S方程，其在雷諾平均N-S (RANS)方程組的基礎上提出密度加權平均概念，即壓強和密度采用時間平均，其他變量采用密度加權平均。FANS方程組解決了RANS方程組處理可壓縮流動較為復雜的問題，通適于不可壓縮流動與可壓縮流動。

二維軸對稱的FANS方程組的表達式為

(2)

式中為流場守恒變量矢量；為對流通量矢量；為粘性擴散通量矢量。

其中，散度的表達式為

(3)

Favre平均N-S方程組中包含湍流粘性系數與湍動能，需通過求解湍流模型得到。本文采用Realizable-湍流模型，通過偏微分方程處理湍動能及其耗散速度。

邊界條件是流場仿真計算的關鍵因素，本文設置入口邊界條件為壓力入口，不同于普通固體火箭發動機的定值壓力入口，喉栓式變推力發動機喉栓運動過程中會引起燃燒室壓力變化。因此，邊界入口壓強隨喉栓的運動發生變化：

(4)

式中為推進劑密度；為特征速度；為燃速系數；為燃面面積；為等效喉部面積；為壓強指數。

出口邊界定義為壓力出口，分別設定環境壓強及環境溫度；噴管、喉栓等壁面采用增強壁面方式處理。

由于噴管及喉栓型面變化較為復雜，本文采用非結構網格進行網格劃分，并對壁面處進行適當加密處理，如圖2所示。網格質量檢測顯示，所劃分網格正交性良好，滿足計算要求。網格尺度劃分驗證結果如表1所示。可以看出，當網格數量大于60 000時，流場結構變化較小。因此，選取60 000網格進行計算。

圖2 計算域網格劃分Fig.2 The grid of calculation domain

表1 網格無關性驗證

2 基于代理模型的優化方法

基于代理模型的優化方法在20世紀80年代被首次提出，并成功應用于計算耗時的復雜模型優化問題。該方法基于少量的高精度采樣點建立近似模型，用以在優化過程中替代復雜的高精度學科模型，通過對近似模型的搜索確定新的采樣點并不斷重構更新近似模型，以期逐漸逼近真實學科模型，預測全局最優解。基于代理模型的優化可以大幅減少真實模型的調用次數，從而提高優化設計效率，降低設計成本。

其主要步驟包括實驗設計采樣、構造代理模型以及序列近似采樣等。基于代理模型的喉栓式變推力發動機型面優化流程如圖3所示。

2.1 試驗設計

實驗設計是基于代理模型優化方法的關鍵步驟，直接影響著初始代理模型的質量，合理的實驗設計方法可以獲得更優的采樣方案，以得到更多的空間信息。本文采用優化拉丁超立方實驗設計方法，在設計空間內生成均勻的初始樣本點。

2.2 Kriging模型

代理模型基于已有的樣本點，采用簡單的數學函數來對計算耗時的高精度仿真模型進行逼近，從而實現快速地對全設計空間函數值進行預測。龍騰等將目前使用較為廣泛的多項式響應面模型、徑向基函數模型、Kriging模型以及人工神經網絡進行綜合分析比較，發現Kriging模型在近似精度、魯棒性等方面有著更好的綜合性能。因此，采用Kriging模型對CFD真實模型進行逼近。

Kriging模型是一種基于隨機過程的最優線性無偏估計。該模型通過高斯隨機過程描述未知函數，并通過已知樣本點對未知函數的響應進行預測，即

(5)

其中，為未知常數，表征變量域的整體變化趨勢；()為局部隨機偏差項，其均值為零，方差為。因此，設計空間內兩點、的協方差可表示為

[(),()]=(,)

(6)

基于無偏性以及預測方差最小準則，可求得模型預測函數為

(7)

式中∈為個樣本點響應值組成的維列向量；()∈為預測點與樣本點的協方差；∈×為樣本點協方差矩陣；∈為單位列向量。

此外，Kriging模型還能夠給出預測方差：

(8)

圖3 型面優化流程圖Fig.3 Optimization design flow chart for the contour

2.3 優化算法及序列采樣準則

在基于代理模型的優化過程中，優化算法以特定采樣準則為目標函數，在代理模型上搜索最優解作為新的樣本點用以更新近似模型，稱為“輔助優化”過程。本文在“輔助優化”過程中采用增強精英保留的多種群協同遺傳算法，通過為不同的種群設置不同的交叉和變異概率豐富總群個體，并基于概率進行中群間的遷移，實現高效的全局收斂。

序列采樣準則是基于代理模型優化過程的關鍵步驟之一，它直接影響著優化效率和結果的準確性。目前，工程實際中常用的采樣準則為開發準則，即將代理模型的全局最優點作為新樣本點來更新近似模型。開發準則忽略近似模型的全局精度，旨在最有潛力的區域進行采樣，以加快算法的收斂速度。但由于樣本量較少，代理模型在優化初期通常精度較低，對不精確的代理模型進行過渡的搜索會導致計算資源浪費，甚至可能將搜索進程引向前景較差的區域，陷入局部最優。

圖4 基于不精確搜索的采樣過程示意圖Fig.4 Schematic diagram of sampling process based on inaccurate search

針對上述問題，學者將非精確搜索算法及外部精英庫引入“輔助優化”過程，采用精英庫動態存儲所有精英個體，當發現比精英庫最優解更為優秀的個體時，“輔助優化”過程立即停止，并采用高精度仿真模型計算該點的響應值。

(9)

其停止準則為

(10)

非精確搜索能夠逐步地引導搜索過程朝著好的方向發展，可以有效陷入局部最優。因此，本文基于非精確搜索的開發準則序列填充新的樣本點。

3 優化算例及結果分析

本文對已有喉栓式變推力發動機型面進行優化設計，基準構型設計變量及其上下限如表2所示。為衡量喉喉栓式噴管多工況條件下的實際性能，基于Fluent分別計算喉栓未介入狀態(下稱工況1)、喉栓部分介入狀態(指等效喉部面積減小50%，下稱工況2)以及喉栓完全介入狀態(下稱工況3)下噴管的實際比沖，并依據其工作時間進行加權平均，以此平均實際比沖作為衡量噴管性能的指標，即

(11)

式中、分別為喉栓未介入狀態及完全介入狀態下的實際比沖；、分別為兩種狀態下的工作時間；為喉栓部分介入狀態下的實際比沖；為該狀態下的工作時間。

采用優化拉丁超立方實驗設計方法生成20個均勻分布的初始樣本點，對其進行CFD仿真計算后構建Kriging模型。基于2.3節所述“輔助優化”過程新增樣本點150個，其收斂曲線如圖5所示。

表2 設計變量取值范圍及基準構型

圖5 收斂曲線Fig.5 The convergence curve

圖6為優化前后噴管型面對比，圖7為優化前后喉栓型面對比。變推力固體發動機喉栓噴管型面一體化優化設計結果如表3所示。通過對相關型面參數的優化，喉栓式噴管的由2264.536(N·s)/kg增加至2312.781(N·s)/kg，提升2.13%。

對比圖8和圖9優化前后流場的馬赫數云圖可以看出，優化后的型面使得喉栓頭部下游的激波明顯減弱，亞音速區域減小，這有助于減小流動過程中的損失，提高實際比沖。

圖6 優化前后噴管型面對比Fig.6 Comparison of nozzle between benchmark shape and optimized shape

圖7 優化前后喉栓型面對比Fig.7 Comparison of pintle between benchmark shape and optimized shape

表3 優化結果對比

(a) Condition 1: the pintle is not inserted

(b) Condition 2: the pintle is partial inserted

(c) Condition 3: the pintle is fully inserted

(a) Condition 1: the pintle is not inserted

(b) Condition 2: the pintle is partial inserted

(c) Condition 3: the pintle is fully inserted

圖10為噴管出口截面的徑向速度分布示意圖。某種程度上，發動機的比沖損失與徑向速度成正相關。由圖10可知，優化后型面的徑向速度在三種工況下都明顯減小。因此，優化后的型面能夠減小流動損失，獲得更高的平均實際比沖。

圖10 噴管出口徑向速度Fig.10 Radial velocity of the nozzle outlet

4 結論

(1) 本文采用優化拉丁超立方實驗設計方法生成初始樣本點，建立Kriging模型后，基于非精確搜索改進的增強精英保留的多種群協同遺傳算法進行序列采樣更新代理模型。結果表明，該代理模型優化框架能夠以較少的計算量實現喉栓式噴管型面的一體化優化設計，對工程實踐具有一定的指導意義。

(2) 優化后的喉栓型面頭部鈍度減弱，三種典型工況下喉栓下游的激波明顯減弱，亞音速區域減小，符合實際物理過程。

(3) 優化后的喉栓式噴管出口截面徑向速度明顯減小，燃氣的流動損失有所減小，發動機平均實際比沖增大，其綜合性能得到改善。