基于Voronoi 序列采樣的加筋壁板優化設計

于寶石，王志祥，王 婕，李海陽，張大鵬

（1.國防科技大學空天科學學院，長沙 410073；2.空天任務智能規劃與仿真湖南省重點實驗室，長沙 410073；3.北京宇航系統工程研究所，北京 100076）

薄壁加筋結構由于其較高的比剛度和比強度而廣泛應用于航空航天結構［1?2］，如運載火箭的燃料貯箱［3］、導彈的級間段［4］等。此類結構在受載時失穩破壞往往先于強度破壞，所以其穩定性問題一直是研究的熱點。薄壁加筋結構直接關系到航空航天結構的服役安全和輕質化水平［5?6］，特別作為主承力結構，軸壓載荷下其結構分析與輕量化設計對降低飛行和發射成本具有重要意義。

對于薄壁加筋結構的輕量優化問題，眾多學者采用傳統的梯度類優化算法求解并獲得了較為理想的優化結果［7］，但對于高維設計變量、模型復雜的結構，傳統優化方法又面臨難以收斂到全局解、計算量過高導致優化效率較低等問題。針對上述問題，具備快速預測和算法適用性強等優點的代理模型技術受到廣泛關注，并被成功應用到結構優化設計中［8?10］。目前，常用的代理模型方法包括多項式響應面、Kriging 模型以及徑向基函數（Radial ba?sis function，RBF）等。Jin 等［11］詳細對比了多種常用的代理模型方法，指出徑向基函數較好兼顧了計算效率與計算精度。Gutmann［12］對徑向基函數進行了更為詳細的介紹。Mullur 等［13］為進一步提升模型近似精度，提出了一種改進的徑向基函數。Wang 等［14?17］基于徑向基函數近似模型對薄壁加筋柱殼后屈曲優化設計開展了大量的研究。郝鵬等［2?3，18?22］結合代理優化方法與等效剛度模型，對加筋柱殼進行優化設計，后續又基于缺陷敏感性分析開展優化研究。但上述研究的優化方法均依賴初始近似建模的全局精度，針對多極值問題，容易陷入局部最優解。龍騰等［23?24］針對飛行器多學科設計優化中傳統靜態代理模型方法難以保證全局近似精度與計算效率較低的問題，提出了一種動態徑向基函數代理模型的優化策略。王志祥等［25］將徑向基函數與多島遺傳、非線性二次規劃算法相結合，提出一種序列近似優化方法，并在蒙皮桁條結構的優化中取得理想結果。上述研究的序列采樣策略均未有效兼顧全局近似精度與局部開發效率，實現局部與全局加點的動態平衡。然而，針對復雜結構，如何進一步提升基于代理模型序列近似優化方法的全局最優性與優化效率，仍是當前研究的重點。

Voronoi 法是計算幾何領域內一種重要的幾何結構，常常被應用于全局區域檢測與路徑搜索［26?27］。Jiang 等［28］將Voronoi 法引入優化流程中，提出了一種分步輔助序列采樣方法，具有較好的全局近似精度。Xu 等［29］提出了基于Voronoi 法和交叉驗證的穩健誤差追蹤序列采樣方法，并對比幾種序列采樣方法，具有更高的優化效率與魯棒性。上述研究均表明，Voronoi 法在改善樣本點空間分布均勻性，進而提高近似模型全局近似精度方面具有顯著優勢，但將Voronoi 法應用到基于代理模型的序列優化中的研究相對較少。

為綜合平衡序列采樣的探索和開發性能，提高近似模型優化效率，本文開展了基于Voronoi 法序列采樣的研究，并成功應用到加筋壁板的優化設計中。首先，對軸壓載荷下加筋壁板結構的屈曲特性進行分析，基于Python 語言建立了結構參數化模型，分析了在不同腹板高厚比下網格規模對結構屈曲特性的影響規律，根據分析規律，結合結構自身特點，依據腹板高厚比動態調整網格劃分，平衡了計算精度和優化效率。然后，基于Voronoi 法與徑向基函數代理模型，提出了一種序列近似優化方法，針對加筋壁板結構，獲得了可行的優化解。最后，對比文獻［19］中的優化方法，結果表明本文方法具有更高的優化效率與全局精度。

1 加筋板有限元建模與優化問題

1.1 加筋壁板有限元建模與分析

如圖1 所示，加筋壁板模型主要由蒙皮、桁條與端框構成。蒙皮外側沿縱向分布3 根“幾”形桁條，內側上下兩端分布“L”形端框，桁條與端框截面形式及參數如圖2 所示?；赑ython 語言對高600 mm、寬600 mm 的加筋壁板結構進行參數化有限元建模，為模擬結構板殼特性，采用殼單元劃分網格。為獲得軸壓作用下結構的極限載荷，分別建立“L”形彈性邊界與端框對接。分別建立參考點與上下彈性邊界節點進行剛性耦合，下參考點固支約束，上參考點約束除軸向位移外的其他自由度。勻速施加10 mm 強制位移，綜合考慮分析效率與計算精度［30］，將加載速度設置為200 mm/s。軸壓載荷作用下，加筋壁板的破壞形式為對應載荷最高點結構的失穩模式，失效載荷即為位移?載荷曲線最高點處的極限載荷。模型采用鈦合金材料，其密度為4.45×106kg/mm3，彈性模量為110 GPa，泊松比為0.3，屈服應力825 MPa，強度極限895 MPa，延伸率為10%?；趥鹘y工程設計方法，加筋壁板初始設計結構參數如表1 所示。表中M、Fcr0分別表示初始設計結構質量和極限載荷，其他變量含義如圖2 所示。精度與效率之間的平衡，成為結構優化設計中不可忽視的問題。針對上述問題，本文對加筋壁板結構進行網格無關性分析，總結分析網格規模對結構承載性能和失穩模式的影響規律，進而實現兼顧計算效率與分析精度的動態有限元建模。

圖1 結構有限元模型Fig.1 Structural finite element model

表1 加筋壁板初始設計及變量取值范圍Table 1 Initial design of stiffened panels and range of variable values

圖2 加筋壁板結構組件截面尺寸Fig.2 Section dimensions of stiffened panel structural assemblies

圖3 與表2 分別給出了不同蒙皮單元大小時位移?載荷曲線、結構失穩云圖、極限載荷以及計算耗時。隨蒙皮單元尺寸的減小，極限承載能力逐漸降低并收斂到一穩定值。當蒙皮單元尺寸小于15 mm×15 mm 時，結構失穩波形基本相同，繼續加密網格對極限載荷的影響小于0.224%，而計算耗時增加

表2 不同蒙皮網格尺寸下結構整體失穩情況Table 2 Overall instability under different skin mesh sizes

圖3 不同蒙皮網格尺寸下結構的位移-載荷曲線Fig.3 Displacement-load curves of structures with different skin mesh sizes

針對薄壁板殼結構，單元網格劃分的疏密程度將嚴重影響結構軸壓載荷與失穩模態的分析精度［31］。如何在大規模的優化計算中保持模型計算71.8%。因此綜合考慮計算精度與效率，在后續的優化計算中蒙皮單元尺寸確定為15 mm×15 mm。

對應結構優化迭代過程中尺寸參數的變化，若采用統一的網格劃分精度將會造成個別樣本參數下網格過疏或過密，相應地導致仿真誤差較大或計算成本較高?；谏鲜鰡栴}，以表1 所示的初始結構尺寸為基準，在不改變桁條腹板橫截面積下（即hfb×tfb不變），分析不同腹板高厚比(αfb=hfbtfb)下網格密度比對加筋壁板承載能力和失穩模式的影響規律。

圖4 為不同網格密度下結構承載能力隨桁條腹板高厚比變化曲線及結構整體失穩時桁條位移云圖。其中U 為桁條位移。由圖可知，在αfb≤52時，不同網格密度下結構承載能力與失穩模態基本相同。但隨αfb繼續增加，結構承載能力下降，桁條腹板出現失穩波形。這主要是因為，在αfb較小時，增大腹板高度即提升了桁條橫向彎曲剛度，桁條臨界失穩應力提高，桁條呈現整體失穩模式。而當αfb逐漸增大后，軸壓下桁條腹板抵抗面外變形的剛度降低，桁條腹板會先于整體發生局部失穩，從而使得結構承載能力下降。且隨著腹板單元層數Nfb增多，承載能力下降趨勢更加明顯并趨于收斂。這進一步說明在αfb較大的情況下，較少的單元層數會導致計算精度嚴重下降，無法反映出極限承載力變化規律與桁條局部失穩模態。

圖4 不同網格密度下結構承載能力隨桁條腹板高厚比變化曲線Fig.4 Curves of structural load-carrying capacity with height-thickness ratio of webs under different mesh densities

基于上述分析，在優化歷程中根據桁條尺寸參數在桁條腹板高度方向上進行網格劃分如下

式中：Nfb表現為隨αfb變換的分段函數。當αfb≥52時，較少的單元劃分存在較高的計算誤差；當αfb≤52 時，3～5 層單元足夠滿足計算精度，過度的劃分又會降低計算效率。對比設定單一網格密度，依據腹板高厚比的網格劃分兼顧了計算效率與分析精度，為后續的優化迭代節約了大量計算成本。

1.2 優化問題

以桁條、端框截面尺寸參數及蒙皮厚度等13個結構尺寸參數為設計變量，以滿足設定承載能力為約束，開展加筋壁板輕質優化設計。其中，設計變量均為連續變量。針對本結構的優化問題，其數學模型可表示為

式中：x表示結構設計參數，取值范圍參考表1。M(x)及Fcr(x)為結構參數取值x時對應的結構質量和極限載荷。F?cr為目標極限載荷，設定為1.13×106N。

2 序列近似優化方法

針對上述加筋壁板的后屈曲輕質優化，考慮到問題具有多變量、高耗時等特性，本文以徑向基函數近似模型為主干，提出了一種基于探索和開發并行采樣的序列近似優化方法。

2.1 徑向基函數近似模型

不失一般性，單位立方體空間Ω=[0，1]m中，給定一組訓練樣本集

具體建模過程可參考文獻［20，35］?；诖藦较蚧瘮到颇Ｐ停奢^為準確地反映加筋壁板結構參數與極限承載性能的近似映射關系，進而開展后續優化設計。

2.2 序列近似優化方法

序列采樣策略的基本思想是根據當前樣本點及近似模型相關信息，獲得下一個（組）采樣點，以期提高近似模型全局/局部近似精度，引導近似優化算法高效快速搜索到最優解?；诖怂枷?，本文提出了基于Voronoi 法的序列采樣策略。其中，基于Voronoi 法的探索策略通過動態識別稀疏樣本點區域并添加采樣點，進而提高近似模型的全局近似精度；基于近似模型與組合優化算法的開發策略通過在潛在最優點處加點，并利用近似最優解加速算法收斂，進而提升近似模型的局部近似精度。

2.2.1 基于Voronoi 法的探索策略

基于Voronoi 法的探索策略可以歸納為以下幾個步驟：

式中：dom(xi，xj)是被點xi與xj之間連線的中垂線分割的封閉區域，點xi稱為多邊形ci的中心點，多邊形ci中任意一點距其中心點xi的距離都是最近的。

（2）識別加點最稀疏的泰森多邊形cs。泰森多邊形的邊界一般不規則，很難準確描述其形狀，Voronoi 法通過Monte Carlo 法構造泰森多邊形。向設計空間Rm內生成NRP個隨機點，NRP取值為

當Voronoi 圖為二維圖形時，識別cs過程如圖5 所示。其中，cs區域隨機點數量最多，覆蓋區域最大。

圖5 二維Voronoi 圖識別cs 示例Fig.5 Two-dimensional Voronoi diagrams identify cs exam?ples

（3）在cs中確定新加點xn+1。如前所述，多邊形cs可以用一組隨機點Ps描述。因此，新加點xn+1可以按照一定的采樣準則從隨機點集Ps中選取。為盡可能多地獲得關于目標函數的信息，cs中采集新點xn+1應距中心點xs盡可能遠，即

2.2.2 基于近似模型與組合優化算法的開發策略

基于近似模型與組合優化算法的開發策略可以歸納為以下步驟：

（1）利用當前訓練樣本集S，采用2.1 節中的數學方法構造RBF 近似模型。

（2）采用多島遺傳算法（Multi?island genetic algorithm，MIGA）與非線性二次規劃算法（Non?linear programming by quadratic Lagrangian pro?gramming，NLPQLP）進行內層迭代優化，對RBF模型進行尋優，迭代收斂后獲得一個近似最優解x?以及對應的近似最優目標函數值f?(x?)，并計算近似最優解x?的真實響應f(x?)，即

式中：M?(x)為結構質量的近似模型；F?cr(x)為極限載荷的近似模型。

（3）將近似最優目標函數值f?(x?)與真實響應值f(x?)代入收斂判斷，滿足收斂條件則結束優化，否則將此近似最優解x?作為新采樣點，與其真實的響應值一起更新到樣本集中，即S=S∪[x?，f(x?)]與n=n+1，進而提高近似模型的局部近似能力。

2.3 算法整體框架

本文提出的序列近似優化方法基于徑向基函數近似模型，在每輪迭代中通過探索與開發策略序列采樣，以期不斷提高近似模型的全局與局部精度，并通過并行計算新采樣點，實現加快算法搜索效率與縮短計算時長的目的。算法整體框架如圖6 所示。

圖6 序列近似優化流程圖Fig.6 Flowchart of sequence approximate optimization method

利用初始樣本點構建的近似模型需進行精度檢驗，通過檢驗模型的均值與R2值等確定滿足精度要求的最小初始樣本規模nmin。不滿足精度要求則增大樣本集規模重新采樣，直至構建出滿足近似模型精度的初始樣本集。在優化策略收斂判斷中，滿足如下任一條件時，停止優化。

（1）第k次與第k-1 次迭代的采樣點處有限元計算值之間的相對誤差，第k次迭代采樣點處近似模型預測值f?(x?k)與有限元計算值f(x?k)之間的相對誤差，均滿足給定的收斂標準ε，如式（12）所示。

（2）序列采樣達到最大迭代次數kmax。在迭代過程中，基于Voronoi 法全局加點側重于改善樣本點在空間中的均布性及填充性，進而提高了近似模型全局近似精度。RBF 近似模型預測潛在最優加點側重于進一步開發局部最優值區域，以期提升近似模型對該區域的局部近似能力。本文所提優化算法綜合利用了二者優勢，實現全局探索與局部開發的平衡，有效地提升了優化效率。

3 加筋壁板結構輕質優化設計

根據2.2 節所述的序列近似優化方法對加筋壁板結構開展輕質優化。本文采用優化拉丁超立方實驗設計方法生成30 個初始樣本點（即nmin=30），然后計算對應樣本點的極限載荷與結構質量，以此構建初始近似模型。考慮到初始近似模型在最優點處預測不夠準確，采用第2 節所述方法進行序列加點。

在相關優化算法的參數設置上，為提高MIGA種群多樣性與收斂速度，多島遺傳算法中設置島數為20，各島種群數為20，進化代數為40，交叉率為1.0，變異率與遷徙率均為0.01；為提高NLPQLP收斂精度，非線性規劃算法設置最小步長為10-4，最大迭代次數為50。

為了驗證本文方法的有效性，將本文方法與Wang 等［18?21］提出的序列近似優化方法相對比。圖7 所示為兩種方法在迭代過程中載荷相對誤差eFcr和質量相對誤差emass。其中，圖7（a）為本文方法誤差收斂曲線，在迭代20 次后，極限載荷與質量的相對誤差均小于0.01，誤差趨于收斂。對比圖7（b）文獻方法誤差收斂情況，在迭代32 次后，相對誤差趨于收斂，說明了本文優化策略能在迭代過程中更快的收斂，具有更高的優化效率。

圖7 優化歷程中的誤差Fig.7 Error in optimization process

圖8 所示為目標函數與性能約束函數的迭代歷程曲線，本文方法在經歷25 輪迭代后，結構極限承載能力滿足設計需求，且結構質量也趨于收斂。對比文獻方法，在得到幾乎相同的優化結果的情況下，本文方法收斂更快，進一步驗證了本文序列近似優化方法具有更高的優化效率。

圖8 加筋壁板結構質量及極限載荷迭代歷程Fig.8 Iteration curve of structural mass and limit load

圖9 為兩種優化方法中近似最優加點與已有樣本點的最小距離隨迭代次數的收斂曲線。對比曲線收斂情況，文獻［19］的優化方法在局部區域內迭代加點并快速收斂，而本文優化方法在每次迭代時先向樣本點分布稀疏區域加點，再于迭代中快速收斂到局部區域。這說明基于Voronoi 法動態序列加點策略有助于跳出局部最優區域，是在全局搜索后再預測最優解，有效提高了模型的預測精度，具有更好的潛在最優區域開發能力。

圖9 兩種優化方法新加點與已有樣本點最小距離迭代曲線Fig.9 Iterative curve of minimum distance between new point and existing sample point in two optimization methods

經本文序列近似優化后，在軸壓載荷加載到4.7 mm 時，結構整體失穩破壞，對應的極限載荷為1.134×106N，得到結構質量為9.810 kg。在承載力滿足設計要求下，結構減重4.742 kg，減輕初始設計結構質量的32.6%。表3 列出了所提方法和文獻方法得到最優設計的比較結果，采用本文優化策略，只經過25 次迭代，結果收斂，較文獻方法優化效率提升了40.5%。在最優解方面，兩種方法最優解相對誤差僅為0.028%。且由于本文方法具有更高的全局近似精度，在初始樣本點的選取上，僅選取了30 個樣本點，相較文獻中的優化方法減少了40%，進一步節約了計算成本。

表3 設計變量初始值與最優值Table 3 Comparison of initial desing and optimal variables

結合圖10 與文獻方法最優解處軸壓位移?載荷曲線及整體失穩時位移云圖對比情況，由圖可知，兩條曲線的線性前屈曲段與非線性后屈曲段幾乎重合，極限載荷近似相同，且結構整體失穩模式也基本相同。表4 對比了兩種方法獲取的優化結構在整體失穩時的應力S云圖與等效塑性應變（Equivalent plastic strain，PEEQ）云圖，本文方法獲得的結構應力分布與塑性形變與文獻方法得到的基本一致，桁條與蒙皮連接處應力較大，最大應力分別為961 MPa 和935 MPa，蒙皮中心與桁條局部區域進入塑性形變，蒙皮最大等效塑性應變分別為0.8%和0.7%，中間桁條最大等效塑性應變分別為2.2% 和2.3%。兩種方法對比結果驗證了本文提出的序列近似優化方法的有效性。

表4 加筋壁板結構整體失穩時應力云圖與等效塑性應變云圖Table 4 Stress and equivalent plastic strain cloud diagram of completely buckling stiffened panel

圖10 加筋壁板結構軸壓位移-載荷曲線及整體失穩時位移云圖Fig.10 Loading-displacement curve and displacement cloud diagram of completely buckling stiffened panel

4 結論

本文以加筋壁板結構設計為研究背景，開展了基于Voronoi 法的加筋壁板高效序列輕質優化方法研究。主要工作與結論可總結如下：

（1）探究了加筋壁板結構蒙皮與桁條腹板網格劃分對結構屈曲特性的影響規律，建立了網格劃分方式與桁條腹板高厚比的對應關系，兼顧了建模精度與計算效率。

（2）提出了一種基于Voronoi 法的序列采樣策略，實現了探索和開發動態平衡的并行采樣。并將此方法應用于加筋壁板的后屈曲輕質優化中，結果表明，相較于初始設計，在滿足極限載荷約束的情況下，優化后結構有效減重了32.6%，驗證本文方法的有效性。

（3）將本文提出的序列近似優化方法與現有的僅添加近似最優點方法進行了對比，結果表明，本文方法提高了模型的預測精度，有效縮減了初始樣本集的規模，降低了迭代次數，提升了優化的效率，并且具有更好的全局探索與潛在最優區域開發能力。