基于等幾何方法的功能梯度Mindlin板結構拓撲優化設計

劉宏亮，姚力銘，王培金

(1.沈陽航空航天大學遼寧省飛行器復合材料結構分析與仿真重點實驗室，沈陽 110136；2.大連理工大學工業裝備結構分析國家重點實驗室，大連 116023)

0 引言

隨著高新技術的發展，飛行器通常需要在惡劣極端的環境下長時間工作[1]，對結構和材料性能的要求越來越高。由于比剛度和比強度等方面的優點，復合材料板在航空航天工業的應用愈發普遍。與傳統的復合材料板不同，功能梯度復合板具有材料相之間平滑過渡且連續分布的特點，可以在很大程度上避免材料界面的出現對結構力學性能的影響，例如降低應力集中和減少殘余應力。這實際上也是功能梯度結構相比傳統復合材料結構的一個優勢[2-3]。

在航空航天產品設計過程中，利用優化方法獲得輕質高承載結構是實現裝備輕量化和高性能化的一個重要且有效途徑。相比尺寸和形狀優化，結構拓撲優化的設計潛力更大，對結構的綜合性能提升也更加明顯，但同時其挑戰性也更大[4-5]。目前，各種拓撲優化方法不斷發展，例如變密度方法[6]、漸進結構優化方法[7]、水平集方法[8]、移動變形組件/孔洞方法[9-10]等。而且，相較于單相材料拓撲優化，考慮多相材料的結構拓撲優化研究近年來受到了更多的關注[11-13]。然而，多材料結構拓撲優化需要一個合適的優化模型，可以有效描述迭代過程中設計域內不同材料的分布，同時應該方便于設計優化過程的靈敏度分析計算[14-16]。

功能梯度結構的拓撲優化可以是單材料優化，也可以是多材料優化。當功能梯度結構的材料漸變方式確定時，只優化結構的拓撲構型屬于單材料優化，例如，材料沿著厚度或者面內方向漸變且其分布確定的功能梯度板結構拓撲優化[17]。當梯度材料漸變與結構拓撲同步優化以實現功能梯度結構性能最優時，其本質上是一種多材料結構拓撲優化[18-19]。與傳統多材料結構拓撲優化不同的是，這樣的設計在材料相之間具有平滑且連續分布的特點。值得指出，現有的功能梯度結構拓撲優化研究主要考慮平面問題，尚未考慮工程中普遍存在的板結構設計。事實上，板結構的多材料拓撲優化設計研究也是近年來才受到關注[20-21]。因此，本文研究的重點是一種沿著面內方向材料漸變的功能梯度板結構拓撲與材料分布協同設計優化問題。

對于板結構的分析和設計，目前主要采用基于有限元方法的Kirchhoff板單元和Mindlin板單元。其中，Kirchhoff板單元適合薄板模型，而Mindlin板單元適合中厚板模型。薄板模型的適用范圍有限，因此Mindlin板單元在工程中的應用范圍更廣。然而，Mindlin板單元進行有限元分析時可能會存在數值閉鎖問題，嚴重影響數值計算的精度。為了實現更具工程應用性的板結構拓撲優化，克服數值求解難題，本文基于等幾何分析[22]發展一種功能梯度Mindlin板結構的拓撲優化設計方法。采用具有高階連續性的非均勻有理B樣條(non-uniform rational B-spline，NURBS)基函數作為功能梯度板結構建模、分析和設計的基礎，獲得了梯度材料與結構拓撲協同優化設計方案。

1 功能梯度Mindlin板等幾何拓撲優化設計

1.1 Mindlin板結構等幾何分析

本文考慮的Mindlin板基于一階剪切變形理論，假設板中面的法向纖維在變形過程中保持平直，長度不變，但不一定保持垂直于中面。因此，橫向剪切變形不可忽略，Mindlin板的應變能由下式給出

(1)

式中，σf和εf是彎曲應力和應變，σc和εc是橫向剪切應力和應變，α為剪切修正系數。線彈性應力應變關系可定義為

σf=Dfεf

σc=Dcεc

(2)

其中，Df和Dc分別表示彎曲項與剪切項的彈性常數矩陣，即

(3)

式中，E是彈性模量，ν是泊松比，G則是剪切模量。將式(2)和式(3)帶入Mindlin板應變能式中可得

(4)

基于等幾何分析，利用NURBS曲面表示結構的幾何模型，即

(5)

(6)

式中，Ni，p和Mj，q分別表示定義在節點矢量{ξ0，ξ1，…，ξm1}和{η0，η1，…，ηm1}上的B樣條基函數，wi，j是控制點對應的權重系數。

在Mindlin板結構等幾何分析中，NURBS基函數不僅用于結構幾何模型的定義，同時也作為等參元的形函數。離散后的Mindlin板線彈性控制方程可表示為

KU=F

(7)

式中，U為位移矢量，K表示結構的整體剛度陣，它由單元剛度陣Ke組裝而成，表示如下

(8)

式中，Bf和Bc分別為彎曲項和剪切項的應變位移矩陣，h表示板的厚度，ui，j為結構的位移場。

本文利用等幾何分析的k細分方案可以實現NURBS基函數跨單元高階的Cp-1次連續，在很大程度上可以緩解閉鎖作用的影響。對于較薄的板模型，將采用NURBS等幾何分析單元結合縮減積分的策略解決數值閉鎖問題。

為了構建一個基于連續密度的梯度漸變結構優化問題，式(8)中的彈性系數矩陣近似為不同材料相的彈性常數，這將在下面討論的插值方案中實現。

1.2 功能梯度拓撲優化的插值方案

在功能梯度結構中，材料之間不存在明確的界面，材料性能平滑變化。借鑒基于SIMP的變密度方法插值方案，假設板由兩種材料構成，彈性模量分別為E1和E2，通過設計域內某一點上兩種材料所占的百分比，即體積分數來表示這種連續變化材料的性質。

功能梯度材料性質的計算通常采用數值均勻化法和混合法。本文基于混合法，假設兩種材料均勻混合，沒有微觀結構，如圖1所示。目前，被廣泛采用的幾種混合模型包括Voigt模型、Halpin-Tsai復合模型、Voigt-Reuss界限和Hashin-Shtrikman界限[23]。本文研究重點是采用等幾何方法實現優化設計，因此選擇使用各向同性復合材料的Hashin-Shtrikman上界和下界的平均值作為材料模型。

圖1 結構中某一點材料混合物

當兩種材料的泊松比等于1/3時，Hashin-Shtrikam的上限和下限可簡化為楊氏模量條件

(9)

當y→1時，Hashin-Shtrikam模型上下限均為E1，當y→0時，Hashin-Shtrikam模型上下限均為E2。假設E1>E2，當y=1時，此時結構中某一處全為強材料，當y=0時即為弱材料。當結構中一點是兩種材料均勻混合時，采用Hashin-Shtrikam模型上下限均值近似表示為

(10)

借鑒SIMP方法，各向同性均質材料的性質被描述為E(x)=xpE0，設計變量為0時表示此處為優化所產生的孔洞。為了實現材料屬性和結構拓撲的協同設計，功能梯度材料的插值模型被定義為

(11)

式中，p為懲罰系數。值得指出，對于式(11)中的設計變量y并不需要懲罰，因為設計目標是獲得光滑漸變的功能梯度結構。為了避免數值奇異，插值格式中采用一個足夠小的量綱為1的數Emin，即

(12)

1.3 功能梯度結構拓撲與材料協同優化列式

考慮體積約束下的功能梯度Mindlin板結構最小柔順度優化問題，當功能梯度材料是由兩種材料(例如金屬與陶瓷)構成時，其對應的拓撲優化列式可以表示為

minimize:J(xi,yi)=FTU(xi,yi)

(13)

式中，J表示目標函數即結構的柔順度，xi和yi表示控制點處的設計變量，分別對應總材料的體分比和強材料的體分比，Ve表示NURBS單元體積，V0表示結構總體積，fvx和fvy分別為材料用量的限制分數，即體積分數。fvx限制結構的總材料用量，fvy則控制有材料處的強材料占比。ym和yM為給定設計變量yi的上下限，由式(9)可知，若ym和yM取0和1則對應為Hashin-Shtrikam模型的上下限，為了保證實現強弱材料的光滑漸變，本文中ym和yM分別設置為0.1和0.9。

對于本文基于等幾何方法的功能梯度Mindlin板結構拓撲優化，用于幾何建模和結構分析的NURBS基函數也用于表示相對密度場，即

(14)

其中，Θ表示設計變量x和y組成的集合。由于NURBS基函數的非負性和規范性，當設計變量在給定的設計區間變化時，所對應的相對密度分布也會滿足規定的區間分布，從而建立起合理的控制點密度設計變量和材料相對密度分布之間的映射關系。此外，NURBS基函數的高階連續性可以避免棋盤格模式等數值問題[24]。

2 設計優化靈敏度分析

本文采用基于梯度的移動漸進線法[25](MMA)來更新設計變量，進而求解建立的優化列式(13)。這里需要計算目標函數和約束函數的靈敏度。功能梯度板的幾何建模、結構分析和材料分布場采用統一的參數化，NURBS基函數的特性保證了解析靈敏度的有效計算。目標函數對于設計變量xe的靈敏度可表示為

(15)

式中單元剛度陣對設計變量xe的導數表達式為

(16)

根據鏈式求導法則可得到

(17)

式(17)中Mindlin板單元剛度陣由兩項組成，分別為彎曲項與剪切項，則?E/?xe的計算表達式為

(18)

由式(15)和式(18)可得

(19)

目標函數對于設計變量ye的靈敏度表達式為

(20)

式中，?Ke/?ye可以參考式(16)和式(17)，而彈性模量對設計變量ye的一階導數為

(21)

式中，E′(y)表示Hashin-Shtrikam界限對設計變量ye求導，可以表示為

(22)

根據式(21)和式(22)可得

(23)

另外，體積約束的靈敏度可以直接求導得出。

3 數值算例與討論

通過兩個數值算例證明本文方法可以有效獲得體積約束下的功能梯度Mindlin板結構拓撲優化設計。優化算例中采用二次NURBS基函數，為避免數值閉鎖現象對Mindlin板結構等幾何分析精度的影響，采用縮減的兩點積分策略。剪切修正系數設置為α=5/6，材料用量限制比值為0.5。當連續兩次迭代的目標值相對變化小于0.01或達到最大迭代次數50時，優化過程終止。算例中材料、幾何以及載荷參數選擇為量綱為1的參數。

3.1 四邊固定的方板

第一個算例考慮一個四邊固定的矩形方板(a×a=1×1)，在板中心處施加一個集中載荷(F=1)，如圖2所示。板的厚度為0.1，參數設置參考了文獻[20]，用紅色表示強材料E1=4，藍色表示弱材料E2=2，兩種材料泊松比為0.3。總體積分數設為0.5，強材料體積分數占比V1從0.4下降到0.1，弱材料體積分數占比V2從0.1上升到0.4，J表示目標函數即柔順度。將設計域細分為100×100的單元時，優化結果如圖3所示。

(a)V1=0.4，V2=0.1，J=30.52

在相同厚度和邊界條件下，可以觀察到材料填充面積與體積分數具有明顯的相關性，隨著硬質材料(紅色)的體積分數逐漸下降，結構的柔順度逐漸增加。可以推斷出，在對結構剛度貢獻較小的地方，軟質材料傾向于替代硬質材料，而強材料總是分配在中心區域(集中載荷作用區)和應力集中區。4種不同強材料體積分數所對應的優化迭代曲線如圖4所示。

圖4 4種不同強材料體積分數下的迭代曲線

圖3中的4種優化結果，拓撲設計構型幾乎一致，不同點在于材料分布，即強弱材料之間的過渡區域。圖3右側顏色欄中的0.1和0.9對應式(13)設定的Hashin-Shtrikam模型上下限。本文協同優化體現在兩種設計變量分別控制結構拓撲與材料配比，如圖5所示。右側“0”表示結構中的孔洞，即圖5黑白設計中的白材料。由于圖3中可視化點的數量有限，功能梯度材料的漸變細節表示得并不是很清楚。由圖4可知，其收斂速度較快，在20多步就已經收斂。在不改變設計優化計算效率的情況下，為了更清楚地顯示梯度漸變細節，在圖3優化結果的基礎上，將用于描述最終優化結果的可視化點擴大100倍，如圖6所示，可見優化結果具有清晰的拓撲設計和平滑的材料過渡區域。

圖5 體積約束下的功能梯度板協同優化設計

(a)強材料體積分數為0.4

3.2 四邊簡支的方板

第2個算例考慮一個四邊簡支的矩形板，厚度為0.01，設計域參考圖2，板的中心仍然施加一個集中載荷。強材料(紅色相)的彈性模量E1=2，弱材料的彈性模量E2=1，兩種材料泊松比為0.3，總體積分數為0.5，其中強材料占比為0.3。該算例的設計參考了文獻[20]在多材料薄板拓撲優化設計中的算例。將設計域離散為50×50的單元，最終的優化結果如圖7所示，目標函數的迭代曲線如圖8所示。

(a)結構拓撲

圖8 四邊簡支功能梯度板的迭代曲線

功能梯度結構拓撲優化設計也可以視作一種多材料結構拓撲優化設計，但是傳統多材料結構拓撲設計的不同材料相之間存在明顯的界面。圖7(a)中的優化設計構型與參考文獻[20]中單材料結構設計類似，但通過功能梯度材料結構協同設計優化，最終結構的目標函數值為137 319.91，在材料參數設置一致的情況下，結構性能優于參考文獻中的雙材料板拓撲優化設計，目標函數下降了11.28%，證明了本文設計方案的有效性和優越性。另外，可視化點的數量不依賴建模、分析和設計變量的數量，可以在保證功能梯度材料相之間連續且平滑漸變的同時兼顧設計優化的效率。

4 結束語

本文研究了一種在等幾何分析框架下的功能梯度Mindlin板結構拓撲優化方法實現材料漸變分布與結構拓撲協同設計。采用基于NURBS的建模、分析和設計參數化方法來實現梯度材料結構的一體化優化設計。等幾何方法中NURBS基函數用于梯度材料分布拓撲優化使得設計結果可以自然避免拓撲優化的常見數值不穩定問題，高階連續性可以促進材料分布的光滑性以及設計優化的靈敏度分析，因此建立了合理的設計變量與材料分布的關系。數值算例表明，本文方法具有較好的收斂性，在不需要顯著增加計算單元和設計參數的情況下能夠有效獲得具有清晰細節特征的結構拓撲與材料漸變的協同優化設計。通過與已有研究結果對比，本文優化設計方法展現了有效性和優越性。