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

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

0 引言

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

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

1 優(yōu)化設(shè)計模型

1.1 型面參數(shù)化

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

(1)

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

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

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

1.2 數(shù)值計算模型

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

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

(2)

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

其中，散度的表達式為

(3)

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

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

(4)

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

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

由于噴管及喉栓型面變化較為復(fù)雜，本文采用非結(jié)構(gòu)網(wǎng)格進行網(wǎng)格劃分，并對壁面處進行適當(dāng)加密處理，如圖2所示。網(wǎng)格質(zhì)量檢測顯示，所劃分網(wǎng)格正交性良好，滿足計算要求。網(wǎng)格尺度劃分驗證結(jié)果如表1所示?？梢钥闯觯?dāng)網(wǎng)格數(shù)量大于60 000時，流場結(jié)構(gòu)變化較小。因此，選取60 000網(wǎng)格進行計算。

圖2 計算域網(wǎng)格劃分Fig.2 The grid of calculation domain

表1 網(wǎng)格無關(guān)性驗證

2 基于代理模型的優(yōu)化方法

基于代理模型的優(yōu)化方法在20世紀(jì)80年代被首次提出，并成功應(yīng)用于計算耗時的復(fù)雜模型優(yōu)化問題。該方法基于少量的高精度采樣點建立近似模型，用以在優(yōu)化過程中替代復(fù)雜的高精度學(xué)科模型，通過對近似模型的搜索確定新的采樣點并不斷重構(gòu)更新近似模型，以期逐漸逼近真實學(xué)科模型，預(yù)測全局最優(yōu)解?；诖砟Ｐ偷膬?yōu)化可以大幅減少真實模型的調(diào)用次數(shù)，從而提高優(yōu)化設(shè)計效率，降低設(shè)計成本。

其主要步驟包括實驗設(shè)計采樣、構(gòu)造代理模型以及序列近似采樣等?；诖砟Ｐ偷暮硭ㄊ阶兺屏Πl(fā)動機型面優(yōu)化流程如圖3所示。

2.1 試驗設(shè)計

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

2.2 Kriging模型

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

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

(5)

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

[(),()]=(,)

(6)

基于無偏性以及預(yù)測方差最小準(zhǔn)則，可求得模型預(yù)測函數(shù)為

(7)

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

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

(8)

圖3 型面優(yōu)化流程圖Fig.3 Optimization design flow chart for the contour

2.3 優(yōu)化算法及序列采樣準(zhǔn)則

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

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

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

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

(9)

其停止準(zhǔn)則為

(10)

非精確搜索能夠逐步地引導(dǎo)搜索過程朝著好的方向發(fā)展，可以有效陷入局部最優(yōu)。因此，本文基于非精確搜索的開發(fā)準(zhǔn)則序列填充新的樣本點。

3 優(yōu)化算例及結(jié)果分析

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

(11)

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

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

表2 設(shè)計變量取值范圍及基準(zhǔn)構(gòu)型

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

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

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

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

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

表3 優(yōu)化結(jié)果對比

(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為噴管出口截面的徑向速度分布示意圖。某種程度上，發(fā)動機的比沖損失與徑向速度成正相關(guān)。由圖10可知，優(yōu)化后型面的徑向速度在三種工況下都明顯減小。因此，優(yōu)化后的型面能夠減小流動損失，獲得更高的平均實際比沖。

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

4 結(jié)論

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

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

(3) 優(yōu)化后的喉栓式噴管出口截面徑向速度明顯減小，燃?xì)獾牧鲃訐p失有所減小，發(fā)動機平均實際比沖增大，其綜合性能得到改善。