厚薄通用四邊形平板殼元在薄壁結(jié)構(gòu)加筋布局優(yōu)化中的應(yīng)用

王 謙，丁曉紅，張 橫

（上海理工大學(xué) 機(jī)械工程學(xué)院，上海 200093）

0 引 言

舵面作為彈體的重要部件，是飛行器主要承力部件之一，具有輕質(zhì)、高承載的設(shè)計(jì)要求［1］。舵面結(jié)構(gòu)由薄壁板殼結(jié)構(gòu)組成，單純薄壁板殼結(jié)構(gòu)的剛度、強(qiáng)度、振動(dòng)等力學(xué)性能不佳，在實(shí)際應(yīng)用中一般采用加筋結(jié)構(gòu)形式。布置加強(qiáng)筋能夠以較小的質(zhì)量增加為代價(jià)，大幅提高板殼結(jié)構(gòu)的靜動(dòng)態(tài)性能。因此，如何布局加強(qiáng)筋成為了加筋設(shè)計(jì)中的重點(diǎn)問(wèn)題。

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

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

綜上所述，利用自適應(yīng)成長(zhǎng)法對(duì)舵面結(jié)構(gòu)的加筋布局進(jìn)行拓?fù)鋬?yōu)化，構(gòu)建一種厚薄通用、高準(zhǔn)確度的四邊形平板殼元十分必要。本文通過(guò)耦合GQ12膜元和MITC4 板元構(gòu)造了四邊形平板殼元并構(gòu)建基結(jié)構(gòu)，隨后引入自適應(yīng)成長(zhǎng)法對(duì)舵面結(jié)構(gòu)進(jìn)行加筋布局拓?fù)鋬?yōu)化，以此實(shí)現(xiàn)結(jié)構(gòu)的輕量化設(shè)計(jì)。

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

1.1 具有旋轉(zhuǎn)自由度的四邊形膜元GQ12

為克服低階位移型單元計(jì)算精度差的問(wèn)題，Long等［12-13］采用廣義協(xié)調(diào)條件構(gòu)造了具有旋轉(zhuǎn)自由度的四邊形膜元GQ12，該單元列式簡(jiǎn)單，能通過(guò)任意四邊形分片檢驗(yàn)收斂單元，計(jì)算精度較高。GQ12引入剛體轉(zhuǎn)角作為單元節(jié)點(diǎn)的旋轉(zhuǎn)自由度，每個(gè)節(jié)點(diǎn)有兩個(gè)線自由度和一個(gè)平面內(nèi)轉(zhuǎn)角自由度，每個(gè)節(jié)點(diǎn)的自由度為

式中：ui、vi為節(jié)點(diǎn)線位移分量；θi為節(jié)點(diǎn)附加平面內(nèi)轉(zhuǎn)角位移分量。

GQ12膜元的位移場(chǎng)由兩個(gè)部分組成：

式中：u0=[u0v0]T，為雙線性協(xié)調(diào)位移場(chǎng)；uθ=[uθvθ]T，為由單元節(jié)點(diǎn)角位移分量θi引起的附加位移場(chǎng)。附加位移場(chǎng)為

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

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

式中：N0i、Nuθi、Nvθi為單元對(duì)不同位移場(chǎng)的形函數(shù)。

式中：ξ、η和ξi、ηi分別為單元和單元節(jié)點(diǎn)的自然坐標(biāo)系坐標(biāo)。

GQ12的應(yīng)變場(chǎng)為

單元?jiǎng)偠染仃嚍?/p>

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

1.2 剪應(yīng)變混合插值板元MITC4

基于一階剪切變形理論的Mindlin 板元對(duì)位移場(chǎng)和轉(zhuǎn)角場(chǎng)各自獨(dú)立插值，由于單元邊界只要求C0類連續(xù)，因此滿足網(wǎng)格靈活性的需求，以及舵面等復(fù)雜外形結(jié)構(gòu)的分析需要。一階剪切變形理論放棄了Kirchhoff 理論的直法線假設(shè)，因此對(duì)于橫向剪切變形明顯的中厚板的分析，Mindlin 板元的預(yù)測(cè)效果相比Kirchhoff板元更好。但是，Mindlin板用于薄板時(shí)存在剪切閉鎖和應(yīng)力振蕩等問(wèn)題，不能滿足幾何顯式加筋拓?fù)鋬?yōu)化對(duì)有限元分析的要求。

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

圖1 板單元位移場(chǎng)Fig.1 Displacement field of plate element

由圖1可得，板單元的橫向剪切應(yīng)變?yōu)?/p>

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

Mindlin板理論中，對(duì)于系統(tǒng)總位能公式中的彎曲應(yīng)變能項(xiàng)，MITC4 單元使用與基于位移方法的Mindlin板單元相同的插值方法，廣義位移通過(guò)相同的形函數(shù)獨(dú)立插值為

式中：ωi、θxi、θyi分別為單元各節(jié)點(diǎn)處橫向位移分量、x軸和y軸轉(zhuǎn)角位移分量。

但是，對(duì)于橫向剪切應(yīng)變能項(xiàng)，如果采用相同的插值形式，會(huì)在薄板情況下出現(xiàn)剪切閉鎖問(wèn)題。Bathe 等［14-15］提出的MITC4 單元方案中，構(gòu)造了一種在單元邊和單元內(nèi)滿足剪應(yīng)變一致的插值函數(shù)，插值形式為

式中：γAξz、γBηz、γCξz、γDηz分別為MITC4 單元在ξ-η-z坐標(biāo)系下A、B、C、D四邊中點(diǎn)的橫向剪應(yīng)變。

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

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

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

則有

同時(shí)，由于剪應(yīng)變表達(dá)式定義在ξ-η-z坐標(biāo)系下，需要轉(zhuǎn)化到x-y-z坐標(biāo)系下，轉(zhuǎn)換形式為

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

1.3 驗(yàn)證算例

為驗(yàn)證GQ12 膜元與MITC4 板元耦合平板殼元的分析精度，設(shè)計(jì)一個(gè)四邊形折板，其幾何尺寸如圖3所示。折板一側(cè)約束，另一側(cè)一點(diǎn)受集中力載荷F=-1 000 N 的作用，取結(jié)構(gòu)E=210 GPa，v=0.3，分別對(duì)5 mm厚度及10 mm厚度折板算例求解結(jié)構(gòu)應(yīng)變。

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

5 mm 厚度折板GQ12 膜元與MITC4 板元耦合平板殼元自編Matlab 程序的分析結(jié)果如圖4（a）所示，采用商業(yè)有限元分析軟件Ansys的分析結(jié)果如圖4（b）所示。不同厚度折板加載節(jié)點(diǎn)和最大位移節(jié)點(diǎn)Z軸方向位移數(shù)值對(duì)照見表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 四邊形折板分析結(jié)果對(duì)比Fig.4 Analysis results comparison of the quadrilateral folded-plate

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

2 自適應(yīng)成長(zhǎng)法

在自適應(yīng)成長(zhǎng)法結(jié)構(gòu)加筋拓?fù)鋬?yōu)化中，采用基結(jié)構(gòu)法建立優(yōu)化幾何模型，其優(yōu)點(diǎn)在于避免優(yōu)化迭代過(guò)程中網(wǎng)格的重劃分。參考承載和約束條件確定“種子點(diǎn)”的位置，如圖5（a）所示。從“種子點(diǎn)”生長(zhǎng)出來(lái)的加強(qiáng)筋根據(jù)一定規(guī)則成長(zhǎng)或退化；生長(zhǎng)到一定尺度可以“分歧”，再?gòu)摹胺制琰c(diǎn)”繼續(xù)生長(zhǎng)，在下步迭代中參與成長(zhǎng)與退化；退化到一定尺度則“消失”，如圖5（b）—圖5（c）所示。這個(gè)過(guò)程在優(yōu)化迭代中反復(fù)進(jìn)行，直至滿足收斂條件。

圖5 加強(qiáng)筋拓?fù)鋬?yōu)化的自適應(yīng)成長(zhǎng)過(guò)程Fig.5 Adaptive growth process of stiffener topology optimization

2.1 數(shù)學(xué)模型及靈敏度分析

靜載荷下，薄壁結(jié)構(gòu)自適應(yīng)成長(zhǎng)法加筋拓?fù)鋬?yōu)化問(wèn)題的數(shù)學(xué)模型可以表述為：在給定的體積約束下，選擇加強(qiáng)筋厚度作為設(shè)計(jì)變量，以結(jié)構(gòu)在載荷條件下的總應(yīng)變能為目標(biāo)函數(shù)，逐步尋求最優(yōu)解。其統(tǒng)一的數(shù)學(xué)模型可以表達(dá)為

式中：X為設(shè)計(jì)變量；xi為第i個(gè)活動(dòng)筋板單元的厚度；n為加強(qiáng)筋單元的總數(shù)；Φ(X)為目標(biāo)函數(shù)，即總應(yīng)變能；V、V0分別為整體結(jié)構(gòu)總體積和結(jié)構(gòu)初始體積；η為體積約束因子；xmin、xmax為設(shè)計(jì)變量xi的取值下限和上限；g(X)為約束函數(shù)。

薄壁結(jié)構(gòu)的柔度關(guān)于加強(qiáng)筋厚度的靈敏度計(jì)算公式為

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

2.2 優(yōu)化設(shè)計(jì)流程

基于優(yōu)化數(shù)學(xué)模型和自適應(yīng)成長(zhǎng)法基本原理，以移動(dòng)漸近線法（method of moving asymptotes，MMA）作為迭代更新準(zhǔn)則對(duì)結(jié)構(gòu)分布進(jìn)行尋優(yōu)，具體流程如圖6所示。

圖6 自適應(yīng)成長(zhǎng)法優(yōu)化設(shè)計(jì)流程Fig.6 Flow chart of optimization design of adaptive growth method

1） 建立基結(jié)構(gòu)，設(shè)置設(shè)計(jì)域。

2） 選取“種子點(diǎn)”，設(shè)置優(yōu)化參數(shù)。參考結(jié)構(gòu)的載荷和支撐情況，選取基結(jié)構(gòu)上若干節(jié)點(diǎn)作為“種子點(diǎn)”，給定加強(qiáng)筋的初始厚度x0、分歧臨界值xb、體積約束因子η、收斂容差ε以及最大迭代次數(shù)N。

3） 靈敏度分析。對(duì)結(jié)構(gòu)進(jìn)行有限元分析，通過(guò)解析法計(jì)算活動(dòng)筋板單元的靈敏度。

4） 成長(zhǎng)和退化。判斷收斂狀態(tài)，如果滿足筋板成長(zhǎng)條件，通過(guò)MMA 算法更新xi；如果筋板厚度達(dá)到分歧臨界值xb，該單元進(jìn)行分歧；如果滿足退化條件，該單元退化。

3 舵面結(jié)構(gòu)加筋布局優(yōu)化

3.1 舵面基結(jié)構(gòu)構(gòu)建

根據(jù)舵面幾何形狀，在優(yōu)化設(shè)計(jì)前構(gòu)建拓?fù)鋬?yōu)化基結(jié)構(gòu)，舵面外層為蒙皮結(jié)構(gòu)，定義厚度為0.5 mm。選擇舵面內(nèi)部空間區(qū)域作為設(shè)計(jì)域，定義基結(jié)構(gòu)初始厚度為0.5 mm，如圖7 所示。由于舵面的幾何特征，構(gòu)成基結(jié)構(gòu)的單元為任意凸四邊形單元，GQ12 膜元和MITC4 板元耦合構(gòu)造的四邊形等參平板殼元適用于舵面基結(jié)構(gòu)的有限元分析。

圖7 舵面優(yōu)化初始基結(jié)構(gòu)Fig.7 Initial ground structure of rudder structure for optimization

3.2 基于自適應(yīng)成長(zhǎng)法的加筋布局

根據(jù)舵面的實(shí)際工作狀態(tài)，對(duì)其約束施加于舵軸，受單側(cè)橫向均布載荷。因此，本文以舵面舵軸處約束點(diǎn)為自適應(yīng)成長(zhǎng)法的“種子點(diǎn)”，以整體結(jié)構(gòu)應(yīng)變能最小為設(shè)計(jì)目標(biāo)，骨架結(jié)構(gòu)體積為約束條件，對(duì)舵面進(jìn)行加筋拓?fù)鋬?yōu)化，尋找最優(yōu)的骨架結(jié)構(gòu)分布。優(yōu)化的數(shù)學(xué)模型為

式中：T為設(shè)計(jì)變量，為每個(gè)活動(dòng)筋板的厚度；U(T)為目標(biāo)函數(shù)，即整體結(jié)構(gòu)應(yīng)變能；V、V0分別為加筋舵面結(jié)構(gòu)總體積和初始實(shí)體舵面結(jié)構(gòu)體積；η=20%為體積約束因子。

經(jīng)過(guò)拓?fù)鋬?yōu)化設(shè)計(jì)得到的舵面內(nèi)部最優(yōu)加筋布局如圖8 所示。圖中呈現(xiàn)出清晰的分枝結(jié)構(gòu)，拓?fù)湫螒B(tài)合理且具有較優(yōu)的可制造性。

圖8 拓?fù)鋬?yōu)化后的舵面加筋布局Fig.8 Rudder stiffener layout obtained by topology optimization

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

圖9 舵面結(jié)構(gòu)加筋布局優(yōu)化迭代曲線Fig.9 Iteration history curve of rudder stiffener layout

3.3 設(shè)計(jì)結(jié)果分析

為了評(píng)估加筋拓?fù)鋬?yōu)化效果，利用耦合四邊形平板殼元建立舵面的傳統(tǒng)蒙皮-網(wǎng)格式骨架設(shè)計(jì)模型，如圖10所示。對(duì)文獻(xiàn)［1］中介紹的傳統(tǒng)骨架設(shè)計(jì)進(jìn)行建模并分析，計(jì)算得到拓?fù)鋬?yōu)化前后舵面的整體結(jié)構(gòu)應(yīng)變能，如圖11所示。

圖10 舵面?zhèn)鹘y(tǒng)骨架設(shè)計(jì)Fig.10 Traditional rudder skeleton design

圖11 舵面結(jié)構(gòu)優(yōu)化前后整體結(jié)構(gòu)應(yīng)變Fig.11 Strain of overall structure before and after rudder structure optimization

對(duì)比傳統(tǒng)網(wǎng)格式骨架設(shè)計(jì)與加筋拓?fù)鋬?yōu)化后舵面整體結(jié)構(gòu)的應(yīng)變能數(shù)值以及結(jié)構(gòu)質(zhì)量，結(jié)果見表3。由表3可見，優(yōu)化后結(jié)構(gòu)的應(yīng)變能為17.28 J，較傳統(tǒng)設(shè)計(jì)25.80 J降低了33%，結(jié)構(gòu)剛度得到增強(qiáng)；傳統(tǒng)骨架設(shè)計(jì)的質(zhì)量為1.19 kg，優(yōu)化后結(jié)構(gòu)的質(zhì)量降低了29%。

表3 舵面結(jié)構(gòu)優(yōu)化前后性能比較Tab.3 Performance comparison before and after rudder structure optimization

4 結(jié)束語(yǔ)

本文利用GQ12 膜元與MITC4 板元耦合構(gòu)造了四邊形等參平板殼元，通過(guò)有限元分析實(shí)現(xiàn)了基于自適應(yīng)成長(zhǎng)法的不規(guī)則薄壁結(jié)構(gòu)的加筋設(shè)計(jì)。隨后以輕質(zhì)、高承載的舵面結(jié)構(gòu)為例，獲得了靜載荷狀態(tài)下舵面內(nèi)部的最優(yōu)骨架布局。經(jīng)過(guò)對(duì)優(yōu)化前后舵面結(jié)構(gòu)的性能對(duì)比，驗(yàn)證了GQ12膜元與MITC4板元耦合平板殼元分析工具的準(zhǔn)確性和實(shí)用性，同時(shí)驗(yàn)證了自適應(yīng)成長(zhǎng)法對(duì)舵面優(yōu)化設(shè)計(jì)的有效性。