基于響應(yīng)面法的喉栓式噴管型面優(yōu)化設(shè)計(jì)①

成 沉，鮑福廷，劉 旸，許 昊

(西北工業(yè)大學(xué)航天學(xué)院，西安 710072)

0 引言

為提高導(dǎo)彈的機(jī)動(dòng)性與靈活性，使得系統(tǒng)滿足多任務(wù)需求，作為動(dòng)力裝置的固體火箭發(fā)動(dòng)機(jī)要具備更強(qiáng)的推力調(diào)節(jié)能力。在噴管中嵌入一根喉栓，就能在不改變傳統(tǒng)固體火箭發(fā)動(dòng)機(jī)結(jié)構(gòu)形式情況下，實(shí)現(xiàn)推力大小的調(diào)節(jié)［1－3］。這種推力調(diào)節(jié)技術(shù)因其推力調(diào)節(jié)比大，理論上能實(shí)現(xiàn)無(wú)級(jí)調(diào)節(jié)而受到越來(lái)越多的關(guān)注。

噴管優(yōu)化設(shè)計(jì)是喉栓式固體火箭發(fā)動(dòng)機(jī)設(shè)計(jì)中的一個(gè)重要方面。由于在噴管中嵌入的喉栓會(huì)擾亂傳統(tǒng)噴管中的氣流路徑，使噴管擴(kuò)張段產(chǎn)生復(fù)雜的激波系，甚至引發(fā)流動(dòng)分離等現(xiàn)象，造成比沖損失。所以，喉栓式噴管的優(yōu)化設(shè)計(jì)需綜合考慮喉栓的型面與噴管型面的相互作用，減小因喉栓插入而引起的比沖損失。

本文通過(guò)CFD模擬對(duì)喉栓式噴管的噴管型面及喉栓型面進(jìn)行綜合氣動(dòng)優(yōu)化設(shè)計(jì)，探尋兩者的設(shè)計(jì)參數(shù)對(duì)比沖損失的影響，為喉栓式噴管的設(shè)計(jì)提供參考。由于CFD本身計(jì)算量較大，若直接進(jìn)行優(yōu)化必須進(jìn)行大范圍的搜索，需消耗巨大的計(jì)算量。本文采用響應(yīng)面法，在設(shè)計(jì)空間內(nèi)對(duì)有限的設(shè)計(jì)試驗(yàn)點(diǎn)進(jìn)行CFD計(jì)算，并在此基礎(chǔ)上建立響應(yīng)面模型，再結(jié)合序列二次規(guī)劃法(NLPQL)，對(duì)噴管型面及喉栓型面進(jìn)行綜合氣動(dòng)優(yōu)化設(shè)計(jì)。結(jié)果表明，影響比沖效率的因素主要為噴管的設(shè)計(jì)參數(shù)，在進(jìn)行喉栓式噴管的優(yōu)化設(shè)計(jì)時(shí)，噴管部分完全可采用傳統(tǒng)噴管的優(yōu)化設(shè)計(jì)方法，再針對(duì)喉栓進(jìn)行匹配性優(yōu)化。應(yīng)用響應(yīng)面法，能以相對(duì)較小的計(jì)算量進(jìn)行快速優(yōu)化設(shè)計(jì)，滿足工程優(yōu)化需求，提高工程優(yōu)化效率。

1 CFD模型

1．1 數(shù)學(xué)模型

采用可壓縮流動(dòng)的二維軸對(duì)稱守恒型N-S方程。對(duì)流項(xiàng)采用二階迎風(fēng)格式進(jìn)行離散，粘性項(xiàng)按照中心差分格式進(jìn)行離散。湍流模型采用RNG k-ε模型，增強(qiáng)壁面函數(shù)。

1．2 邊界條件

邊界處理是流場(chǎng)數(shù)值模擬的重要環(huán)節(jié)之一，當(dāng)喉栓從完全打開(kāi)位置向著發(fā)動(dòng)機(jī)喉部運(yùn)動(dòng)時(shí)，喉部面積減小，燃燒室壓力升高，發(fā)動(dòng)機(jī)入口質(zhì)量流率隨燃燒室壓力升高而增加。因此，入口邊界條件定義為自適應(yīng)的質(zhì)量入口，在計(jì)算過(guò)程中讀入發(fā)動(dòng)機(jī)燃燒室壓力，自動(dòng)調(diào)整入口質(zhì)量流率，使其滿足藥柱的燃速變化規(guī)律［3］。單位面積質(zhì)量流率為

式中 Ab為藥柱截面積;Ain為噴管入口面積;ρp為推進(jìn)劑密度;pc為燃燒室壓強(qiáng);n為燃速指數(shù);a為燃速系數(shù)。

發(fā)動(dòng)機(jī)噴管出口面上為超音速出口，不需給定任何邊界條件，全部氣流參數(shù)二階外推即可。

2 響應(yīng)面法

響應(yīng)面法是試驗(yàn)設(shè)計(jì)與數(shù)理統(tǒng)計(jì)相結(jié)合的優(yōu)化方法，當(dāng)實(shí)驗(yàn)結(jié)果與已知參數(shù)間的函數(shù)關(guān)系為隱式時(shí)，在試驗(yàn)測(cè)量、經(jīng)驗(yàn)公式或數(shù)值分析的基礎(chǔ)上，對(duì)指定的設(shè)計(jì)點(diǎn)集合進(jìn)行連續(xù)的試驗(yàn)，可設(shè)計(jì)實(shí)驗(yàn)結(jié)果與參變量間的函數(shù)關(guān)系，并求得參變量的系數(shù)，最終建立響應(yīng)與參變量間的函數(shù)關(guān)系，在設(shè)計(jì)空間構(gòu)造測(cè)定量的全局逼近。響應(yīng)面方法計(jì)算簡(jiǎn)單，因而可僅采用簡(jiǎn)單的代數(shù)形式，依靠目標(biāo)函數(shù)本身的性質(zhì)確定最優(yōu)解。通過(guò)回歸模型的選擇，可擬合復(fù)雜的響應(yīng)關(guān)系，具有良好的魯棒性［5］。本文用設(shè)計(jì)試驗(yàn)點(diǎn)的CFD數(shù)值計(jì)算結(jié)果替代實(shí)驗(yàn)結(jié)果生成響應(yīng)面，可大大減少尋求最優(yōu)解的優(yōu)化計(jì)算次數(shù)。

2．1 響應(yīng)面初始樣本點(diǎn)選取方法

幾乎所有的基于響應(yīng)面的優(yōu)化算法，第一步就是構(gòu)造初始響應(yīng)面，而初始響應(yīng)面構(gòu)造的首要問(wèn)題就是初始樣本點(diǎn)的選取，科學(xué)合理的樣本點(diǎn)在樣本空間的分布應(yīng)是均勻的。本文采用中心復(fù)合實(shí)驗(yàn)設(shè)計(jì)(CCD)，生成初始試驗(yàn)樣本點(diǎn)。CCD的主要特點(diǎn)是可評(píng)估因素的非線性影響。

2．2 響應(yīng)面回歸模型

選擇適合的響應(yīng)面回歸模型是至關(guān)重要的。這關(guān)系到響應(yīng)面模型是否能真實(shí)反映目標(biāo)函數(shù)與設(shè)計(jì)變量之間的關(guān)系。

常用的響應(yīng)面模型有多項(xiàng)式(polynomial regression surrogate，PRS)、徑向基函數(shù) (radial basis functions，RBF)、支持向量回歸(support vector regression，SVR)、多元自適應(yīng)回歸樣條(multivariate adaptive regression spline，MARS)以及 Kriging等。本文選擇多項(xiàng)式模型及Kriging模型分別生成響應(yīng)面，并進(jìn)行了對(duì)比分析。

多項(xiàng)式模型:通常選擇一階多項(xiàng)式為響應(yīng)面模型，但一階模型很難反映真實(shí)的響應(yīng)情況;而選擇大于二階的多項(xiàng)式雖然有較高的擬合精度，但它由于包含較多的項(xiàng)，需付出較大的計(jì)算代價(jià)，尤其是在多變量情況下，擬合響應(yīng)面要花費(fèi)的計(jì)算時(shí)間將是無(wú)法承受的。對(duì)于氣動(dòng)優(yōu)化問(wèn)題，相對(duì)來(lái)說(shuō)二階模型形式比較靈活，對(duì)真實(shí)響應(yīng)近似程度較好。本文基于CFD的計(jì)算結(jié)果，采用完全二次多項(xiàng)式構(gòu)造響應(yīng)面模型。其一般公式為

式中 k為變量個(gè)數(shù);ε為統(tǒng)計(jì)誤差。

Kriging模型:Kriging方法是一種通過(guò)已知點(diǎn)來(lái)預(yù)測(cè)未知觀察點(diǎn)的一種插值方法。Kriging方法利用方差的變化來(lái)表達(dá)空間的變化，且可保證由空間分布得到的預(yù)測(cè)值的誤差最小。

由于Kriging方法生成的響應(yīng)面經(jīng)過(guò)所有的試驗(yàn)樣本點(diǎn)，無(wú)法直接通過(guò)試驗(yàn)樣本點(diǎn)與響應(yīng)面計(jì)算值的擬合程度來(lái)評(píng)估響應(yīng)面模型有效性。文中采取額外隨機(jī)抽取測(cè)試樣本點(diǎn)的方法，用測(cè)試樣本點(diǎn)與響應(yīng)面計(jì)算值的擬合程度來(lái)更好地評(píng)估響應(yīng)面模型有效性。

2．4 優(yōu)化算法

本文采用序列二次規(guī)劃法在響應(yīng)面上進(jìn)行尋優(yōu)。這種算法假設(shè)目標(biāo)函數(shù)是連續(xù)可微的?；舅枷胧菍⒛繕?biāo)函數(shù)以二階拉氏方程展開(kāi)，并把約束條件線性化，使其轉(zhuǎn)化為一個(gè)二次規(guī)劃問(wèn)題。二階方程通過(guò)Quasi-Newton公式得到了改進(jìn)，而且加入了直線搜索，提高了算法的穩(wěn)定性。

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

3．1 優(yōu)化模型

本文優(yōu)化的噴管型面采用雙圓弧法設(shè)計(jì)，而喉栓型面為頭部帶小圓弧的圓錐形，噴管及喉栓構(gòu)型如圖1所示。一般設(shè)計(jì)喉栓式噴管時(shí)，要滿足噴管長(zhǎng)度限制的要求。為了滿足發(fā)動(dòng)機(jī)的內(nèi)彈道要求，還需保證噴管的最大、最小等效喉部截面積(喉栓插入噴管中所形成的最小截面面積)一定。所以，本文在噴管長(zhǎng)度和喉栓直徑一定、且保證等效喉部面積恒定的條件下，為了達(dá)到比沖損失最小目標(biāo)，選擇了圖1所示的設(shè)計(jì)參數(shù)，對(duì)噴管型面及喉栓頭部型面進(jìn)行了優(yōu)化。

圖1 噴管及喉栓型面構(gòu)型圖Fig．1 The contour of the pintle nozzle

其中，噴管設(shè)計(jì)參數(shù)為噴管收斂段與水平方向的夾角α，喉部小圓弧半徑R1，擴(kuò)張段大圓弧半徑R2，噴管出口半徑Re;喉栓設(shè)計(jì)參數(shù)為頭部圓錐半角θ，頭部小圓弧半徑r。

本文以與喉栓式噴管等效喉部截面積相等的最優(yōu)標(biāo)準(zhǔn)噴管(無(wú)喉栓的傳統(tǒng)噴管)的比沖為標(biāo)準(zhǔn)，評(píng)價(jià)喉栓式噴管的比沖損失。比沖損失定義如下:

式中 Isp為喉栓式噴管的比沖;Iss為最優(yōu)標(biāo)準(zhǔn)噴管的比沖。

最優(yōu)標(biāo)準(zhǔn)噴管的比沖是基于響應(yīng)面法對(duì)相同的噴管設(shè)計(jì)參數(shù)，在同樣條件下，對(duì)標(biāo)準(zhǔn)噴管進(jìn)行優(yōu)化后計(jì)算得到的。

3．2 結(jié)果及分析

基于CCD方法，生成了45個(gè)初始試驗(yàn)樣本點(diǎn)(DOE points)，并隨機(jī)抽取了38個(gè)測(cè)試樣本點(diǎn)(Verification points)進(jìn)行CFD計(jì)算。分別用多項(xiàng)式模型及Kriging模型生成響應(yīng)面。圖2和圖3分別顯示了2種模型的擬合效果。圖中橫坐標(biāo)為樣本點(diǎn)CFD計(jì)算值，縱坐標(biāo)為響應(yīng)面預(yù)測(cè)值;方塊表示初始試驗(yàn)樣本點(diǎn)，圓圈表示測(cè)試點(diǎn)，樣本點(diǎn)越接近中線，說(shuō)明擬合程度越好?？煽闯?，對(duì)于測(cè)試樣本點(diǎn)，2種模型的預(yù)測(cè)效果都不錯(cuò)，多項(xiàng)式模型的均方根誤差(Root Mean Square Error)RMSE=0．009%，Kriging 模型的 RMSE=0．005%。均方根誤差RMSE(%)的定義如下:

式中 yi為CFD計(jì)算的準(zhǔn)確值;y^i為響應(yīng)面模型預(yù)測(cè)值;N為評(píng)估測(cè)試點(diǎn)數(shù)。

圖2 多項(xiàng)式模型擬合圖Fig．2 The response surface created by the PRS

圖3 Kriging模型擬合圖Fig．3 The response surface created by the Kriging

可見(jiàn)，2種模型的精度都足以滿足要求，而Kriging模型初始試驗(yàn)樣本點(diǎn)全部通過(guò)響應(yīng)面，多項(xiàng)式模型中存在偏差較大的初始試驗(yàn)樣本點(diǎn)。相比而言，Kriging模型表現(xiàn)較好，較適用于喉栓式噴管的優(yōu)化設(shè)計(jì)問(wèn)題。本文優(yōu)化時(shí)，選擇Kriging模型。用序列二次規(guī)劃法在Kriging響應(yīng)面上優(yōu)化的結(jié)果如表1所示。在響應(yīng)面上尋優(yōu)的結(jié)果，使比沖損失從9．01%降到了4．52%。對(duì)此優(yōu)化點(diǎn)進(jìn)行CFD計(jì)算，得到的比沖損失為4．52%。響應(yīng)面誤差為1．6%，精度滿足要求。對(duì)優(yōu)化點(diǎn)整化之后，再進(jìn)行 CFD計(jì)算，得到最終的比沖損失為4．50%。圖4和圖5是優(yōu)化前后流場(chǎng)的對(duì)比?？煽闯?，優(yōu)化后喉栓頭部的亞音速區(qū)明顯減小，從喉栓頭部延伸出的斜激波強(qiáng)度減弱。從圖6可看出，噴管出口處的徑向速度明顯減小。由于只有軸向速度對(duì)比沖有貢獻(xiàn)，而徑向速度越大，比沖損失越大。優(yōu)化后，喉栓對(duì)流場(chǎng)的擾動(dòng)作用減弱，流動(dòng)損失自然減小了。

表1 優(yōu)化結(jié)果(Kriging模型)Table 1 Optimization results(Kriging model)

圖4 優(yōu)化前流場(chǎng)馬赫數(shù)云圖Fig．4 Mach number contour before the optimization

圖5 優(yōu)化后流場(chǎng)馬赫數(shù)云圖Fig．5 Mach number contour after the optimization

圖6 噴管出口徑向速度Fig．6 The radial velocity of the outlet

3．3 關(guān)鍵參數(shù)分析

各設(shè)計(jì)參數(shù)的敏感性如圖7所示。其中，噴管設(shè)計(jì)參數(shù)擴(kuò)張段大圓弧半徑R2、噴管出口半徑Re及喉栓設(shè)計(jì)參數(shù)頭部圓錐半角θ對(duì)比沖損失的影響較大。

圖7 參數(shù)敏感性Fig．7 The sensitivity of parameters

以優(yōu)化點(diǎn)為基礎(chǔ)，再進(jìn)一步針對(duì)這幾個(gè)參數(shù)進(jìn)行單獨(dú)的優(yōu)化，得到的響應(yīng)面如圖8所示?？煽闯?，R2小于140 mm時(shí)，比沖損失較小;超過(guò)140 mm時(shí)，比沖損失會(huì)急劇增加。比沖損失隨Re的增大幾乎呈線性減小，在允許范圍內(nèi)，Re越大越好。θ對(duì)比沖損失的影響相對(duì)前兩個(gè)參數(shù)較小，小于70°時(shí)，隨θ的增大，比沖損失增大的較明顯;超過(guò)70°時(shí)，比沖損失幾乎不變。

再對(duì)這幾個(gè)參數(shù)進(jìn)行單獨(dú)尋優(yōu)，組成新的優(yōu)化點(diǎn)，比沖損失下降為3．63%。優(yōu)化點(diǎn)為 α=44°，R1=3 mm，R2=120 mm，Re=26 mm ，θ=10°，r=1．3 mm。

3．4 優(yōu)化方式改進(jìn)

總體而言，噴管設(shè)計(jì)參數(shù)比喉栓設(shè)計(jì)參數(shù)的影響要明顯一些，而喉栓的設(shè)計(jì)參數(shù)對(duì)比沖損失的影響不是很大，影響比沖效率的主要因素還是噴管的設(shè)計(jì)參數(shù)。

圖8 關(guān)鍵參數(shù)響應(yīng)面Fig．8 The response surface of the key parameters

基于上述分析，改變優(yōu)化方式，先對(duì)噴管進(jìn)行優(yōu)化，再對(duì)喉栓進(jìn)行優(yōu)化。這樣試驗(yàn)點(diǎn)從45個(gè)，減小到了34個(gè)，減小了計(jì)算量。優(yōu)化結(jié)果與綜合優(yōu)化結(jié)果的對(duì)比如表2所示，分開(kāi)優(yōu)化得到的優(yōu)化結(jié)果與綜合優(yōu)化的差別不大，這種優(yōu)化方式能在達(dá)到同樣優(yōu)化效果的前提下，減小計(jì)算量。

在進(jìn)行喉栓式噴管的優(yōu)化設(shè)計(jì)時(shí)，噴管部分完全可采用傳統(tǒng)噴管的優(yōu)化設(shè)計(jì)方法，再針對(duì)喉栓進(jìn)行匹配性優(yōu)化。

表2 改進(jìn)的優(yōu)化結(jié)果Table 2 Improving optimization results

4 結(jié)論

(1)喉栓式噴管的型面優(yōu)化問(wèn)題需綜合考慮喉栓型面和噴管型面的相互影響，以得到最優(yōu)氣動(dòng)型面。

(2)對(duì)于文中的問(wèn)題，Kriging模型比完全二階多項(xiàng)式模型擬合程度高，預(yù)測(cè)能力好，適用于喉栓式噴管的優(yōu)化設(shè)計(jì)。

(3)基于響應(yīng)面法進(jìn)行喉栓式噴管的型面優(yōu)化設(shè)計(jì)，計(jì)算量相對(duì)較小，計(jì)算結(jié)果的精度和可靠度相對(duì)較高，可快速進(jìn)行氣動(dòng)優(yōu)化設(shè)計(jì)，具有一定的工程實(shí)用價(jià)值。

(4)影響比沖效率的因素主要為噴管的設(shè)計(jì)參數(shù)，在進(jìn)行喉栓式噴管的優(yōu)化設(shè)計(jì)時(shí)，噴管部分完全可采用傳統(tǒng)噴管的優(yōu)化設(shè)計(jì)方法，再針對(duì)喉栓進(jìn)行匹配性優(yōu)化。

［1］ 張淑慧，胡波，孟雅桃．推力可控固體火箭發(fā)動(dòng)機(jī)應(yīng)用及發(fā)展［J］．固體火箭技術(shù)，2002，25(4)．

［2］ Randall Smith-kent，Hai-tien Loh，Pawel Chwalowski．Analytical contouring of pintle nozzle exit cone using computational fluid dynamics［R］．AIAA 1995-2877．

［3］ 李娟，李江．喉栓式固體火箭發(fā)動(dòng)機(jī)噴管性能影響研究［J］．彈箭與制導(dǎo)學(xué)報(bào)，2007，27(3):154-160．

［4］ 鄒林君．基于Kriging模型的全局優(yōu)化方法研究［D］．武漢:華中科技大學(xué)，2011．

［5］ 熊俊濤．基于響應(yīng)面方法的氣動(dòng)優(yōu)化設(shè)計(jì)［D］．西安:西北工業(yè)大學(xué)，2005．

［6］ 王筱蓉，周長(zhǎng)省，鞠玉濤，等．固體火箭發(fā)動(dòng)機(jī)特型噴管的型面設(shè)計(jì)［J］．彈道學(xué)報(bào)，2008，20(4):77-80．