面向航程能力的固體火箭發動機方案設計優化

宋少倩，陳永信，任鵬飛，周文勇，李偉喆

(1.中國運載火箭技術研究院,北京 100076；2.上海機電工程研究所,上海 201109；3.北京宇航系統工程研究所,北京 100076)

為適應未來高威脅作戰環境的迫切需要，對飛行器的航程提出了越來越高的要求[1]。由于發射平臺的約束，飛行器的質量與尺寸往往有著嚴苛的限制，要實現高性能的戰術技術指標，除了先進總體布局及新理論、新材料、新技術的應用外，必須有高性能的推進系統與之相匹配[2]。

固體火箭發動機(Solid rocket motor,SRM)具有結構簡單、快速反應能力強、工作可靠性高等優點，廣泛應用于運載火箭、導彈等飛行器系統中[3]。傳統上在飛行器設計初期，發動機設計單位按給定質量與尺寸要求給出初步的性能參數，供總體設計評估，在不滿足航程要求時再迭代設計，由于傳統的設計方法沒有直接引入航程能力等總體性能的影響[4]，故難以得到系統最優需求的發動機設計方案。因此在飛行器設計初期，除了不斷提升SRM自身性能外，基于總體性能驅動的一體化設計優化是提高飛行器航程的有效途徑之一。后者需要結合飛行器總體布局合理設計發動機的內彈道特性，并通過飛行器軌跡合理分配與使用發動機的能量。

早期的一體化設計優化方法受限于計算資源和時間成本，往往采用大量近似解析公式、經驗數據，適用范圍受限，對于新型高性能飛行器，其設計結果往往與后續詳細設計的性能差別較大，容易導致設計的反復[4-5]。為保證設計優化結果的可信度，需要在分析模型中采用大量的高精度數值計算，這導致進行一次方案評估的計算時間成本很高。若直接采用傳統參數優化方法調用精確分析模型進行設計優化，則優化總時間難以承受。因此，有必要采用近似優化(Approximate optimization,AO)方法[6-11]進行求解，即采用計算代價較小的代理模型來代替復雜且耗時的數值分析，從而以較小的時間成本獲得較優的設計方案。

AO方法的策略主要分為：直接法(One-shot approach,OA)[6-8]和序列迭代法(Updating approach, UA)[7-11]。OA策略僅構建一次代理模型，并通過一定準則校驗代理模型的全局精度，而后在該代理模型上進行優化。UA策略則采用較少的初始樣本點構建低精度代理模型，并根據一定加點準則，在序列優化過程中不斷更新代理模型的信息，直至滿足優化停止條件。相比于OA策略，UA策略無需校驗代理模型的全局精度，特別是對于高維優化問題，UA策略調用耗時模型的次數更少。近年來，基于UA策略的序列近似優化(Sequential approximate optimization,SAO)方法得到了廣泛發展[7-11]。例如，基于Kriging代理模型的高效全局優化(Efficient global optimization,EGO)算法[8]，其定義了目標函數值改善期望(Expected improvement,EI)，并選擇EI最大值點作為新的樣本點，該算法具有較好的全局搜索性，但在最優解附近收斂速度較慢[9]；基于目標函數預測值最小(Minimizing the predicted objective function,MP)準則的優化算法，其屬于貪婪算法，收斂速度較快，但容易陷入局部最優解；結合多種加點準則的并行加點算法[7,9,11]，其能夠彌補單一加點準則的不足，目前仍是研究的熱點。

本文以固體動力助推滑翔式飛行器為背景，針對其規模約束嚴格、飛行環境嚴酷及高性能需求的特點，提出面向航程能力的SRM方案設計優化方法，建立復雜三維裝藥發動機性能分析的參數化模型，以及基于自適應偽譜法的多約束飛行器航程能力評估模型，提出基于Kriging代理模型的改進混合整數序列優化算法，開展發動機與飛行器軌跡一體化設計優化。

1 模型建立

1.1 參數優化問題

不失一般性，參數優化問題可表示為

minf(z)

(1)

式中：z為設計變量，f(z)為目標函數，gi(z)為第i個不等式約束，ng為不等式約束數量，hj(z)為第j個等式約束，nh為等式約束數量。

對于SRM總體方案設計優化問題，本文以單級助推滑翔式飛行器為例，在不改變飛行器氣動外形的條件下，通過設計優化發動機的裝藥參數和噴管參數，即設計變量z，來提升飛行器的航程能力，即目標函數f(z)，其中要求發動機總質量不大于要求值，即不等式約束gi(z)，發動機的外徑、總長度和裝藥質量為固定值，即等式約束hj(z)。

(2)

式中：s為飛行器的航程；dM、dMd分別為發動機外徑和要求值；lM、lMd分別為發動機總長度和要求值；mp、mpd分別為裝藥質量和要求值；mM為發動機總質量；下標max表示最大值或上限值；|z表示設計方案z所對應的結果。

值得注意的是，設計變量z中除了實數型變量外，往往還包括整數型變量，例如翼柱型裝藥的翼槽數量等。因此，式(1)所示的參數優化問題為混合整數優化問題。

1.2 發動機性能分析模型

1.2.1 裝藥模型

考慮到翼柱形裝藥結構完整性好，體積裝填分數高，燃面可調范圍大，目前廣泛應用于各類SRM中[12-13]。本文發動機裝藥采用前、后三角翼柱形構型，如圖1所示。

1.2.2 殼體模型

發動機前、后封頭均采用橢球比相同的橢球封頭，殼體圓筒段的厚度δc1采用最大應變能強度理論確定，橢球封頭厚度沿橢球表面從邊緣向中心逐漸加厚，頂點處的厚度δc2采用最大應力強度理論確定，其計算式[13]為

(3)

式中：pc為燃燒室壓強；kc、φc、ξc分別為安全系數、壓力波動系數和焊縫強度系數，本文分別取為1.25、1.10和0.90；[σc]為殼體材料許用應力；d1為圓筒的內徑；km為前、后封頭的橢球比，取為2。

圖1 翼柱形裝藥幾何示意

1.2.3 噴管模型

噴管采用特型噴管，入口段采用橢圓形型面，擴張段采用優化3次多項式型面，入口段與擴張段之間采用圓弧過渡。其中擴張段型面的表達式為[13]

(4)

式中：xn為擴張段型面上的點到噴管入口段前端面的距離；yn為擴張段型面上的點到旋成軸線的距離；lnb、lnt、lne分別為噴管收斂段、過渡段和擴張段的長度；rnt、rne分別為喉部半徑和擴張段出口半徑；rnm為過渡段圓弧半徑，本文取rnm=rnt；αnm、αne分別為擴張段入口和出口的擴張半角，本文分別取值為26°和22°。

圖2 噴管構型

1.2.4 質量模型

發動機總質量mM可以表示為

(5)

式中：ρp、Vp分別為裝藥的密度和體積；ρc、Vc分別為燃燒室殼體的材料密度和體積；ρq、Vq分別為燃燒室絕熱層的材料密度和體積；ρnb、Vnb分別為喉襯的材料密度和體積；ρne、Vne分別為擴張段絕熱層的材料密度和體積；ma為其他裝置質量，包括：點火裝置、發火裝置、連接裝置等，按經驗估算。其中Vnb和Vne采用文獻[14]方法計算，其余部件體積由幾何模型計算獲得。不同時刻發動機相對飛行器頭部頂點處的質心位置xGM亦由幾何模型計算獲得。

1.2.5 推力模型

對于長徑比較小的SRM，采用零維內彈道模型進行性能預示可滿足工程設計的精度要求[12]。零維內彈道模型假設燃燒室內的燃速、壓強、溫度等參數隨空間分布均勻，根據質量守恒關系，則可建立燃燒室內的流動控制方程為

(6)

式中：Vce為燃燒室空腔自由容積；Rg、Tg分別為燃氣的氣體常數和燃燒溫度；ep為燃燒厚度；Ap為裝藥在燃燒厚度為ep時的燃面面積；a0、n0分別為裝藥在標準設計狀態下的燃速系數和燃速壓強指數；C*為裝藥的特征速度；ant為喉部燒蝕速率，取為0.24 mm/s。

則發動機推力大小FT為

(7)

式中：FT0為動推力，γg為燃氣比熱比，pe為噴管出口壓力，pa為環境壓力，εe為噴管膨脹比，其中εe與pe存在關系為

(8)

可采用牛頓迭代法求解式(8)以獲得pe。

對于復雜三維裝藥的燃面計算問題，最小距離函數(Minimum distance function,MDF)法能夠自動處理燃燒過程中裝藥幾何特征消失/增加等拓撲變化，具有較好的穩定性、通用性和計算精度[13,15]。該方法首先建立裝藥體網格到初始燃面面網格的最小距離函數，根據平行層燃燒定律，距離值等于燃燒厚度ep的等值面即為當前燃面Ap，從而為式(6)提供輸入數據。詳細的計算步驟見文獻[15]，本文不再贅述。

1.3 航程能力評估模型

1.3.1 飛行器運動模型

為有效進行航程能力分析，在方案設計優化階段可忽略地球扁率和自轉影響，建立縱向平面有動力無側滑的質點運動方程。同時，在方程中引入迎角變化率，以便考慮控制系統限制的影響。則運動方程可表示為

(9)

式中：v、θ、h、s、α分別為飛行器的速度、航跡角、高度、航程和迎角；ud為迎角變化率控制指令；mT為飛行器有效載荷質量，其值為500 kg；μ、RE分別為地球引力常數和半徑，其值分別為3.986 005×1014m3/s2和6 371 004 m；FD、FL分別為飛行器的阻力和升力，其可表示為

(10)

式中：CA、CN分別為軸向力系數和法向力系數；CD、CL分別為阻力系數和升力系數；SR為飛行器的參考面積；Ma為馬赫數；q為動壓；ρ、a分別為當地大氣密度和當地聲速；δ為俯仰舵偏角，根據“瞬時平衡假設”有

(11)

式中：CMZ、xGT分別為相對飛行器頭部頂點處的俯仰力矩系數和有效載荷質心位置，xGT= 1.8 m；LR為飛行器的參考長度，值為6 m。由于式(10)難以解析求解，故本文采用牛頓迭代法來獲得當前Ma和α所對應的δ。

1.3.2 飛行約束條件

助推滑翔式飛行器在大氣層內高速飛行，其力熱環境嚴酷，因此，在航程評估時需要考慮防熱、載荷、控制的限制約束，即

(12)

式中：RN為駐點曲率半徑，值為0.03；γ為完全氣體比熱比，其值為1.4；γh為高溫氣體比熱比，其值為1.2；VQ為駐點熱流密度；ny為法向過載；g0為海平面引力加速度，其值為9.806 6 m/s2；下標min為最小值或下限值。

此外，飛行末端還需要進行能量管理，即

(13)

式中：tf、vf、θf分別為末段的時刻、速度和航跡角。

1.3.3 基于最優控制方法的航程能力評估

當給定一組設計方案z時，根據發動機性能分析模型方法可獲得所對應的發動機推力FT、發動機質量mM、發動機質心位置xGM隨時間t的變化關系。至此，對于式(2)中目標函數的計算，則可轉化為對連續時間最優控制問題的求解，即

f(z)=smaxz=-s(tf)

(14)

式中：狀態變量x=[vθhsα]T；控制變量u=ud；C為路徑約束，見式(12)；φ為端點約束，見式(13)。

上述航程能力評估的實質為求解最優控制問題，即為內層優化。為了保證式(1)所示的參數優化問題能夠穩健有效進行，在保證式(14)的求解精度前提下，需要對求解方法的快速性和穩定性提出要求。自適應偽譜法將連續時間最優控制問題轉化為非線性規劃問題，根據一定準則來保證問題轉化的精度，這類方法對于初值不敏感，近年來廣泛應用于飛行器軌跡優化中[16-19]。

本文采用改進自適應Legendre-Gauss-Radau (LGR)偽譜法[19]作為內層優化求解器。文獻[19]驗證了該方法對于多約束飛行器軌跡優化問題的求解快速性和可靠性，該方法根據當前解的相對誤差來調整配點區間數量，在精度不滿足的區間內根據Legendre近似理論增加配點或進行分割，在精度滿足的區間內根據多項式系數分析減少配點或進行合并[18-19]，直至計算精度滿足要求。

1.4 發動機方案評估模型

綜上所述，進行單次發動機總體方案評估的流程如圖3所示。

圖3 單次SRM設計評估流程

2 改進混合整數序列近似優化算法

2.1 Kriging代理模型

Kriging代理模型能夠提供任意位置處函數的預測值以及該預測值的統計學誤差估計。該誤差估計信息可用于指導新樣本點的選擇。

(15)

(16)

本文中回歸基函數采用一次多項式模型[9,20]，相關基函數采用3次樣條模型[9,20]，即

(17)

式中：θ(k)為模型控制參數，可根據極大似然原則[20]，求解參數優化問題式(18)獲得，即

(18)

2.2 混合整數LHD初始采樣

提高投影均勻性和空間均布性，能夠有效豐富設計空間的信息獲取[7]。拉丁超立方設計(Latin hypercube design,LHD)采樣自然滿足投影均勻性，已廣泛應用于各類近似優化算法的初始采樣中。LHD采樣點可表示為

(19)

由于LHD采樣具有隨機性，難以保證空間均布性，故本文將πi和ρi作為優化參數以改善空間均布性，其中優化準則采用Maximin準則[7]，即

(20)

對于整型變量，可將其視作實數型變量處理，并通過四舍五入取整以確保變量始終為整數[21]。

2.3 局部增強的EI加點準則

(21)

則I(z)的數學期望E[I(z)]為

(22)

式中：Φ[·]為標準正態累積分布函數，φ[·]為標準正態概率密度函數。

EI加點準則即選擇E[I(z)]值(稱為EI值)最大的樣本點。對于混合整數參數優化問題，可在滿足等式約束、不等式約束以及整數約束[21]的可行集內搜索EI值最大點作為新的采樣點，即

min(-E[I(z)])

(23)

本文針對上述問題提出一種局部增強的EI加點準則(簡稱為LEEI加點準則)。其基本思想是：在當前最優解附近限定一個小的局部搜索區域，在該區域內采用EI加點準則，以提升局部搜索性能；同時在全局搜索區域內采用EI加點準則，以保留原算法的全局收斂性；引入比例參數對局部搜索范圍進行動態調整，以控制算法收斂。

則引入局部搜索范圍的EI加點優化問題(23)可改寫為

min(-E[I(z)])

(24)

式中：round(·)為四舍五入取整算子，rD為搜索空間縮放因子，z*(n)為當前最優解z*為的第n維值。

2.4 算法流程

基于LEEI加點準則的改進混合整數SAO算法的具體計算步驟為：

1)確定初始采樣點數量(可取ns=10nz)，求解參數優化問題式(20)獲得初始樣本點集{z1,z2,…,zns}，計算各樣本點真實的目標函數值和約束函數值；初始化迭代次數k=0，并設置最大迭代次數kmax(可取kmax=200)和rD的初始值(可取rD=0.2)。

(25)

3)令k=k+1，基于當前樣本點集構建目標函數和約束函數的Kriging代理模型，分別求解參數優化問題式(23)和式(24)獲得全局EI點zG和局部EI點zL。

4)若zG在局部搜索區域內(即zG滿足式(24)中的約束3、4)，則將zL加入樣本點集中,計算zL真實的目標函數值和約束函數值；否則將zG和zL同時加入樣本點集中，分別計算zG和zL真實的目標函數值和約束函數值。

6)判斷收斂條件式(26)，其中εi(i=1,2,3,4)為收斂精度(可根據目標函數和約束函數的值域對應選擇εi的大小)。若滿足則停止；否則返回步驟3)。

(26)

3 仿真驗證

本文燃燒室殼體材料采用40SiMnCrMoV，絕熱層材料采用EPDM，推進劑采用HTPB/AP/Al，噴管喉襯和擴張段材料采用碳/碳復合材料。共優化17個設計變量，包括前、后翼槽數量共2個整型設計變量，以及裝藥參數、噴管參數、燃速系數共15個實數型設計變量，為保證幾何參數化模型的拓撲正確，定義設計變量與幾何參數的關系為

(27)

設計變量搜索范圍見表1。設計約束、飛行條件定義見表2，其中上標P和H分別為發動機工作和不工作時的參數，其余均為全程參數。

表1 設計變量搜索范圍

表2 設計約束、飛行條件定義

本文基于PYTHON腳本語言驅動ABAQUS 6.14實現發動機殼體、絕熱層、裝藥和噴管的幾何參數化建模和四面體/三角形網格自動劃分；采用文獻[15]方法實現燃面計算；采用文獻[19]方法實現改進自適應LGR偽譜法；采用經典的DACE工具箱實現Kriging代理模型的構建；基于MATLAB R2019b平臺編寫混合整數LHD初始采樣程序，以及基于LEEI加點準則的改進混合整數SAO算法(簡稱為LEEI算法)程序，其中混合整數參數優化問題均采用MATLAB自帶的混合整數遺傳算法(Genetic algorithm, GA)進行求解，LEEI算法中的收斂精度分別設置為：ε1=ε2=10-3，ε3=0，ε4=10-1)；仿真計算環境為Windows 7 2.10 GHz，32.0 GB內存，優化結果如圖4所示。

紅色為發動機工作段，黑色為發動機不工作段，實線為參數曲線，虛線為約束邊界

為說明設計優化的有效性，以初始采樣中的最優方案作為優化前的方案(由算法流程步驟2)獲得)進行對比，其中優化前的最大航程為554.90 km，裝藥前、后翼槽數nw1、nw2分別為8和6，裝藥質量mp為593.82 kg，不滿足等式約束要求；采用LEEI算法優化后，最大航程為626.42 km(如圖4(a)所示)，相比優化前提升了12.89%，裝藥的前、后翼槽數nw1、nw2均為8個，發動機外徑dM、總長度lM和裝藥質量mp分別為0.44，2.90，589.97 kg，發動機總質量mM為737.40 kg，其中等式約束函數違反程度為0.03(≤ε4=10-1)，不等式約束函數違反程度為0(≤ε3=0)，均滿足式(2)中的等式和不等式約束要求，驗證了LEEI算法對于等式和不等式約束處理的有效性。

由圖4(b)～(e)可知，飛行器的迎角、迎角變化率、駐點熱流密度、動壓、法向過載和配平舵偏角均全程滿足相應的約束要求。其中迎角在初始爬升時達到了約束上界，即最大爬升迎角是航程的限制因素之一；法向過載在發動機工作時達到了約束上界，而駐點熱密度全程均未達到約束邊界，可知對于本文飛行器，法向過載為主要限制因素，其制約了可用迎角和發動機推力大??；最高點動壓達到了約束下限，由于控制系統對動壓的限制，導致最大飛行高度受到制約，進而影響了航程；配平舵偏角全程變化平穩，僅在滑翔末段達到了約束下界，即限制了可用迎角。由圖4(d)可知，為提高航程，飛行末端的飛行速度達到了約束下界，而航跡角達到了上界，兩者均滿足相應約束要求。由此可知，優化結果滿足全部飛行約束要求，驗證了LEEI算法結果的正確性。

下面為對比分析優化結果，設置2個整型變量為若干固定值，并分別對剩余的15個實數型變量進行優化，優化對比結果如圖5(a)～(d)所示。

圖5 不同翼槽數對應的優化結果

由圖5(a)可知，優化方案的最大航程優于其他情況；由圖5(b)～(d)可知，前8后8翼槽方案的推力接近于雙推力的形式，即采用短時間大推力輔助飛行器爬升，而后采用長時間小推力持續加速，由于發動機長度、外徑以及裝藥質量被嚴格限制，因此當翼槽數較少時，則需要翼面和翼厚更大，以滿足尺寸質量的約束條件，從而導致發動機的初始推力相對更大，工作時間相對更短；對于翼槽數較多時，結果相反；此外，可以發現優化方案對應的關機點高度和速度也相對更大，從而其最大航程也更大。

為進一步驗證LEEI算法的性能，分別采用基于目標函數預測值最小加點準則的SAO算法(簡稱為MP算法)、基于傳統EI加點準則的SAO算法(簡稱為EI算法)和GA算法對本文問題進行求解與對比，其中MP算法來自文獻[14]，EI算法來自MATLAB軟件包Surrogate Model Toolbox。為剝離初始采樣對SAO算法的性能影響，LEEI、MP、EI算法均采用本文混合整數LHD初始采樣算法進行初始采樣，采樣數量為170個，算法性能對比結果如圖6、7和表3所示。

圖6 目標函數的收斂過程

圖7 約束函數的收斂過程

表3 算法性能對比

由圖6、7可知，隨著序列加點與迭代優化，3種SAO算法均得到收斂。MP算法收斂速度最快，滿足全部約束要求，但收斂時的目標函數值較大，即優化的最大航程較小。LEEI和EI算法的結果與GA算法結果一致，均滿足全部約束要求，且優于MP算法的結果，其中EI算法初始收斂速度較快，但后期收斂速度很慢，耗時函數調用次數約為MP算法的2倍；GA算法雖然優化結果較好，但耗時函數調用次數過大，約為SAO算法的30～60倍；本文LEEI算法的耗時函數調用次數多于MP算法，但相比于EI算法顯著減少，由于增加了局部搜索控制策略，從而有效改善了收斂速度，同時保留了原算法的全局搜索性，驗證了LEEI算法的有效性。

4 結 論

1)LEEI算法能夠有效求解SRM方案設計優化的混合整數參數優化問題，優化后飛行器的最大航程提升了12.89%，且設計結果滿足全部約束條件。

2)LEEI算法的優化結果與GA算法的優化結果相一致，但耗時函數的調用次數僅為后者的2%。

3)LEEI算法具有一定全局搜索能力，且相比于EI算法有效提高了局部收斂速度。