基于等幾何分析的參數多孔結構拓撲優化*

胡傳豐, 任靖雯, 胡 慧, 藺宏偉,2

(1. 浙江大學 數學科學學院, 杭州 310027;2. 浙江大學 計算機輔助設計與圖形學國家重點實驗室, 杭州 310058)

多孔結構是一種由大量孔洞組成的實體結構, 在自然界和人工制品中廣泛存在, 如木材、骨骼、 珊瑚、 海綿等, 可長期承受較大的靜態載荷和周期載荷. 與傳統結構相比, 多孔結構具有質量輕、 比表面積大、 高滲透性、 高比強度等優點, 以及抗沖擊性[1]、 阻尼增強[2]、 缺陷容忍性[3]等特性. 這些優良特性使其應用范圍遠超出單一功能材料, 因而廣泛應用于組織工程、 輕量化設計及能量吸收等領域. 在組織工程領域, 高滲透性和高比表面積的多孔結構, 有助于建立一個適宜細胞附著、 遷移繁殖、 營養運輸和新陳代謝等的生物微環境, 常被作為組織支架移植到人體組織缺損部位, 輔助組織修復再生[4-6]. 在輕量化設計領域, 由于多孔結構質量輕、 相對密度低, 因此利用多孔結構進行飛機機翼內部結構設計可有效降低機翼重量, 同時提高機翼抗彎剛度[7]. 在能量吸收領域, 由于多孔結構具有較大的壓縮應變, 因此在受到外力沖擊時, 可借助自身結構特性將動能轉變為壓縮能, 從而提高能量吸收能力[8]. 多孔結構的力學性能不僅與材料相關, 還與自身分布相關, 所以需對多孔結構進行結構分析和優化設計, 以提高其力學相關性能.

多孔結構的設計方法主要包括CAD(management software computer aided design)造型設計方法和隱式曲面造型設計方法, 其中CAD造型設計方法適用于設計簡單規則的多孔結構. Cheah等[9]通過對多面體形狀的研究, 設計了基于多面體的多孔單元庫; Lal等[10]提出利用微球填充方法設計多孔支架. 隱式曲面造型以三向周期極小曲面(triply periodic minimal surface, TPMS)為研究熱點. Yoo首先利用TPMS設計了多孔單元庫, 同時提出了利用六面體單元映射的方法構建多孔結構[11]; 之后, Yoo通過對TPMS與實體進行求交并構建多孔結構, 求交運算中引入了距離場算法替代Boole操作, 極大減少了時間的消耗[4]，并利用徑向基函數進行空間插值控制孔徑大小分布, 構建了非均質多孔結構[5]. 為在多孔結構設計中充分利用多孔單元庫, Yang等[12]利用Sigmoid函數和Gauss徑向基函數以任意形狀的過度邊界融合兩種不同類型的多孔單元, 生成了形狀更復雜的多孔結構； Shi等[13]結合TPMS和Sigmoid函數從CT數據中重建多孔支架結構； Feng等[6]利用T樣條函數表示幾何模型, 通過分析TPMS的相關參數與多孔結構的孔隙率、 比表面積之間的關系, 設計了孔隙率、 比表面積可控的多孔結構； Savio等[14]基于TPMS提出了CAD環境下變厚度多孔結構幾何建模的方法.

拓撲結構優化設計以力學原理和數學規劃算法為基礎, 通過優化方法改變工程結構的尺寸、 形狀和拓撲, 在給定的設計域和約束條件下, 實現結構的最佳性能設計. 目前的拓撲優化方法主要分為4類: 變密度法[15-16]、 水平集法[17-18]、 拓撲導數法[19]和相場法[20]. 相比于傳統有限元分析方法, 等幾何分析方法緊密結合幾何模型信息, 避免網格劃分過程, 具有高階連續性, 在保證幾何精確性的同時, 可有效降低求解問題的自由度, 提高計算模擬精度和效率[21-22]. 由于拓撲優化中的變密度法具有直觀的數學模型, 且實現簡單、 計算高效, 因此, 可將等幾何分析與變密度法融合發展形成基于變密度的等幾何拓撲優化方法. Hassani等[23]提出了結合優化準則法的等幾何拓撲優化方法, 并通過二維平面優化問題算例表明該方法可有效抑制棋盤格現象； Qian[24]提出了一種基于B樣條函數的變密度框架下拓撲優化方法, 將密度分布引入B樣條函數空間, 并將控制點對應的相對密度值作為優化變量進行拓撲優化； Liu等[25]利用變密度框架下的等幾何拓撲優化方法分析了全局應力約束下的拓撲優化問題.

目前, 針對多孔單元的拓撲優化研究已取得許多成果[26-29], 對多孔結構的拓撲優化, 主要包含以下兩類方法： 第一類方法[30-32]先優化設計域的材料密度分布, 然后根據密度分布將分析單元替換為對應材料密度下的多孔單元； 第二類方法[33-34]預先設計晶格和單元結構, 然后在一個單元內優化結構尺寸或壁厚. 這兩類方法各有其優缺點, 例如： 第一類方法在結構層次上, 而不是在每個多孔單元中優化拓撲結構; 第二類方法在結構層和單元層上都得到了優化, 但僅適用于具有一定壁厚的規則多孔單元; 上述方法均為基于有限元分析的優化方法, 從而導致在處理相互連通的多孔單元時, 無法保證多孔單元密度分布或壁厚分布的連續性, 進而相鄰多孔單元不能光滑連接. 特別地, Li等[26]基于TPMS提出了功能梯度周期曲面, 構建多孔單元結構與密度建立映射關系, 利用有限元分析方法求解柔順度和熱傳導問題, 尋找最優的單元密度分布, 最后對單元密度進行空間插值獲得節點對應密度值, 并生成功能梯度多孔結構. 本文提出的變密度框架下多孔結構等幾何拓撲優化方法可有效結合上述兩類方法的優勢, 通過TPMS多孔單元特征參數改變單元壁厚, 同時建立與材料相對密度間的映射關系, 以多孔單元特征參數作為優化變量進行拓撲優化, 保證優化構建的多孔結構中多孔單元光滑連接.

本文利用TPMS設計參數多孔單元, 并分析多孔單元特征參數與材料分布間的關系, 再利用變密度框架下的等幾何拓撲優化模型對多孔單元特征參數分布進行結構分析與優化, 以實現最佳的多孔單元特征參數分布, 提高多孔結構的力學性能. 相比基于傳統有限元分析的結構優化, 等幾何分析的物理模型和多孔單元特征參數采用B樣條表示, 提高了仿真計算的精度, 同時多孔結構中多孔單元個數可調控, 無需與分析單元數保持一致, 且多孔結構中多孔單元間連接的連續性得到了保證. 本文方法主要有以下創新點： 1) 提出了基于等幾何分析的參數多孔結構拓撲優化方法; 2) 多孔結構中多孔單元個數可調控; 3) 多孔結構中多孔單元間光滑連接.

1 預備知識

1.1 B樣條

p次B樣條曲線定義[35]為

(1)

其中{Pi}是控制點，Ni,p(u)是定義在節點向量U=(u0,u1,…,um+p+1)上的p次B樣條基函數, 且滿足u0≤u1≤…≤um+p+1. B樣條基函數為

(2)

由此延伸定義出B樣條曲面和B樣條體, 分別為

(3)

(4)

特別地, 由于本文研究薄板多孔結構優化問題, 實際薄板以B樣條曲面形式表示, 但針對平面應力問題, 薄板以二維B樣條曲面表示.

1.2 三向周期極小曲面

TPMS是在歐氏空間中沿3個獨立方向周期性無限延伸的隱式曲面, 具有平均曲率為零的特點, 并將空間平滑而連續地一分為二, 產生連通性優異的孔結構, 是多孔結構設計領域中一種較好的設計工具. 由于TPMS參數表達形式相對復雜, 因此通常采用Fourier級數定義的周期曲面對其逼近[36]:

(5)

其中Ak為振幅,hk為倒空間的格矢量,r為空間位置矢量,λk為波長,Pk為相位,C為等值面閾值常數. 表1列出了P,D,G,IWP 4種類型TPMS的表達式, 其中閾值C的有效范圍確保TPMS連通.αu,αv,αw表示曲面在空間3個方向上的周期. 特別地, 相比于其他類型的TPMS, 這4種曲面具有更大的表面積, 在多孔結構設計中應用廣泛[37]. 此外, 本文利用移動四面體方法[38]提取TPMS, 并構建多孔結構, 如圖1所示, 均為一個完整周期內的曲面.

表1 TPMS三角函數表達式Table 1 Trigonometric function expression of TPMS

圖1 4類TPMSFig.1 Four types of TPMS

1.3 多孔單元

本文基于一個完整周期內4種類型TPMS(圖1)設計了如圖2所示的4種類型多孔單元.

圖2 4類TPMS多孔單元Fig.2 Four types of TPMS porous elements

在多孔結構拓撲優化設計中, 定義C值為多孔單元的特征參數, 用于調控多孔單元孔徑的大小以及多孔單元內材料的相對密度. 通過數據統計與分析, 發現4種多孔單元的特征參數C與相對密度ρ間存在如圖3所示的關系. 數據擬合結果表明, 在一定誤差范圍內, 特征參數C與材料相對密度ρ存在如下關系：

圖3 多孔單元特征參數C與材料相對密度ρ的關系Fig.3 Relationship between characteristic parameter C of porous element and relative density ρ of material

ρ=k1C+k2.

(6)

表2列出了不同多孔單元對應的系數k1和k2.

表2 多孔單元特征參數C與材料相對密度ρ的函數系數Table 2 Function coefficients of characteristic parameter C of porous element and relative density ρ of material

2 基于等幾何分析的設計參數化

等幾何分析是一種新型的數值計算方法, 其將工程結構的幾何表示、 分析及拓撲優化過程緊密結合. 相比于傳統有限元分析方法, 等幾何分析方法具有幾何精確性和單元間高階連續性的優點, 提高了結構分析的精確性和可信性. 結合等幾何分析的拓撲優化方法, 可有效提高結構分析和拓撲優化的效率.

2.1 等幾何結構分析

等幾何分析以等參思想為基礎, 利用幾何模型的樣條基函數和控制點替代有限元分析中的單元形函數和節點. 圖4為雙線性(p=q=1)和雙二次(p=q=2)等幾何單元組成的等幾何分析網格.

圖4 等幾何分析網格Fig.4 Iso-geometric analysis grid

進行結構分析時, 位移場表示為

(7)

其中,Ni,p(u),Nj,q(v)為式(2)中的B樣條基函數，Uij為控制點處的位移系數. 本文利用等幾何方法求解平面應力問題, 每個控制點處對應兩個方向上的位移量Ux,Uy.

線彈性連續體的靜力平衡方程可表示為

KU=F,

(8)

其中K為整體剛度矩陣，U為控制點處位移向量，F為外部載荷向量. 整體剛度矩陣K可由單元剛度矩陣裝配得到, 表示為

(9)

其中Ωd為設計域Ω的離散域. 單元剛度矩陣為

(10)

式中的下標i1,i2表示影響該單元Ωe的控制點索引,λ,μ為Lamé常數, 與材料屬性相關, 計算公式為

E為材料彈性模量,ν為材料Poisson比.

2.2 離散多孔結構特征參數分布

本文采用B樣條函數建立多孔單元特征參數分布, 在二維設計域內任意一點的多孔單元特征參數值可由下式計算:

(11)

其中Cij表示控制點處的多孔單元特征參數變量. 為簡化計算過程, 本文以單元內恒定的特征參數進行結構分析. 以等幾何單元中心點的特征參數值作為該單元的特征參數值, 表示為

(12)

(13)

基于B樣條的多孔單元特征參數分布是光滑連續分布函數, 這種表示方法相當于對特征參數分布施加了光滑效果, 可有效避免在有限元分析中以單元進行拓撲優化出現的棋盤格或孤島現象, 這是基于B樣條參數化設計的一個顯著優勢.

3 多孔結構的拓撲優化

3.1 最小柔度優化理論

結構拓撲優化最早是設計輕量化結構[39], 其中固體各向同性材料懲罰模型(solid isotropic material with penalization, SIMP)[40]應用最廣泛, 其通過在每個材料單元引入相對密度, 再優化材料密度的分布, 所以也稱為變密度法. 這種優化算法對目前多孔結構的優化非常有用, 通過建立多孔單元特征參數分布與材料相對密度間的映射關系, 然后在變密度法框架下基于等幾何分析求解優化問題, 實現最佳性能的多孔單元特征參數分布, 提高多孔結構的性能.

3.2 優化函數

下面介紹多孔結構優化設計中的拓撲優化算法, 結構的總勢能[39]表示為

(14)

本文基于最小柔順度優化問題, 確定多孔單元特征參數在設計結構中的分布, 使得柔順度最小化, 即最小化總勢能, 以使柔順度達到最小. 優化問題的數學表達式為

(15)

其中:k1,k2為多孔單元特征參數與材料相對密度函數系數(見表2);Cij為多孔單元特征參數分布控制點, 即函數優化變量；Ce為等幾何單元中心特征參數(式(12))；t為密度懲罰因子, 本文取值為3； 第一個約束條件KU=F為彈性平衡方程, 表示位移向量U在外部載荷向量F下滿足平衡方程條件； 第二個約束條件為材料體積約束,V(·)表示多孔結構相對體積函數,Ve表示等幾何單元體積,V0表示設計域總體積,τ為給定體積分數, 表示優化后的多孔結構材料體積必須滿足給定材料體積量； 第三個約束條件為特征參數的約束, 表示特征參數約束在有效范圍內, 確?？蓸嫿ㄍ暾亩嗫讍卧?

3.3 靈敏度分析

下面推導在變密度框架下等幾何拓撲優化的靈敏度計算公式, 優化變量為樣條控制點處對應的多孔單元特征參數變量. 首先, 柔順度關于優化變量的偏導數計算公式為

(16)

結合式(13)可得

(17)

其次, 材料體積關于優化變量的偏導數計算公式為

(18)

本文采用優化準則法求解優化問題(15), 參考文獻[39], 利用啟發式迭代更新優化變量:

(19)

其中:m是移動步長, 本文取值0.2;η(0<η<1)是阻尼系數, 本文取值0.5;Bij表達式為

式中γ是Lagrange乘子, 可利用二分法求解.

3.4 多孔結構生成

通過上述結構優化后, 實現了在設計域中的最佳多孔單元特征參數分布, 該多孔單元特征參數分布同樣可定義在參數域中. 由于本文研究平面應力問題, 因此在求解結構優化時設計域及其對應的參數域均為二維, 而實際薄板及多孔結構均為三維. 由于平面應力分析不受第三個維度的影響, 因此可直接將多孔單元特征參數分布和幾何模型控制點映射到第三個維度上, 采用三變量B樣條函數進行表示：

(20)

此時多孔單元特征參數分布映射三元函數C(u,v,w)定義在幾何體參數域內, 則隱式曲面的函數表達式修改為

ψ(u,v,w;αu,αv,αw)=C(u,v,w),

(21)

其中(αu,αv,αw)表示周期參數, 即多孔單元在對應方向的個數. 本文固定αw=1, 使得薄板厚度對應一個多孔單元; 再利用移動四面體方法提取等值面, 構建參數域多孔結構; 最后利用薄板B樣條體函數(式(4))將參數域空間中的多孔結構映射到物理域中, 即得到經過等幾何拓撲優化后的薄板多孔結構. 特別地, 在參數域中生成多孔結構再投影回物理空間中, 可保證所有的多孔單元都是完整的, 這是由于參數域是規則的, 而實際幾何模型并非規則模型.

在懸臂梁多孔結構的拓撲優化設計中, 多孔單元為G型單元. 首先設計優化問題, 并通過拓撲結構優化求解獲得最佳多孔單元特征參數分布, 再將該多孔單元特征參數分布投影回二維參數域中; 然后將特征參數分布拓展至三維獲得三維的多孔單元特征參數分布, 再基于式(21)利用移動四面體方法提取等值面, 構建參數域多孔結構； 最后經過三元B樣條函數的映射將參數域多孔結構映射到物理域中, 實現一定體積分數下最佳性能薄板多孔結構的構建. 多孔結構優化設計算法流程如圖5所示.

圖5 多孔結構優化設計算法流程Fig.5 Workflow of porous structure optimization design algorithm

4 實驗結果與分析

下面給出一些實驗案例, 并進行分析. 為簡化計算過程, 實驗部分數據均采用無量綱量, 包括單位集中載荷為1, 彈性模量E=1, Poisson比ν=0.3. 實驗相關數據列于表3.

4.1 懸臂梁

圖5已給出了基于懸臂梁優化算法的流程, 其中懸臂梁左側固定, 在右下角施加單位集中載荷, 在體積分數為0.4的條件約束下, 基于G型多孔單元進行結構優化. 下面通過與傳統有限元分析方法進行對比, 分析基于等幾何分析進行多孔結構拓撲結構優化的優勢. 特別地, 進行有限元分析的單元個數與等幾何分析方法保持一致. 分別統計基于有限元方法和等幾何方法求解拓撲優化問題的耗時, 結果表明, 利用等幾何分析方法可有加速優化問題求解.

4.1.1 棋盤格現象

本文利用有限元方法對參數多孔結構進行拓撲優化, 其中每個分析單元對應一個特征參數變量, 作為優化變量進行拓撲優化, 最后得到一個在離散設計域的多孔單元特征參數分布, 如圖6(A)所示. 對比如圖6(B)所示的用等幾何分析方法的懸臂梁多孔結構特征參數分布, 利用有限元分析方法進行參數多孔結構的拓撲優化會產生嚴重的棋盤現象, 而等幾何分析方法可有效抑制發生棋盤現象.

圖6 多孔單元特征參數分布對比Fig.6 Comparison of porous element characteristic parameter distribution

4.1.2 多孔單元個數

傳統有限元分析利用離散網格近似幾何模型, 對離散網格單獨進行優化, 多孔單元個數需與分析單元個數保持一致, 如圖7(A)所示. 而基于等幾何分析方法, 通過對等幾何單元對控制點進行優化, 最后得到以B樣條函數表示的多孔單元特征參數分布, 在構建多孔結構時多孔單元個數無需與等幾何單元個數保持一致, 可通過修改式(21)中周期參數αu,αv調控多孔單元個數. 圖7(B)和圖7(C)分別為基于等幾何分析方法進行拓撲優化后構建的兩種不同多孔單元個數的懸臂梁多孔結構.

圖7 多孔單元個數對比Fig.7 Comparison of number of porous elements

4.1.3 多孔單元連接

有限元方法以單元為分析目標, 在優化得到多孔單元特征參數分布后, 分別以分析單元的特征參數值構建多孔單元, 最后通過單元映射的方法將多孔單元映射到分析單元中生成多孔結構, 由于多孔單元特征參數分布存在嚴重的棋盤現象, 因此兩兩相鄰多孔單元間無法保證光滑連續. 如圖8(A)所示, 多孔結構中相鄰的多孔單元間不能光滑連接, 通常需要對多孔結構進行光滑處理, 這會改變拓撲優化的結果.

而基于等幾何分析的方法, 以多孔單元特征參數變量控制點為優化目標, 最終獲得一個光滑分布函數, 再將該函數代入式(21), 可構建一個光滑連續的等值面, 基于該光滑連續等值面構建的多孔結構, 可保證多孔單元間光滑連續, 如圖8(B)所示.

圖8 多孔單元間連接性對比Fig.8 Comparison of connectivity between porous elements

4.2 MBB

圖9為平面應力狀態下MBB(miniature bending beam)的幾何模型, 其中圖9(A)的MBB左下角受滾動鉸鏈約束, 右下角固定支撐, 梁上部中點處受單位集中載荷, 在體積分數為0.45的約束條件下, 基于P型多孔單元進行拓撲結構優化, 優化后多孔單元特征參數分布在設計域中, 如圖9(B)所示, 最后生成的多孔結構如圖9(C)所示.

圖9 MBB多孔結構拓撲優化(P型多孔單元)Fig.9 Topology optimization for porous structure of MBB (P type porous element)

4.3 1/4圓環

圖10為平面應力狀態下1/4圓環結構的幾何模型, 其中圖10(A)的1/4圓環模型底邊固定, 在左上角處受單位集中載荷, 在體積分數為0.45的約束條件下, 基于D型多孔單元進行拓撲結構優化, 優化后多孔單元特征參數分布在設計域中, 如圖10(B)所示, 最后生成的多孔結構如圖10(C)所示.

圖10 1/4圓環多孔結構拓撲優化(D型多孔單元)Fig.10 Topology optimization for porous structure of quarter annulus (D type porous element)

4.4 多載荷懸臂梁

對多載荷的情形本文也進行了實驗分析, 圖11為多載荷下懸臂梁的幾何模型, 其中圖11(A)的懸臂梁左側固定, 在右側上下兩端點處分別受反向的單位集中載荷, 在體積分數為0.4的約束條件下, 基于IWP型多孔單元進行拓撲結構優化, 優化后多孔單元特征參數分布在設計域中, 如圖11(B)所示, 最后生成的多孔結構如圖11(C)所示.

圖11 多載荷懸臂梁多孔結構拓撲優化(IWP型多孔單元)Fig.11 Topology optimization for porous structure of cantilever beam under multiple loads (IWP type porous element)

綜上所述, 本文提出了一種基于變密度框架下的參數多孔結構等幾何拓撲優化方法. 首先, 基于TPMS設計多孔單元, 并分析多孔單元特征參數與材料相對密度分布間的關系; 其次, 考慮薄板多孔結構平面應力的靜態平衡問題, 以多孔單元特征參數為優化對象, 利用變密度框架下的等幾何分析方法在材料均勻分布的設計空間實現最佳的特征參數分布; 最后, 利用優化后的多孔單元特征參數分布構建多孔結構, 實現了同等材料體積下力學性能最佳的多孔結構設計. 由于等幾何單元間具有高階連續性, 因此保證了應力函數的連續性, 提高了計算精度與效率. 此外, 多孔單元特征參數分布用樣條形式表示, 方便調控多孔結構中的多孔單元數, 且保證了多孔單元間的光滑連接.