厚薄通用四邊形平板殼元在薄壁結構加筋布局優化中的應用

王 謙，丁曉紅，張 橫

（上海理工大學 機械工程學院，上海 200093）

0 引 言

舵面作為彈體的重要部件，是飛行器主要承力部件之一，具有輕質、高承載的設計要求［1］。舵面結構由薄壁板殼結構組成，單純薄壁板殼結構的剛度、強度、振動等力學性能不佳，在實際應用中一般采用加筋結構形式。布置加強筋能夠以較小的質量增加為代價，大幅提高板殼結構的靜動態性能。因此，如何布局加強筋成為了加筋設計中的重點問題。

近年來，研究人員將結構拓撲優化引入舵面結構加強筋分布的設計中，提升了舵面結構的各項性能［2-4］。但是基于變密度法（density method）尋求最優材料布局進而間接確定加強筋位置的方法需要特殊的后處理手段，才能判定加強筋的實際位置、走向、厚度等，但后處理過程存在一定的主觀性，會影響優化設計結果的最優性。丁曉紅等［5-6］提出了一種啟發式幾何顯式加筋拓撲優化方法“自適應成長法（adaptive growth method）”，該方法基于自然界分枝系統成長、分歧與退化機理發展而來。相較變密度法等幾何隱式方法，自適應成長法直接以加強筋截面幾何參數作優化設計變量，具有優化效率高、優化結果清晰、可制造性強等優點，近年來被廣泛運用于各類板殼結構加筋靜動態優化設計中。

自適應成長方法通過有限元方法進行結構分析，在迭代中以筋板截面參數變化表示“成長”和“退化”。季學榮［7］利用自適應成長法對汽車發動機罩板加筋布局進行了設計。侯劍云［8］同樣受自然界生物自適應成長機理啟發，提出基于連續階躍參考應力的軟殺法（soft kill optimization，SKO）。這兩項研究的結構分析部分利用商業有限元分析軟件，采用差分法計算單元靈敏度，因此，計算精度受到限制且計算效率較低；Shen 等［9-10］基于自適應成長法提出了3D 箱型結構的剛度、阻尼協同優化方法，結構分析部分利用Long 等［11］提出的基于廣義協調理論的平板型矩形殼元GCR24 處理的，該單元由平面薄膜元GR12 和薄板彎曲元GPL-R12 單元耦合構建。GCR24 單元用于幾何顯式板殼加筋拓撲優化時具有列式簡單、計算效率高的優勢，同時在筋板厚度變化過程中避免了剪切閉鎖問題，但GCR24 單元是矩形單元，不能滿足舵面等復雜結構的有限元分析對網格靈活性的需求。

綜上所述，利用自適應成長法對舵面結構的加筋布局進行拓撲優化，構建一種厚薄通用、高準確度的四邊形平板殼元十分必要。本文通過耦合GQ12膜元和MITC4 板元構造了四邊形平板殼元并構建基結構，隨后引入自適應成長法對舵面結構進行加筋布局拓撲優化，以此實現結構的輕量化設計。

1 GQ12 膜元與MITC4 板元耦合的四邊形平板殼元

1.1 具有旋轉自由度的四邊形膜元GQ12

為克服低階位移型單元計算精度差的問題，Long等［12-13］采用廣義協調條件構造了具有旋轉自由度的四邊形膜元GQ12，該單元列式簡單，能通過任意四邊形分片檢驗收斂單元，計算精度較高。GQ12引入剛體轉角作為單元節點的旋轉自由度，每個節點有兩個線自由度和一個平面內轉角自由度，每個節點的自由度為

式中：ui、vi為節點線位移分量；θi為節點附加平面內轉角位移分量。

GQ12膜元的位移場由兩個部分組成：

式中：u0=[u0v0]T，為雙線性協調位移場；uθ=[uθvθ]T，為由單元節點角位移分量θi引起的附加位移場。附加位移場為

式中：(ξ，η)為單元的自然坐標系。

將附加剛體轉角位移而產生的單元邊界位移沿單元邊界設為三次分布，并對每一邊引入廣義協調條件。待定參數αi、βi(i=1，2，3，4)的詳細推導過程參見文獻［12］，代入式（3）可得GQ12 單元的總位移場，表示為

式中：N0i、Nuθi、Nvθi為單元對不同位移場的形函數。

式中：ξ、η和ξi、ηi分別為單元和單元節點的自然坐標系坐標。

GQ12的應變場為

單元剛度矩陣為

式中：B為單元應力-應變矩陣；D為彈性系數矩陣；|J|為雅可比矩陣行列式；t為單元厚度。

1.2 剪應變混合插值板元MITC4

基于一階剪切變形理論的Mindlin 板元對位移場和轉角場各自獨立插值，由于單元邊界只要求C0類連續，因此滿足網格靈活性的需求，以及舵面等復雜外形結構的分析需要。一階剪切變形理論放棄了Kirchhoff 理論的直法線假設，因此對于橫向剪切變形明顯的中厚板的分析，Mindlin 板元的預測效果相比Kirchhoff板元更好。但是，Mindlin板用于薄板時存在剪切閉鎖和應力振蕩等問題，不能滿足幾何顯式加筋拓撲優化對有限元分析的要求。

縮減積分和選擇積分方法列式簡單，可解決剪切閉鎖的問題，但根據約束類型等的不同，會帶來零能模式問題?；旌喜逯凳且环N比較理想的構造厚薄通用板元的方法，通過對單元面內應變和橫向剪切應變使用不同的插值函數，避免剪切閉鎖和零能模式問題。Bathe等［14-15］基于Mindlin板理論提出剪應變混合插值板元（mixed interpolated tensorial component，MITC4），計算結果較好，滿足復雜結構的工程分析需要，如圖1所示。

圖1 板單元位移場Fig.1 Displacement field of plate element

由圖1可得，板單元的橫向剪切應變為

式中：θx、θy為x軸和y軸轉角位移分量；ω為橫向位移。

Mindlin板理論中，對于系統總位能公式中的彎曲應變能項，MITC4 單元使用與基于位移方法的Mindlin板單元相同的插值方法，廣義位移通過相同的形函數獨立插值為

式中：ωi、θxi、θyi分別為單元各節點處橫向位移分量、x軸和y軸轉角位移分量。

但是，對于橫向剪切應變能項，如果采用相同的插值形式，會在薄板情況下出現剪切閉鎖問題。Bathe 等［14-15］提出的MITC4 單元方案中，構造了一種在單元邊和單元內滿足剪應變一致的插值函數，插值形式為

式中：γAξz、γBηz、γCξz、γDηz分別為MITC4 單元在ξ-η-z坐標系下A、B、C、D四邊中點的橫向剪應變。

如圖2 所示，對于x-y平面內任意凸四邊形單元I-J-K-L，γξz、γηz的計算式為

圖2 任意四邊形單元Fig.2 Quadrilateral element

式中：detJ為雅可比矩陣行列式，

則有

同時，由于剪應變表達式定義在ξ-η-z坐標系下，需要轉化到x-y-z坐標系下，轉換形式為

式中：α為ξ軸與x軸夾角；β為η軸與y軸夾角。

1.3 驗證算例

為驗證GQ12 膜元與MITC4 板元耦合平板殼元的分析精度，設計一個四邊形折板，其幾何尺寸如圖3所示。折板一側約束，另一側一點受集中力載荷F=-1 000 N 的作用，取結構E=210 GPa，v=0.3，分別對5 mm厚度及10 mm厚度折板算例求解結構應變。

圖3 四邊形折板算例Fig.3 Example of quadrilateral folded-plate

5 mm 厚度折板GQ12 膜元與MITC4 板元耦合平板殼元自編Matlab 程序的分析結果如圖4（a）所示，采用商業有限元分析軟件Ansys的分析結果如圖4（b）所示。不同厚度折板加載節點和最大位移節點Z軸方向位移數值對照見表1和表2。

表1 5 mm厚四邊形折板Z軸方向位移Tab.1 Z-axis displacement of the quadrilateral foldedplate with 5 mm thickness

表2 10 mm厚四邊形折板Z軸方向位移Tab.2 Z-axis displacement of the quadrilateral foldedplate with 10 mm thickness

圖4 四邊形折板分析結果對比Fig.4 Analysis results comparison of the quadrilateral folded-plate

對以上算例結果進行對比，考察加載節點與最大位移節點Z軸方向位移值，相同約束和載荷條件下，5 mm 厚度折板與10 mm 厚度折板自編程序GQ12+MITC4 耦合平板殼元與商業有限元軟件分析結果誤差均＜4.5%，說明GQ12+MITC4耦合平板殼元具備厚薄通用性，計算精度良好。

2 自適應成長法

在自適應成長法結構加筋拓撲優化中，采用基結構法建立優化幾何模型，其優點在于避免優化迭代過程中網格的重劃分。參考承載和約束條件確定“種子點”的位置，如圖5（a）所示。從“種子點”生長出來的加強筋根據一定規則成長或退化；生長到一定尺度可以“分歧”，再從“分歧點”繼續生長，在下步迭代中參與成長與退化；退化到一定尺度則“消失”，如圖5（b）—圖5（c）所示。這個過程在優化迭代中反復進行，直至滿足收斂條件。

圖5 加強筋拓撲優化的自適應成長過程Fig.5 Adaptive growth process of stiffener topology optimization

2.1 數學模型及靈敏度分析

靜載荷下，薄壁結構自適應成長法加筋拓撲優化問題的數學模型可以表述為：在給定的體積約束下，選擇加強筋厚度作為設計變量，以結構在載荷條件下的總應變能為目標函數，逐步尋求最優解。其統一的數學模型可以表達為

式中：X為設計變量；xi為第i個活動筋板單元的厚度；n為加強筋單元的總數；Φ(X)為目標函數，即總應變能；V、V0分別為整體結構總體積和結構初始體積；η為體積約束因子；xmin、xmax為設計變量xi的取值下限和上限；g(X)為約束函數。

薄壁結構的柔度關于加強筋厚度的靈敏度計算公式為

式中：ui為第i加強筋單元的位移矢量；ki為第i加強筋單元的剛度矩陣。

2.2 優化設計流程

基于優化數學模型和自適應成長法基本原理，以移動漸近線法（method of moving asymptotes，MMA）作為迭代更新準則對結構分布進行尋優，具體流程如圖6所示。

圖6 自適應成長法優化設計流程Fig.6 Flow chart of optimization design of adaptive growth method

1） 建立基結構，設置設計域。

2） 選取“種子點”，設置優化參數。參考結構的載荷和支撐情況，選取基結構上若干節點作為“種子點”，給定加強筋的初始厚度x0、分歧臨界值xb、體積約束因子η、收斂容差ε以及最大迭代次數N。

3） 靈敏度分析。對結構進行有限元分析，通過解析法計算活動筋板單元的靈敏度。

4） 成長和退化。判斷收斂狀態，如果滿足筋板成長條件，通過MMA 算法更新xi；如果筋板厚度達到分歧臨界值xb，該單元進行分歧；如果滿足退化條件，該單元退化。

3 舵面結構加筋布局優化

3.1 舵面基結構構建

根據舵面幾何形狀，在優化設計前構建拓撲優化基結構，舵面外層為蒙皮結構，定義厚度為0.5 mm。選擇舵面內部空間區域作為設計域，定義基結構初始厚度為0.5 mm，如圖7 所示。由于舵面的幾何特征，構成基結構的單元為任意凸四邊形單元，GQ12 膜元和MITC4 板元耦合構造的四邊形等參平板殼元適用于舵面基結構的有限元分析。

圖7 舵面優化初始基結構Fig.7 Initial ground structure of rudder structure for optimization

3.2 基于自適應成長法的加筋布局

根據舵面的實際工作狀態，對其約束施加于舵軸，受單側橫向均布載荷。因此，本文以舵面舵軸處約束點為自適應成長法的“種子點”，以整體結構應變能最小為設計目標，骨架結構體積為約束條件，對舵面進行加筋拓撲優化，尋找最優的骨架結構分布。優化的數學模型為

式中：T為設計變量，為每個活動筋板的厚度；U(T)為目標函數，即整體結構應變能；V、V0分別為加筋舵面結構總體積和初始實體舵面結構體積；η=20%為體積約束因子。

經過拓撲優化設計得到的舵面內部最優加筋布局如圖8 所示。圖中呈現出清晰的分枝結構，拓撲形態合理且具有較優的可制造性。

圖8 拓撲優化后的舵面加筋布局Fig.8 Rudder stiffener layout obtained by topology optimization

優化過程迭代曲線如圖9 所示。優化迭代過程中，應變能與初始應變能比值在前10個迭代步中下降較快，隨后迭代過程逐漸平緩直至收斂，優化效果理想，驗證了GQ12膜元和MITC4板元耦合四邊形等參平板殼元分析的準確性，自適應成長法的有效性及優化結構的合理性。

圖9 舵面結構加筋布局優化迭代曲線Fig.9 Iteration history curve of rudder stiffener layout

3.3 設計結果分析

為了評估加筋拓撲優化效果，利用耦合四邊形平板殼元建立舵面的傳統蒙皮-網格式骨架設計模型，如圖10所示。對文獻［1］中介紹的傳統骨架設計進行建模并分析，計算得到拓撲優化前后舵面的整體結構應變能，如圖11所示。

圖10 舵面傳統骨架設計Fig.10 Traditional rudder skeleton design

圖11 舵面結構優化前后整體結構應變Fig.11 Strain of overall structure before and after rudder structure optimization

對比傳統網格式骨架設計與加筋拓撲優化后舵面整體結構的應變能數值以及結構質量，結果見表3。由表3可見，優化后結構的應變能為17.28 J，較傳統設計25.80 J降低了33%，結構剛度得到增強；傳統骨架設計的質量為1.19 kg，優化后結構的質量降低了29%。

表3 舵面結構優化前后性能比較Tab.3 Performance comparison before and after rudder structure optimization

4 結束語

本文利用GQ12 膜元與MITC4 板元耦合構造了四邊形等參平板殼元，通過有限元分析實現了基于自適應成長法的不規則薄壁結構的加筋設計。隨后以輕質、高承載的舵面結構為例，獲得了靜載荷狀態下舵面內部的最優骨架布局。經過對優化前后舵面結構的性能對比，驗證了GQ12膜元與MITC4板元耦合平板殼元分析工具的準確性和實用性，同時驗證了自適應成長法對舵面優化設計的有效性。