陳超磊, 王志祥*, 雷勇軍, 王 婕
(1.國防科技大學空天科學學院, 長沙 410073; 2.空天任務智能規劃與仿真湖南重點實驗室, 長沙 410073;3.北京宇航系統工程研究所, 北京 100076)
隨著深空探測和載人登月等航天任務逐步深入,研制大型運載火箭成為各國競相追逐的熱點[1]。加筋圓錐殼結構作為大型運載火箭不同直徑艙段的連接部分因具有高的軸向、彎曲以及扭轉承載效率被廣泛應用。在研制階段,加筋圓錐殼結構的優化設計作為一項重要技術影響火箭的承載能力和結構穩定性。在運載火箭發射過程中,由于自身的慣性力、相鄰部段的反作用力和空氣動力等因素[2],加筋圓錐殼艙段結構承載性能對確保運載火箭安全發射升空起到重要作用。
國內外研究人員對加筋薄殼圓柱結構開展了大量研究。王博等[3]考慮了初始缺陷對薄壁結構承載力的影響,提出了一種雙層蒙皮加筋柱殼結構優化設計方法,與傳統網格加筋結構相比,大幅提高了結構的承載能力;趙振等[4]開展了加筋肋圓柱殼穩定性優化設計,給出了加筋肋圓錐殼結構優化設計的準則;為提高大型運載火箭框桁加強薄壁圓柱殼在加工誤差影響下的結構性能,王志祥等[5]通過基于探索策略和開發策略的并行序列近似優化方法,開展了考慮加工誤差框架加筋圓柱殼的優化設計,提高了加強薄壁圓柱殼后屈曲優化效率;Ye 等[6]基于試驗設計方法,采用多目標遺傳算法對復合加筋板結構進行優化設計;Zarei 等[7]研究了加筋方向、蒙皮層數和半頂角對復合材料加筋圓錐殼結構整體屈曲行為的影響;Hao 等[8]考慮結構缺陷對加筋圓錐殼承載性能的影響,并開展了3.35 m 復合材料加筋圓錐殼結構輕量化設計。現有研究多針對加筋圓柱殼開展結構優化設計,而運載火箭加筋圓錐殼艙段結構優化設計相關工作有待豐富。
在結構優化設計中,針對大型復雜結構有限元計算耗時長的問題,可將結構有限元分析當作黑盒處理,建立結構參數-極限載荷的輸入輸出非線性映射關系,提高大型復雜結構的分析優化效率。基于代理模型的優化設計可以選取較少數量的試驗設計,獲得較多的變量響應信息[9],因此,基于代理模型的優化方法被廣泛應用[10-12]。王文竹等[13]基于Kriging 代理模型對鼓式制動器穩定性進行優化設計;胥磊等[14]對加筋柱殼結構開展結構優化,采用Kriging 代理模型得到了11.91%減重的最優結構;郝鵬等[15]在加筋柱殼優化設計中給出了一種基于代理模型和等效剛度模型的混合優化策略;在運載火箭蒙皮桁架結構輕質化研究上,王志祥等[16-20]建立了大型運載火箭柱殼模型,提出了基于矩估計的ARBF 近似建模方法和基于近似模型和組合優化算法的序列近似方法。移動可變形組件[21](Moving Morphable Component,MMC)在蒙皮點陣等結構輕量化設計中具有顯著優勢,該方法以具有顯式幾何參數的、可在設計空間自由移動和變形的組件作為結構拓撲的單元,優化得到的結構具有無灰度單元、顯式幾何信息、設計變量少等優勢[22],可為加筋圓錐殼輕量化設計提供一定理論參考。Meng 等[23]提出了一種基于甲蟲翅鞘的仿生雙曲面點陣設計,通過設計具備壓-扭特性的點陣單胞,獲得了具有優異緩沖和吸能效果的抗沖擊點陣結構。
本文以運載火箭輕質優化設計為背景,以加筋圓錐殼級間段為研究對象,基于Python 語言建立加筋圓錐殼結構參數化模型,采用非線性顯式動力學算法[24]分析結構后屈曲特性,通過分析確定影響大型圓錐殼結構承載能力的關鍵設計參數,并定量分析不同加載速度以及網格劃分規模對結構有限元分析精度和速度的影響;分別采用基于徑向基函數代理模型和Kriging 代理模型序列近似優化方法開展結構優化設計,獲得最優化加筋圓錐殼結構設計參數。
加筋圓錐殼結構由蒙皮、桁條和中間框組成。蒙皮用于維持和支撐桁條(圖1);縱向承載桁條一般要求具有比較大的截面慣性矩,以提高圓錐殼結構抗失穩能力和整體承載能力;蒙皮內側沿高度方向等距分布3 個“Ω”形截面中間框,2 個“L”形截面端框,并且在蒙皮的外側沿環向均勻分布若干個“工”形截面的豎向桁條,端框、中間框、桁條的截面構型如圖2 所示。

圖2 上下端框、中間框以及桁條截面構型Fig.2 Method for calculating the effective response acceleration
采用Python 語言開展加筋圓錐殼結構參數化建模。蒙皮采用2Al2-T4 鋁合金材料,桁條、中間框及端框等其他結構采用TC4 鈦合金材料,相關材料參數如表1 所示。為了更準確地模擬蒙皮和桁條的失穩狀態,采用殼單元劃分網格[25]。加筋圓錐殼結構有限元模型劃分的總單元數為92 700,總節點數為100 692。加筋圓錐殼模型采用如下邊界條件和加載條件[17]:在加筋圓錐殼模型上彈性邊界面以及下彈性邊界面中心各建立一個參考點,并與相應端框邊界面剛性耦合,對下參考點進行固支約束,對上參考點約束除軸向位移外的其余自由度,同時在上參考點勻速施加35 mm 軸向位移。初始設計結構參數如表2 所示。

表1 材料相關參數表Table 1 Material related parameters

表2 初始設計結構參數表Table 2 Structure parameter list of the initial design
相比非線性隱式動力學法,非線性顯式動力學法能夠更穩健地描述加筋圓錐殼后屈曲狀態[26],與試驗一致性較好,但分析時長與網格大小以及加載速度有關。本文采用非線性顯式動力學算法開展結構有限元分析。
在采用顯式動力學方法進行加筋圓錐殼結構后屈曲分析時,加載速度會對結構失穩載荷與失穩模態產生影響,本文分別采用不同的加載速度對模型進行收斂性分析[25]。使用8 核3 GHz 主頻CPU 及32 GB 內存的計算機,表3 是加載速度分 別 為 700 mm/s、350 mm/s、175 mm/s、116.7 mm/s(即分別在0.05 s、0.1 s、0.2 s、0.3 s內將加筋圓錐殼結構的軸向位移從0 mm 加載到35 mm)結構的失穩云圖、極限載荷和計算耗時。由表3 可知,加載速度小于175 mm/s 時,減小加速度對極限載荷的影響不大,結構失穩波形基本穩定,但計算耗時成倍增加,故在后續結構優化中加載速度選取175 mm/s。

表3 不同加載速度下結構整體失穩情況Table 3 Overall structure instability at different loading speeds
圖3 與圖4 分別是對應加載速度的動能/內能比值隨加載位移變化曲線以及結構載荷-位移曲線。從圖3 可以看出,加載時間不同,即對結構采用不同軸壓加載速度,結構內部的動能/內能具有如下特點:①在加載的初始階段,由于結構從無載荷狀態突然進入有載荷狀態,曲線產生峰值,而后逐漸加載,動能/內能比值逐漸減小,且曲線趨于平穩;②在加載過程中,結構內部的動能/內能比值均在0.05 以內,故可認為是準靜態加載[27];③加載末期的峰值是由于載荷達到結構承載極限,結構突然發生整體失穩而導致的。從圖4 可以看出,隨著加載速度的增大,結構極限承載能力(即曲線的峰值)逐漸增大,對應的加載位移也逐漸增大,且在峰值點過后,曲線立即下探,說明結構在極限載荷后迅速發生整體失穩。

圖3 動能/內能隨加載位移變化曲線Fig.3 Kinetic energy/internal energy change curve with loading displacement

圖4 載荷位移曲線Fig.4 Load displacement curve
劃分不同的網格數量對結構后屈曲狀態的計算效率有影響。為了準確反映加筋圓錐殼結構整體失穩模態,需要分析蒙皮和桁條網格規模對顯式非線性動力學算法的影響,即分析得到結構的網格收斂性。
由表4 和圖5 可以看出,隨著劃分網格數量的增加,計算耗時逐漸增大,但極限載荷(圖5 中各曲線峰值)逐漸減小,且后屈曲現象出現時間逐漸前移。根據不同網格數的仿真結果,增大圓錐殼周向網格數可以更準確地反映蒙皮環向的失穩波形,相比增大沿母線方向網格數對極限載荷減小的影響程度更大,但耗時增加更明顯。

表4 筋格內不同規模網格計算結果Table 4 Calculation results of different scale grids in the reinforcement

圖5 不同規模網格載荷位移曲線Fig.5 Load displacement curves of grids of different scales
以15×30 筋格網格為標準,綜合考慮計算耗時和極限載荷計算結果,最終選擇如圖6 所示筋格內網格數為10×20 的圓錐殼網格。

圖6 筋格內劃分10×20 網格圖Fig.6 10×20 grid diagram divided in the ribs
加筋圓錐殼結構優化研究中,設計參數眾多,對于不同設計參數,在變化相同程度條件下,結構極限承載能力的變化程度各不相同。為深入了解各設計參數的權重,提高設計人員對各參數的認識,本文開展加筋圓錐殼結構極限承載能力影響因素分析,即對結構參數進行靈敏度分析。
采用優化拉丁超立方試驗設計方法在表5 所示結構參數設計空間中選取400 組樣本點,計算極限載荷,得到結構質量-極限載荷分布散點圖(圖7)。基于最小二乘法建立回歸模型,采用將對應系數進行歸一化處理的靈敏度分析方法[17],得到靈敏度分析指標(圖8)。

表5 加筋圓錐殼結構設計參數取值范圍Table 5 Value range of design parameters of stiffened conical shell structure

圖7 結構質量-極限載荷分布散點圖Fig.7 Structure mass-ultimate load distribution scatter diagram

圖8 靈敏度指標Fig.8 Sensitivity index
圖8 可知,結構參數對極限載荷和結構質量的靈敏度指標前8 項主要為桁條參數及部分中間框參數。影響極限承載力的前5 項參數為:腹板高度、桁條數量、下翼緣厚度、下翼緣厚度以及腹板厚度;影響結構質量的前5 項參數為:桁條數量、下翼緣厚度、腹板厚度、腹板高度及上翼緣厚度。相比其他結構參數,腹板高度、桁條數量、翼緣厚度結構參數及相互耦合效應是對結構承載性能和結構質量影響更為顯著的因素。
圓錐殼與圓柱殼二者的幾何構型不同,存在2種中間框安裝方式:①垂直圓錐母線加筋;②沿圓錐徑向加筋。對中間框2 種不同安裝方式進行計算,結果見表6。由表6 可知,2 種中間框加筋方式失穩波形大體一致,采用沿圓錐徑向加筋方式失穩波形覆蓋面積大,分布稀疏,局部失穩更為明顯。

表6 不同中間框加筋方式對結果整體失穩的影響Table 6 Influence of different stiffening methods on the overall instability of the result
為了方便后續分析,定義承載效率η為極限載荷Fcr與結構質量M的比值,如式(1)所示。
為更具一般性,分別計算400 組中間框采用沿圓錐徑向加筋方式與400 組中間框采用垂直母線加筋方式的樣本(其他各項參數對應一致),承載效率散點分布圖如圖9 所示,二者沒有顯著差異,即2 種不同中間框加筋方式對結構的承載效率無顯著影響。

圖9 承載效率分布散點圖Fig.9 Scatterplot of load-bearing efficiency distribution
傳統優化方法效率低下,成本高昂,采用代理模型優化是一種高效且經濟的方法[28]。代理模型的結構優化設計可有效減少加筋圓錐殼結構有限元分析次數,極大地提高優化效率[29]。
3.1.1 基于徑向基函數代理模型建模方法
徑向基函數[30](Radical Basis Function,RBF)是一種離散多元數據插值模型,具有結構簡單、計算效率高的特點,適用于復雜的高維非線性問題。
常用的徑向基函數有Gauss 基函數、三次函數、逆多二次函數以及多二次函數。本文鑒于Gauss 基函數在全局近似能力方面的明顯優勢[31],選用Gauss 基函數。
由核函數矩陣Φ及真實響應值向量Y可計算得到徑向基函數權系數向量W,如式(2)所示。
3.1.2 Kriging 代理模型建模方法
Kriging 代理模型是一種基于方差最小的無偏估計方法[32],利用已知樣本點函數值線性疊加可計算得到未知點函數值。基于統計學假設,Kriging 模型將未知函數y(x) 看成高斯靜態隨機過程的形式,如式(3)所示。
式中,β0為未知常數,表示隨機高斯過程y(x) 的期望;Z(x) 是均值為0,方差為σ的隨機高斯過程。采用拉格朗日乘數法求解最小化均方差優化問題,并基于KKT 條件,可以求解得到Kriging 代理模型表達式[33]。
3.1.3 基于代理模型序列加點的輕量化設計
基于代理模型序列加點的加筋圓錐殼結構輕量化設計流程如圖10 所示。

圖10 基于代理模型的優化流程圖Fig.10 Optimization flow chart based on proxy model
首先,根據結構優化設計的實際問題建立數學模型,確定設計變量、設計約束、設計空間以及設計目標;其次,在設計空間內選取初始樣本點并計算初始樣本點的目標函數和約束函數值;然后,構造代理模型后,優化求解最優點,并計算最優點的目標函數和約束函數值[34],將其加入到初始樣本集,再更新代理模型;最后,循環迭代,直到滿足預先設定的收斂條件時終止計算,輸出優化結果。收斂條件定義如式(4)所示。
式中,^fout為代理模型最優解,fout為最優解處的真實響應,N為當前迭代步數,Nmax為最大迭代步數,δ為最大容許誤差(取1× 10-3)。
代理模型近似精度評估指標通過R2關系數值衡量,見式(5)。R2接近1,則表示模型精度越高[35]。
式中,n為樣本總數,yi表示真實響應值,^yi表示代理模型預測值,ˉyi表示真實值的均值。
采用優化拉丁超立方設計方法[36],在表5 所示結構參數設計空間中選取400 組參數作為初始樣本集,在400 組初始樣本點集中隨機選取380組樣本點分別構建徑向基代理模型和Kriging 代理模型,用剩余20 組樣本點對代理模型的全局近似精度進行檢驗。為避免隨機因素的影響,該過程重復15 次,統計代理模型近似精度評估指標如圖11 所示。從計算結果來看,徑向基代理模型和Kriging 代理模型近似精度指標R2均在0.8 以上,表明2 種代理模型都具有良好的全局近似能力,可開展后續序列近似優化工作。

圖11 兩種典型代理模型近似精度Fig.11 Approximate precision of two typical proxy models
對構建代理模型建立優化模型如式(6)所示。
式中,x表示加筋圓錐殼結構設計參數,M(x) 及Fcr(x) 分別表示加筋圓錐殼結構設計參數取值為x時,對應的結構質量和極限載荷,表示設定的加筋圓錐殼結構最小極限承載力,xmax與xmin分別表示結構設計參數上下邊界取值,具體參數范圍詳見表5。
基于2 種典型代理模型的最優設計參數如表7 所示,迭代計算結果歷程如圖12 及圖13 所示。在初始訓練樣本集,優化模型及優化算法等其他條件均一致的前提下,在滿足極限承載力大于5×107N 的約束條件下,基于Kriging 代理模型的序列近似優化方法經過2 次迭代優化后得到的最優結構質量為3149.87 kg,與迭代前的初始結構最優化設計質量相比減重120.97 kg(3.70%);基于徑向基代理模型的序列近似優化方法經過23 次迭代優化得到的最優結構質量為2963.37 kg,與迭代前的初始結構最優化設計質量相比減重307.47 kg(9.40%)。相比采用基于Kriging 代理模型近似優化方法,基于徑向基代理模型近似優化方法能得到結構質量更輕、承載效率更高的加筋圓錐殼結構設計。

表7 加筋圓錐殼結構優化結果Table 7 Optimization results of stiffened conical shell structure

圖12 加筋圓錐殼結構質量迭代歷程與位移云圖Fig.12 Mass iterative process and displacement cloud map of stiffened conical shells

圖13 加筋圓錐殼結構極限承載力迭代歷程Fig.13 Iterative process of ultimate bearing capacity of stiffened conical shell structure
通過表7 所示參數,建立2 種質量最優加筋圓錐殼結構,得到如圖14 所示的加筋圓錐殼結構軸壓載荷-位移曲線及徑向位移云圖。基于徑向基代理模型優化得到的加筋圓錐殼結構在加載至30.97 mm 時進入后屈曲狀態,主要徑向位移集中在靠近結構上下端框區域,沿環向規則分布;相比前者,基于Kriging 代理模型優化得到的加筋圓錐殼結構則在加載至29.06 mm 時就進入了后屈曲狀態,主要徑向位移沿環向大致分布在4 塊區域,集中在桁條、中間框以及結構下端框部分。

圖14 加筋圓錐殼結構軸壓載荷-位移曲線及位移云圖Fig.14 Axial load-displacement curve and displacement nephogram of stiffened conical shells
隨著軸向位移的繼續加載,2 個結構軸向載荷迅速降低,結構發生整體的壓潰破壞。綜合圖15 和圖16(a),基于徑向基代理模型優化得到的加筋圓錐殼結構在發生壓潰破壞時最大應力達到1.290×103MPa,最大塑性應變為12.2%,主要分布在“工”字型截面桁條的上翼緣與腹板位置,位于加筋圓錐殼靠近上端框的兩中間框之間的部分,桁條發生彎曲和扭轉組合形式的破壞;綜合圖16(b)和圖17,基于Kriging 代理模型優化得到的加筋圓錐殼結構在發生壓潰破壞時最大應力達到1.286×103MPa,最大塑性應變為13.1%,主要分布在“工”字型截面桁條的上翼緣與腹板位置,位于加筋圓錐殼靠近下端框的兩中間框之間的部分,桁條發生彎曲和扭轉組合形式的破壞。

圖15 基于徑向基代理模型優化后的加筋圓錐殼結構等效塑性應變云圖Fig.15 Equivalent plastic strain nephogram of stiffened conical shell based on RBF surrogate model

圖16 基于2 種典型代理模型優化后的加筋圓錐殼結構應力云圖Fig.16 Stress nephograms of stiffened conical shell based on two typical surrogate models

圖17 基于Kriging 代理模型優化后的加筋圓錐殼結構等效塑性應變云圖Fig.17 Equivalent plastic strain nephogram of stiffened conical shell based on Kriging surrogate model
本文針對加筋圓錐殼結構輕質化設計中變量種類和數量眾多且有限元耗時長的問題,分別采用基于徑向基函數代理模型和Kriging 代理模型的序列近似優化方法開展研究。
1) 加筋圓錐殼承載性能分析時長隨網格劃分規模增加而增加,綜合考慮計算精度和耗時,采用筋格內網格劃分為10×20 的蒙皮網格。
2) 腹板高度、桁條數量和翼緣厚度結構參數及耦合效應是極限載荷和結構質量的關鍵影響因素,而中間框分別采用垂直母線和沿圓錐徑向安裝方式對加筋圓錐殼結構的軸壓承載性能的影響不顯著。
3) 基于徑向基代理模型比Kriging 代理模型得到的結構質量更輕、承載效率更高,其最優結構質量為2963.37 kg,與迭代前結構質量相比減重307.47 kg(9.40%)。