基于序列徑向基函數的運載火箭蒙皮桁條結構輕質優化*

王志祥，歐陽興，王 斌，張大鵬，雷勇軍

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

薄壁加筋圓柱殼在飛行器結構中得到了廣泛的應用，現有的運載火箭型號有80%的箭體艙段采用整體加筋和桁條加強的薄壁殼體結構，如運載火箭的級間段蒙皮桁條和貯箱等結構[1-2]。作為運載火箭主要承力部件，薄壁加筋柱殼輕量化設計可大幅提高火箭運載能力和節約發射成本，但目前我國箭體加筋柱殼結構的設計手段仍以工程手冊和經驗設計為主，再輔以有限元校核分析及地面試驗，缺乏具有針對性的高效優化設計方法[2-4]。特別地，對于我國目前正在研制的長征九號重型運載火箭，其箭體艙段直徑達9.5 m，起飛質量達4 000 t級，起飛推力達5×107N級[5-6]，大直徑大載荷的結構特點將給箭體結構輕量化設計帶來嚴峻挑戰。

對于承受軸向載荷的薄壁加筋圓柱殼結構，整體失穩往往先于結構強度破壞發生，是結構破壞的主要形式，這就使得軸壓穩定性成為設計該類結構的重點[7]。Donnell等[8-9]從求解非線性大撓度方程出發，針對薄殼柱殼結構后屈曲分析方法開展了相關研究。李慶亞等[10-11]從理論上對比研究了隱式弧長法、隱式動力學、顯式動力學三種常用分析后屈曲問題的數值方法，并基于顯式動力學方法對軸壓作用下薄壁圓柱殼結構開展后屈曲分析研究。王博等[12-17]采用非線性顯式動力學分析方法，對5 m及以下直徑運載火箭中的加筋殼結構進行分析，并基于代理模型和等效剛度模型進行優化設計。范書群等[4]基于非線性顯式算法對蒙皮桁條結構的中間框尺寸進行了優選設計，理論分析結果與實驗結果一致性較好，驗證了該分析方法的有效性。從公開文獻來看，目前主要是針對5 m及以下直徑的運載火箭薄壁加筋柱殼開展研究，當直徑達10 m級甚至以上時，薄壁加筋柱殼徑厚比將是以往結構2～3倍，其帶來的穩定性問題將更為突出，針對該問題的研究尚少。蒙皮桁條結構設計變量眾多，涉及離散的拓撲和連續的尺寸優化，其輕質優化問題是典型的混合整數非線性規劃問題，且隨著結構尺寸的跨越式提高，大直徑蒙皮桁條結構后屈曲分析與優化效率將受到極大挑戰[17]，因此亟須針對大直徑蒙皮桁條結構開展輕質優化設計研究。

本文以我國未來大型或超大型運載火箭薄壁加筋結構設計為背景，針對大直徑大載荷蒙皮桁條結構開展結構輕量化研究。基于Python語言建立蒙皮桁條結構全參數化模型，研究了不同加載速度對大直徑蒙皮桁條結構后屈曲行為的影響規律。針對蒙皮桁條結構輕質優化涉及離散的拓撲和連續的尺寸優化問題，將該優化問題轉化成含非線性性能約束的參數優化問題，提出了基于徑向基近似模型和多島遺傳算法加非線性二次規劃算法的混合優化策略求解轉化后的優化問題，獲得了可行的優化解。

1 軸壓后屈曲分析及優化問題描述

1.1 非線性顯式后屈曲分析

對于一個顯式動力學分析，運動方程如式(1)所示。

(1)

采用中心差分法對控制方程進行顯式的時間積分，應用一個增量步的動力學條件去計算下一個增量步的動力學條件。

at=(Ut-Δt-2Ut+Ut+Δt)/Δt2

(2)

Vt=(Ut+Δt-Ut-Δt)/(2Δt)

(3)

其中，Δt為時間增量。

將式(2)和式(3)代入式(1)中，則原運動方程可改寫成：

(4)

從式(4)可以看出，增量步結束時的狀態僅取決于該增量步開始時的位移、速度和加速度，在時間上“顯式地”向前計算位移、速度和加速度，因此不存在收斂性問題。

典型加筋柱殼結構軸壓下位移-載荷曲線如圖1所示，隨著載荷逐步增大，結構可能呈現出線性前屈曲—非線性后屈曲—壓潰破壞行為，結構進入線性前屈曲后仍可繼續承載，直至整體壓潰破壞，壓潰破壞后結構承載力急劇下降后趨于穩定[18]。相比于弧長法和隱式動力學方法，顯式動力學方法可以穩健地跟蹤軸壓薄壁結構的后屈曲及后壓潰路徑及行為，并能與試驗結果吻合較好[19]。

圖1 典型加筋柱殼軸壓位移-載荷曲線Fig.1 Schematic of load versus end-shortening curve for typical stiffened shell

1.2 蒙皮桁條結構有限元建模及其極限載荷計算

作為運載火箭典型的加筋柱殼承力結構，蒙皮桁條結構主要由端框、中間框、桁條和蒙皮組成。蒙皮內側沿高度方向等間距對稱分布4個“Ω”形截面的中間框，2個“L”形截面的端框，同時，蒙皮外側沿環向均勻分布“工”形截面的豎向桁條，端框、中間框及桁條截面形式如圖2所示。

蒙皮桁條結構主要通過桁條來提高結構軸壓承載能力，而蒙皮的作用主要是維形和支持桁條。且蒙皮在較小的載荷下會局部失穩和局部進入塑性狀態，結構仍能繼續承載，直至結構發生整體壓潰破壞，其極限承載力由結構整體失穩和后屈曲狀態決定。

本文基于Python語言對蒙皮桁條結構進行參數化建模[2]。由于蒙皮桁條結構常使用較薄的蒙皮和桁條，結構上體現為板殼特性，為了能夠準確模擬桁條局部截面平動和轉動，采用殼單元劃分網格，且桁條腹板高度方向劃分兩層單元[20]。模型邊界條件及加載條件如下：模型上下端面中心點各建立一個參考點，并分別與上下端面節點進行剛性耦合，在下參考點處進行固支約束，在上參考點處約束除軸向位移外的其余自由度，軸向勻速施加35 mm強迫位移。模型選用鋁合金材料，彈性模量為70 GPa，泊松比為0.3，密度為2.78×10-6kg/mm3，屈服應力440 MPa，強度極限為550 MPa，延伸率為6%。參照以往火箭蒙皮桁條結構設計方法，通過初步結構優化設計，確定蒙皮桁條結構初始設計結構參數，如表1所示。

表1 蒙皮桁條結構初始設計結構參數(1)表格中NHT表示桁條數量，其余設計變量名含義如圖2所示。Tab.1 Initial design parameters of skinned purlin structure

采用顯式非線性屈曲算法求解該結構極限承載能力時，失穩載荷和失穩模態均與加載速度相關，因此分別采用不同的加載速度對該模型進行分析，并觀察結構內部動能與內能的比值，確保加載過程為準靜態加載。計算采用4核2.9 GHz主頻CPU及8 GB內存的計算機。圖3和表2分別給出了加載速度為600 mm/s、300 mm/s、150 mm/s、100 mm/s(即分別在0.05 s、0.1 s、0.2 s、0.3 s內將位移值從0 mm施加到30 mm)時的載荷位移曲線、動能與內能的比值隨加載位移的變化曲線、結構失穩位移云圖、極限載荷以及計算耗時。

從圖3(b)可以看出，對于不同的加載速度，結構內部動能與內能的比值在加載過程中均小于5%，加載之初的峰值是由結構從無載荷狀態突然進入到有載荷狀態引起的，隨著加載繼續，該峰值逐漸減小，曲線趨于平穩；加載后期的峰值是由載荷達到結構承載極限后，結構突然發生整體失穩引起的，可以認為是準靜態加載。從圖3(a)及表2可以看出，隨著加載速度變小，結構極限承載能力(曲線的峰值)逐漸減小，對應的加載位移值也逐漸減小，且加載速度越小，峰值點后載荷位移曲線越陡，說明結構達到極限載荷后便迅速發生整體失穩。

(a) 端框(a) End frame

(a) 載荷位移曲線(a) Load-displacement curves

(a) 結構質量靈敏度(a) Sensitivity analysis of the structural mass

表2 不同加載速度下結構整體失穩情況Tab.2 Overall instability in different loading speed

經大量試算可知，當加載速度小于150 mm/s時，減小加載速度對極限載荷的影響小于4%，結構失穩波形相同，但計算耗時成倍增長。因此，綜合考慮計算精度和計算效率，在后續優化中加載速度可取150 mm/s。

1.3 蒙皮桁條結構優化問題描述

通常，結構設計優化包含拓撲、形狀和尺寸優化。蒙皮桁條結構優化同時涉及拓撲優化和尺寸優化。其中，拓撲優化與桁條數量相關，不同桁條數量決定了蒙皮桁條結構的拓撲形式，屬于離散結構拓撲優化；端框、中間框、桁條截面尺寸以及蒙皮厚度是在結構拓撲形式和形狀固定的情況下，搜索最優的截面尺寸，屬于結構尺寸優化。

考慮到蒙皮桁條結構中桁條沿環向分布的對稱性，本文分別將拓撲優化變量和尺寸優化變量轉化成整數變量和連續變量，從而將離散拓撲優化和尺寸優化問題轉化成混合整數非線性規劃問題。因此，蒙皮桁條結構優化問題可以描述為：

findx=[xc,xd]

minf(x)

s.t.g(x)≤0

(5)

式中，xc表示n個連續變量的向量集合。在蒙皮桁條結構優化中，xc表示連續的截面尺寸變量，而xd表示桁條數目變量。

2 基于徑向基函數近似模型的序列近似優化方法

2.1 徑向基函數近似模型方法

徑向基函數(Radial Basis Function,RBF)近似模型[21-22]是由三層結構構成的前向神經網絡：第一層為輸入層，節點個數等于輸入的樣本點個數；第二層為隱含層，節點個數取決于問題的復雜程度；第三層為輸出層，節點個數等于輸出變量的維數。徑向基函數網絡模型結構如圖4所示。

圖4 RBF網絡模型結構Fig.4 Model structure of radial basis function

基本徑向基函數的數學模型為：

(6)

研究表明[22-25]，Gauss分布函數在全局近似能力方面優勢明顯。因此將采用Gauss分布函數作為基函數，其表示如式(7)所示。

(7)

ω=Φ-1Y

(8)

2.2 基于近似模型和組合優化算法的序列近似優化方法

空間均布性及填充性高的樣本點能夠以更大的概率捕獲近似對象的特征信息，提高近似模型的近似精度[26]?；趦灮〕⒎皆囼炘O計方法選取初始樣本點，由于初始樣本點難以保證近似模型對最優解的準確預測，因此在優化的過程中需要通過加入新的采樣點對近似模型進行更新，逐步提高近似模型的近似精度。

提出的基于徑向基函數近似模型的序列優化方法如圖5所示。首先利用有限元計算模型獲得初始樣本點輸出值，形成訓練樣本點集，然后基于訓練樣本點集構造RBF近似模型并開展優化。基于代理模型和組合優化算法的序列近似優化設計分為：內層迭代優化和外層迭代優化。在內層優化中，采用多島遺傳算法(Multi-Island Genetic Algorithm,MIGA)和非線性二次規劃(Non-Linear Programming by Quadratic Lagrangian Programming,NLPQLP)算法對RBF近似模型進行尋優，收斂后進入外層優化；在外層優化中，將內層優化獲得的近似最優解作為新的采樣點，并調用有限元計算模型進行非線性后屈曲分析，如果有限元計算結果與RBF網絡預測的誤差滿足收斂條件，即小于0.1%，則優化結束，否則將該采樣點及有限元計算結果加入樣本點集中，更新近似模型，提高近似模型的局部近似能力，并再次進入內層優化，直至內外兩層都收斂。

圖5 基于徑向基函數近似模型優化流程Fig.5 Flowchart of sequence approximate optimization method based on RBF model

在基于近似模型的優化研究中，主要耗時集中在初始樣本點的計算上，利用銀河超級計算機對構建近似模型的初始樣本點進行多節點并行計算，近似模型的構建和基于近似模型優化的耗時僅為單次后屈曲分析的1/10左右，這意味著不僅允許選取更多的初始樣本點以提高近似模型全局近似精度，而且可以進行后續采樣以更新近似模型，提高其局部近似精度，大大縮短了優化設計周期。

3 蒙皮桁條結構輕質優化設計

3.1 結構參數靈敏度分析

蒙皮桁條結構輕質優化設計影響因素眾多，為了從中找出影響結構質量和極限載荷的主要因素，進而實現設計空間縮減并指導后續優化設計，針對蒙皮桁條結構進行參數靈敏度分析。采用優化拉丁超立方實驗設計方法選取2 400個初始樣本點，利用并行計算資源對不同結構參數下蒙皮桁條結構進行后屈曲分析。設計變量取值范圍見表3，圖6為樣本點在結構質量和極限載荷空間內的分布散點圖。

表3 設計變量取值范圍Tab.3 Design spaces of design variables

圖6 初始樣本點分布散點圖Fig.6 Distribution scatter diagram of initial sample points

根據上述2 400個樣本點的計算結果，首先將各個變量歸一化到[-1,1]中，然后基于最小二乘法建立二次回歸模型。

(9)

式中：β0，βi，βj和βij分別為對應項系數。

通過將系數歸一化，得到不同項對響應的貢獻率百分比，即對響應的靈敏度。

Ni=100Si/∑Si

(10)

式中：Si為對應項系數，且∑Ni=100。

基于上述靈敏度分析方法和工程設計經驗，對結構質量和極限載荷影響較大的前9個因素分別如圖7所示。由結果可知，桁條截面參數以及部分中間框參數對結構極限載荷影響較大，在后續結構優化設計中應重點考慮。

3.2 基于序列近似優化方法的蒙皮桁條后屈曲輕質優化

本節將基于如第2節所述的序列近似優化方法對蒙皮桁條結構進行后屈曲輕質優化。根據3.1節分析結果，并綜合工程設計經驗，選取如式(1)所示的參數作為優化變量，其他設計變量均取如表1所示的初始設計值，開展蒙皮桁條結構輕質優化設計。蒙皮桁條結構輕質優化數學模型可表示為：

findx

minM

xmin

(11)

為求解如式(11)所示的優化問題，擬采用2.2節所述的基于近似模型和組合優化算法的序列近似優化方法進行研究。通過優化拉丁超立方試驗選取1 000個樣本點，并調用并行資源計算，獲得相應的結構質量和極限載荷。為驗證所選樣本點空間均布性及徑向基函數代理模型全局近似能力，分別從1 000個樣本點中隨機選取980個樣本點構建徑向基函數代理模型，并用剩余的20個樣本點對代理模型全局近似能力進行檢驗，該過程獨立重復10次，分別評估每次近似精度指標R2值[23,26]，計算結果如圖8所示。由計算結果可知，極限載荷和質量的R2值均接近于1，表明所選樣本點具有較好的均布性，且構建的徑向基函數代理模型具有較高的全局近似能力，可用于后續的結構輕質優化研究。

圖8 徑向基函數代理模型近似精度Fig.8 Approximate accuracy of RBF model

基于序列近似模型的優化平臺如圖9所示，其中，經多次試算分析，為提高多島遺傳算法種群多樣性和收斂速度，多島遺傳算法中設置島數為20、每個島種群數為20，進化代數為40，交叉率為1.0，變異率為0.01，遷徙率為0.01；為提高NLPQLP收斂精度，非線性規劃算法設置最小步長為10-4，最大迭代次數為50。

圖9 基于序列近似模型的優化平臺Fig.9 Optimization platform based on sequence surrogate model

3.3 結果分析

圖10為優化后樣本點在結構質量和極限載荷空間內的分布散點圖，可以看出，每次迭代后近似最優點的極限載荷均在目標極限載荷附近，且結構質量逐漸減小。

圖10 優化后樣本點分布散點圖Fig.10 Distribution scatter diagram of optimal sample points

經過27次迭代，優化結果趨于收斂，目標函數及性能約束的外層迭代優化曲線分別如圖11(a)及圖11(b)所示，每輪迭代后結構整體失穩模態如圖11(a)所示。優化后結構質量為3 530 kg，極限載荷為7.028×104kN，相比于初始設計結構，在承載力滿足設計要求的情況下，結構減重273 kg。優化前后結構尺寸如表4所示。

(a) 結構質量迭代曲線(a) Iteration curve of structural mass

表4 初始設計與優化結果Tab.4 Comparison of initial design and optimal variables

4 結論

本文以大型運載火箭結構輕質設計為研究背景，針對大直徑大載荷蒙皮桁條結構開展后屈曲輕質優化研究。

基于Python語言對運載火箭級間段蒙皮桁條結構建立參數化模型，為獲得結構極限承載能力，采用顯式動力學方法對蒙皮桁條結構進行后屈曲分析計算，并分析加載速度對計算精度和效率的影響，確定用于后續優化分析的加載速度。

針對軸壓作用下蒙皮桁條結構優化變量多、計算量大的問題，首先，基于優化拉丁超立方實驗設計方法對結構參數靈敏度進行了分析，根據靈敏度分析結果，合理地選擇了對結構極限承載能力及質量影響顯著的設計變量，有效縮減了優化變量維數；然后，提出基于近似模型和多島遺傳及非線性二次規劃算法的序列近似優化方法，并對蒙皮桁條結構開展后屈曲輕質優化，獲得了工程可行的較優的蒙皮桁條結構設計方案。優化結果表明：相對于初始設計結構，優化后的結構有效減重273 kg，驗證了方法的有效性。

后續將在本文的研究基礎上開展基于不同桁條截面形式的蒙皮桁條結構輕質優化研究，并針對最優結構形式開展試驗驗證研究，為大型運載火箭蒙皮桁條結構工程設計提供參考。