基于Kriging 代理模型的大直徑運載火箭蒙皮桁條結(jié)構(gòu)分步優(yōu)化方法研究

吳 棟，王志祥，劉觀日，張大鵬*，雷勇軍

(1. 國防科技大學(xué)空天科學(xué)學(xué)院，長沙410073； 2. 北京宇航系統(tǒng)工程研究所，北京100076)

1 引言

薄壁加筋圓柱殼結(jié)構(gòu)是航空航天結(jié)構(gòu)中常見的結(jié)構(gòu)形式，蒙皮桁條結(jié)構(gòu)是其中一種典型的結(jié)構(gòu)。 現(xiàn)有運載火箭型號中約有80%的箭體艙段采用了整體加筋和桁條加強的薄壁殼體結(jié)構(gòu)，如運載火箭的級間段蒙皮桁條和貯箱等結(jié)構(gòu)[1-2]。針對蒙皮桁條結(jié)構(gòu)開展優(yōu)化設(shè)計能減輕運載火箭箭體的結(jié)構(gòu)質(zhì)量，并且對提升運載火箭的可靠性、安全性和發(fā)射成功率有著重要意義。

近年來，研究人員基于層級設(shè)計理念[1-3]開展了大量多級加筋結(jié)構(gòu)承載能力和優(yōu)化研究。Wang 等[4]提出了一種多級加筋板結(jié)構(gòu)，并對該結(jié)構(gòu)開展了軸壓作用下極限承載能力和缺陷敏感性的研究；Wang 等[5]和郝鵬[6]提出了考慮初始幾何缺陷的薄壁加筋圓柱殼結(jié)構(gòu)在軸壓作用下極限載荷的預(yù)測方法，并基于該方法提出了考慮幾何缺陷的薄壁加筋圓柱殼結(jié)構(gòu)的優(yōu)化設(shè)計方法；范書群等[7]基于有限元方法，提出一種半硬殼結(jié)構(gòu)框-桁匹配設(shè)計方法；張鐵亮等[8]基于代理模型對加筋板的布局進行優(yōu)化；張柱國等[9]基于進化Kriging 模型的優(yōu)化方法對金屬加筋板的尺寸和布局進行設(shè)計；郝鵬等[10]基于徑向基函數(shù)代理模型，提出一種蒙皮桁條結(jié)構(gòu)新型減輕孔設(shè)計。

10 m 級重型運載火箭對結(jié)構(gòu)減重有較大的需求，同時尺度效應(yīng)又增加了設(shè)計難度。 目前，針對10 m 級重型運載火箭薄壁加筋圓柱殼結(jié)構(gòu)的輕質(zhì)優(yōu)化研究較少。 本文以10 m 級重型運載火箭的結(jié)構(gòu)輕質(zhì)優(yōu)化設(shè)計為背景，以大直徑運載火箭級間段蒙皮桁條結(jié)構(gòu)為研究對象，基于Python語言建立蒙皮桁條結(jié)構(gòu)參數(shù)化模型，采用顯式非線性動力學(xué)算法分析結(jié)構(gòu)的后屈曲狀態(tài)，分析不同蒙皮網(wǎng)格規(guī)模對計算精度和計算效率的影響規(guī)律；針對蒙皮桁條結(jié)構(gòu)優(yōu)化中涉及連續(xù)的尺寸優(yōu)化以及離散的拓撲優(yōu)化問題，對優(yōu)化變量進行合理分組，提出基于Kriging 代理模型和組合優(yōu)化算法對蒙皮桁條結(jié)構(gòu)進行分步優(yōu)化，以獲得不同截面形式桁條下的最優(yōu)結(jié)構(gòu)形式。

2 蒙皮桁條結(jié)構(gòu)后屈曲分析

2.1 蒙皮桁條結(jié)構(gòu)有限元建模

蒙皮桁條結(jié)構(gòu)由蒙皮、桁條、中間框和端框組成。 蒙皮內(nèi)側(cè)沿著軸向均勻分布4 個“Π”型截面的中間框，2 個“L”型截面的端框。 蒙皮外側(cè)沿著環(huán)向均勻分布若干根桁條，蒙皮桁條結(jié)構(gòu)的幾何模型如圖1 所示，端框、中間框截面如圖2 所示，桁條截面如圖3 所示。

圖1 蒙皮桁條結(jié)構(gòu)幾何示意圖Fig.1 Geometric diagram of skinned truss structure

圖2 端框和中間框截面示意圖Fig.2 Schematic diagram of the end frame and the middle frame

圖3 桁條截面示意圖Fig.3 Schematic diagram of truss

基于Python 語言建立蒙皮桁條結(jié)構(gòu)的參數(shù)化模型。 由于蒙皮桁條結(jié)構(gòu)使用較薄的蒙皮和桁條，結(jié)構(gòu)上體現(xiàn)為板殼特性，為了能夠準確模擬桁條局部截面平動和轉(zhuǎn)動，采用殼單元劃分網(wǎng)格[11]，桁條腹板高度方向劃分兩層單元。 在蒙皮與桁條、蒙皮與端框、蒙皮與中間框之間施加粘接約束。 為獲得軸壓作用下蒙皮桁條結(jié)構(gòu)的極限載荷，分別在上下端面建立參考點，并與相應(yīng)端面節(jié)點建立多點耦合約束單元，對下參考點施加固支約束，上參考點約束除軸向位移外的其他自由度，勻速施加35 mm 強制位移。 模型選用鋁合金材料，鋁合金材料的彈性模量為70 GPa，泊松比為0.3，密度為2.78×10-6kg/mm3，屈服應(yīng)力為440 MPa，強度極限為550 MPa，延伸率為6%。

2.2 蒙皮桁條結(jié)構(gòu)后屈曲分析

對于軸壓作用下的薄壁加筋圓柱殼結(jié)構(gòu)，一般情況下，允許蒙皮在較小的載荷下發(fā)生局部失穩(wěn)。結(jié)構(gòu)的極限承載能力是由結(jié)構(gòu)的后屈曲狀態(tài)決定的，相較于弧長法和隱式非線性算法，顯式非線性動力學(xué)算法能很好地跟蹤軸壓載荷下薄壁加筋圓柱殼結(jié)構(gòu)的后屈曲行為[12]。 本文采用顯式非線性動力學(xué)算法分析蒙皮桁條結(jié)構(gòu)的后屈曲狀態(tài)。

使用較粗的單元劃分網(wǎng)格會使計算結(jié)果誤差較大，較細的單元劃分網(wǎng)格會導(dǎo)致結(jié)構(gòu)分析計算量加大[13]，需要合理地設(shè)置網(wǎng)格規(guī)模。 為了準確反映蒙皮失穩(wěn)模態(tài)，需分析蒙皮網(wǎng)格規(guī)模對顯式非線性動力學(xué)算法的影響。 采用如表1 中所示的結(jié)構(gòu)初始設(shè)計參數(shù)，得到3 種蒙皮網(wǎng)格規(guī)模對結(jié)構(gòu)極限載荷的影響規(guī)律，如表2 所示。 隨著蒙皮網(wǎng)格規(guī)模的提升，結(jié)構(gòu)的極限載荷減少，表明網(wǎng)格的細化導(dǎo)致計算誤差減小，計算精度提高，同時計算耗時相應(yīng)增加。 筋格內(nèi)網(wǎng)格數(shù)為4×23 與4×31相比，極限載荷計算結(jié)果相差0.56%，但計算耗時減少了1.3 倍。 綜合考慮計算精度和計算效率，在后續(xù)計算中采用筋格內(nèi)網(wǎng)格數(shù)量為4×23的蒙皮網(wǎng)格，如圖4 所示。

表1 蒙皮桁條結(jié)構(gòu)初始設(shè)計參數(shù)【注】Table 1 Initial design parameters of skinned truss structure

3 基于Kriging 代理模型和組合優(yōu)化算法的分步優(yōu)化策略

3.1 優(yōu)化問題描述

蒙皮桁條結(jié)構(gòu)的輕質(zhì)優(yōu)化包含拓撲優(yōu)化、尺寸優(yōu)化和幾何優(yōu)化[14]。 拓撲優(yōu)化是優(yōu)化桁條的數(shù)量；尺寸優(yōu)化是優(yōu)化單元的截面尺寸；幾何優(yōu)化是優(yōu)化各結(jié)構(gòu)的截面形狀[14]，使其在滿足承載能力需求的前提下通過改變截面形狀提升結(jié)構(gòu)的性能。 具體可以描述如式(1)所示：

其中，xHT為桁條參數(shù)， xHTmin為桁條各參數(shù)的下限，xHTmax為桁條各參數(shù)的上限；xZJK為中間框參數(shù)，xZJKmin為中間框各參數(shù)的下限， xZJKmax為中間框各參數(shù)的上限； xe＝0 表示“Π”型截面桁條，xe＝1 表示“I”型截面桁條，f(x) 為結(jié)構(gòu)的質(zhì)量表達式，F(xiàn)cr為結(jié)構(gòu)的極限載荷，F(xiàn)cr*為結(jié)構(gòu)的設(shè)計載荷。

3.2 Kriging 代理模型

在構(gòu)建代理模型前，需選取樣本點，通過后屈曲分析得到樣本點的結(jié)構(gòu)質(zhì)量和極限載荷。 試驗設(shè)計決定了構(gòu)建代理模型的樣本點規(guī)模和空間分布情況[9]。 本文采用試驗設(shè)計中的優(yōu)化拉丁超立方抽樣方法選取樣本點。 該方法具有良好的空間填充性和均勻性。 大多數(shù)試驗設(shè)計方法中試驗次數(shù)與因素個數(shù)成比例或指數(shù)增加，而優(yōu)化拉丁超立方抽樣方法的試驗次數(shù)等于水平數(shù)，適用于影響因素較多的情況，能顯著減少試驗次數(shù)[9]。

表2 不同蒙皮網(wǎng)格規(guī)模下結(jié)構(gòu)整體失穩(wěn)位移云圖Table 2 Collapse displacement fringes in different skinned grid scale

圖4 筋格內(nèi)的網(wǎng)格數(shù)為4×23Fig.4 The number of grids in the rib is 4×23

代理模型利用已知樣本點的輸出值來預(yù)測未知點的輸出值。 目前常用的代理模型有響應(yīng)面模型、徑向基函數(shù)代理模型、Kriging 代理模型和神經(jīng)網(wǎng)絡(luò)代理模型。 當問題的非線性程度較高時，多使用Kriging 代理模型[9]。 Kriging 模型是一種空間局部插值的方法，其基本思想是對已知函數(shù)進行加權(quán)來估計未知函數(shù)，并使得估計值與真實值之間的方差最小且數(shù)學(xué)期望相同的一種方法[15]。

設(shè)x0為未觀測的需要估值的點，x1，x2，…，xN為周圍的觀測點，觀測值相應(yīng)為y(x1)，y(x2)，…，y(xN)。 未觀測點的估計值記為~y(x0)，它由相鄰觀測點的已知觀測值加權(quán)求和得到[16]，如式(2)所示：

其中，λi為待定的加權(quán)系數(shù)。

采用R2檢驗作為代理模型近似能力的評價標準，R2的表達式如式(3)所示：

其中，[xi，y(xi)] 為待驗證樣本點，s(xi) 為代理模型xi處的預(yù)測值， ˉy 為待驗證樣本點輸出的平均值。 R2的取值范圍在[0，1]之間，且值越大(越接近于1)，擬合的精度越高，當R2＞0.9 時，可認為代理模型具有足夠高的擬合精度。

3.3 分步優(yōu)化策略

在蒙皮桁條結(jié)構(gòu)中，桁條是承受軸向力的主要結(jié)構(gòu)，它提升了蒙皮的受壓、受剪失穩(wěn)臨界應(yīng)力；中間框通過支撐縱向桁條以提高桁條穩(wěn)定性并增加結(jié)構(gòu)徑向剛度，從而提高蒙皮桁條結(jié)構(gòu)的承載能力。 范書群等[7]對大直徑火箭半硬殼結(jié)構(gòu)框-桁匹配性研究表明，桁條在兩框中間發(fā)生彎曲失穩(wěn)時，中間框也同時失穩(wěn)，此時結(jié)構(gòu)為最優(yōu)設(shè)計形式。 因此在第一步對桁條結(jié)構(gòu)參數(shù)進行優(yōu)化設(shè)計時，中間框設(shè)計得足夠強，以確保桁條在兩框之間發(fā)生失穩(wěn)破壞。 在第二步優(yōu)化時，基于桁條結(jié)構(gòu)參數(shù)優(yōu)化結(jié)果，對中間框結(jié)構(gòu)參數(shù)進行優(yōu)化，具體流程如圖5 所示。 其總體思想是利用試驗設(shè)計方法選取樣本點，為縮短設(shè)計周期，利用高性能計算機群對樣本點進行后屈曲分析得到其結(jié)構(gòu)質(zhì)量和極限載荷，根據(jù)樣本點的輸入-輸出關(guān)系構(gòu)建Kriging 代理模型。 對蒙皮桁條結(jié)構(gòu)參數(shù)進行合理分組后，采用基于多島遺傳算法(Multi-Island Genetic Algorithm， MIGA)和修正可行方向算法(Modified Method of Feasible Directions， MMFD)的組合優(yōu)化算法對蒙皮桁條結(jié)構(gòu)進行分步優(yōu)化。組合優(yōu)化算法能充分發(fā)揮全局優(yōu)化算法和梯度優(yōu)化算法的優(yōu)勢，在全局空間內(nèi)高效、準確地找到最優(yōu)解。

本文研究主要耗時集中在對構(gòu)建代理模型所需的初始樣本點的有限元計算上，采用高性能計算機對單個樣本點進行單次有限元計算耗時約20 分鐘。 為縮短計算耗時，提高優(yōu)化效率，本文基于高性能計算機群并行計算多個初始樣本點，對初始樣本點的后屈曲分析時間減少了約3/4，能極大地縮減設(shè)計周期。

圖5 優(yōu)化流程圖Fig.5 Flow chart of optimization

4 蒙皮桁條結(jié)構(gòu)輕質(zhì)優(yōu)化算例

表3 蒙皮桁條結(jié)構(gòu)參數(shù)變化范圍Table 3 Parameter range of the skinned truss structure

4.1 桁條選型及參數(shù)優(yōu)化

優(yōu)化后的桁條結(jié)構(gòu)參數(shù)相對初始值的變化量如圖6 所示，橫坐標中Mass 代表蒙皮桁條結(jié)構(gòu)的質(zhì)量。 可以看出，進行優(yōu)化后，相較于初始設(shè)計結(jié)構(gòu)，“Π”型截面桁條的結(jié)構(gòu)質(zhì)量減輕了6.3%，而“I”型截面桁條的結(jié)構(gòu)質(zhì)量減輕了2.9%，說明“Π”型截面結(jié)構(gòu)的承載效率更高。 因此，選取“Π”型截面桁條進行后續(xù)的優(yōu)化研究。

4.2 中間框參數(shù)優(yōu)化

圖6 優(yōu)化后的桁條結(jié)構(gòu)參數(shù)相對變化量Fig.6 Relative change of structural parameters of optimized trusses

對中間框結(jié)構(gòu)參數(shù)進行優(yōu)化后，優(yōu)化后的蒙皮桁條結(jié)構(gòu)參數(shù)、質(zhì)量和極限載荷相對初始值的變化量如圖7 所示。 可以看出，優(yōu)化后的蒙皮桁條結(jié)構(gòu)相較于初始設(shè)計結(jié)構(gòu)，質(zhì)量減輕16%，減重效果明顯。 優(yōu)化后蒙皮桁條結(jié)構(gòu)位移-載荷曲線如圖8 所示，結(jié)構(gòu)極限載荷為7.01×104kN，大于結(jié)構(gòu)的設(shè)計載荷， 說明優(yōu)化后的蒙皮桁條結(jié)構(gòu)的承載能力滿足要求。

圖7 優(yōu)化后的蒙皮桁條結(jié)構(gòu)參數(shù)相對變化量Fig.7 Relative changes of structural parameters of the skinned trusses

圖8 蒙皮桁條結(jié)構(gòu)位移-載荷曲線Fig.8 Displacement-loading curve of skinned trusse structure

5 結(jié)論

本文針對蒙皮桁條結(jié)構(gòu)輕質(zhì)優(yōu)化設(shè)計中變量類型眾多且有限元分析耗時長的問題，采用基于Kriging 代理模型和組合優(yōu)化算法的分步優(yōu)化策略開展優(yōu)化研究。 主要結(jié)論如下：

1) 計算耗時隨著蒙皮網(wǎng)格規(guī)模增加而增加。當蒙皮網(wǎng)格達到一定規(guī)模時，提升網(wǎng)格規(guī)模，極限載荷的計算結(jié)果相差較小。 綜合考慮計算精度和效率，采用筋格內(nèi)網(wǎng)格數(shù)量為4×23 的蒙皮網(wǎng)格；

2) 采用高性能計算機群對初始樣本點進行并行計算，計算耗時減少了約3/4，極大地縮減了設(shè)計周期；采用組合優(yōu)化方法能充分發(fā)揮全局優(yōu)化算法和梯度優(yōu)化算法的優(yōu)勢，便于在全局空間更快找到最優(yōu)解，二者綜合使優(yōu)化效率得到提升；

3) 相較于“I”型截面桁條的蒙皮桁條結(jié)構(gòu)，“Π”型截面結(jié)構(gòu)承載效率更高；

4) 相較于初始設(shè)計，優(yōu)化后的蒙皮桁條結(jié)構(gòu)減重16%，且滿足承載能力要求，驗證了所提方法的有效性。

參考文獻(References）

[ 1] 梅勇， 雷勇軍， 馮韶偉， 等. 超靜定捆綁火箭傳力路徑的組合優(yōu)化策略[J]. 南京航空航天大學(xué)學(xué)報， 2016， 48(06)： 909-916.

Mei Y， Lei Y J， Feng S W， et al. Combination optimization strategy for load trans-path of hyper-static strap-on launch vehicle [J]. Journal of Nanjing University of Aeronautics ＆ Astronautic， 2016， 48(06)： 909-916.(in Chinese)

[ 2] Hao P， Wang B， Li G， et al. Hybrid optimization of hierarchical stiffened shells based on smeared stiffener method and finite element method[J]. Thin-Walled Structures， 2014，82： 46-54.

[ 3] Quinn D， Murphy A， Glazebrook C. Aerospace stiffened panels initial sizing with novel skin sub-stiffening features[J]. International Journal of Structural Stability and Dynamics，2013， 12(5)： 531-542.

[ 4] Wang B， Tian K， Hao P， et al. Load-carrying capacity and imperfection-sensitivity analysis of hierarchical stiffened panels[J]. Solid Rocket Technology， 2014， 37： 408-412.

[ 5] Wang B，Hao P，Li G，et al. Generatrix shape optimization of stiffened shells for low imperfection sensitivity[J]. Science China Technological Sciences， 2014， 57(10)： 2012-2019.

[ 6] 郝鵬. 面向新一代運載火箭的網(wǎng)格加筋柱殼結(jié)構(gòu)優(yōu)化研究[D]. 大連：大連理工大學(xué)，2013.

Hao P. Optimum Design of Stiffened Shell Structures for New Generation Launch Vehicle [D].Dalian： Dalian University of Technology， 2013. (in Chinese)

[ 7] 范書群，戴政，黃誠，等. 大直徑火箭半硬殼結(jié)構(gòu)框-桁匹配性設(shè)計[J]. 強度與環(huán)境， 2015， 42(3)： 34-41.

Fan S Q， Dai Z， Huang C， et al. Matching design of framestringer for large scale semigrid launch vehicle cylindrical structure[J]. Structure ＆ Environment Engineering， 2015，42(3)： 34-41. (in Chinese)

[ 8] 張鐵亮， 丁運亮. 復(fù)合材料加筋壁板的結(jié)構(gòu)布局優(yōu)化設(shè)計[J]. 南京航空航天大學(xué)學(xué)報， 2010， 42(1)： 8-12.

Zhang T L， Ding Y L. Structural layout optimization of composite stiffened panel[J]. Journal of Nanjing University of Aerona Utics ＆ Astronautics， 2010， 42(1)： 8-12. (in Chinese)

[ 9] 張柱國，姚衛(wèi)星，劉克龍. 基于進化Kriging 模型的金屬加筋板結(jié)構(gòu)布局優(yōu)化方法[J]. 南京航空航天大學(xué)學(xué)報，2008(04)： 497-500.

Zhang Z G， Yao W X， Liu K L. Configuration optimization method for metallic stiffened panel structure based on updated kriging model[J]. Journal of Nanjing University of Aerona Utics ＆ Astronautics， 2008(04)： 497-500. (in Chinese)

[10] 郝鵬， 王博， 鄒威任， 等. 基于RBF 模型的蒙皮桁條結(jié)構(gòu)減輕孔優(yōu)化[J]. 固體火箭技術(shù)， 2015， 38(05)： 717-721.

Hao P， Wang B， Zou W R， et al. Optimum design of lightening holes for skin-stringer structures based on RBF model[J]. Journal of Solid Rocket Technology， 2015， 38(05)：717-721.(in Chinese)

[11] 張希. 基于顯式有限元方法的蒙皮桁條結(jié)構(gòu)穩(wěn)定性分析方法[C]/ /北京力學(xué)會第十九屆學(xué)術(shù)年會論文集， 北京，2013： 495-496.

Zhang X. Stability analysis method of skin truss structure based on explicit finite element method [C]/ / Proceedings of the 19thAnnual Academic Meeting of the Beijing Institute of Mechanics， Beijing， 2013：495-496. (in Chinese)

[12] 郝鵬， 王博， 李剛， 等. 基于代理模型和等效剛度模型的加筋柱殼混合優(yōu)化設(shè)計[J]. 計算力學(xué)學(xué)報， 2012， 29(4)： 541-548.

Hao P， Wang B， Li G，et al. Hybrid optimization of grid-stiffened cylinder based on surrogate model and smeared stiffener model[J]. Chinese Journal of Computational Mechanics，2012， 29(4)： 541-548. (in Chinese)

[13] 李洋， 趙斌， 龍連春. 網(wǎng)格疏密對加筋圓柱殼有限元計算的影響分析[C]/ /北京力學(xué)會第17 屆學(xué)術(shù)年會論文集，北京，2011： 392-393.

Li Y， Zhao B， Long L C. Influence of reinforced cylindrical shell’s grid density on finite element calculation[C]/ /Proceedings of the 17thAnnual Conference of Beijing Mechanics Society， Beijing， 2011： 392-393. (in Chinese)

[14] Dominguez A， Stiharu I，Sedaghati R. Practical design optimization of truss structures using the genetic algorithms[J]. Research in Engineering Design， 2006(17)： 73-84.

[15] 安俊虎. 空間燃料站貯箱結(jié)構(gòu)優(yōu)化設(shè)計[D]. 大連：大連交通大學(xué)， 2018.

An J H. Optimum Structural Design for Tank of Space Fuel Station[D]. Dalian：Dalian Jiaotong University， 2018. (in Chinese)

[16] 韓忠華. Kriging 模型及代理優(yōu)化算法研究進展[J]. 航空學(xué)報， 2016， 37(11)： 3197-3225.

Han Z H. Kriging surrogate model and its application to design optimization：A review of recent progress[J]. Acta Aeronautica ET Astronautica Sinica， 2016， 37(11)： 3197-3225. (in Chinese)