隱式曲面梯度多孔結構拓撲優化設計方法

孫鵬飛,張躍,尹鵬,劉宏磊,李寶童

(1.西安交通大學機械工程學院,710049,西安;2.西安交通大學現代設計及轉子軸承系統教育部重點實驗室,710049,西安)

多孔結構是結構/功能一體化的優良載體,具有低密度、高比表面積[1]、高比力學性能[2-3]及優良的吸能特性等特點,在航空航天、汽車和醫學等領域具有廣泛的應用前景。然而,多孔結構的功能特性與其幾何構型存在復雜的耦合關系,導致多孔結構設計的復雜度急劇上升。因此,需進一步研究多孔結構的優化設計方法,實現對其功能特性的調控,以滿足復雜工程應用的需求。

多孔結構的功能特性取決于其多孔單胞構型與宏觀材料分布形式。隨著結構設計方法的快速發展,例如計算機輔助設計法[4]、圖像法[5]、隱式函數法[6]和拓撲優化法[7-8]等,豐富了具有良好力學性能的多孔單胞構型。相比之下,由隱式函數法設計的多孔單胞具有參數化、設計便捷、可設計性強等優點。作為典型的隱式參數化模型,具有零平均曲率的極小曲面引起了相關領域的廣泛關注[9-10]。根據極小曲面在空間周期延伸的維度,可將其分為單周期極小曲面、雙周期極小曲面和三周期極小曲面(TPMS)。其中,三周期極小曲面廣泛存在于自然界中,例如蝴蝶翅膀、甲蟲外骨骼等[11],因其幾何構型呈現出獨特的對稱性,具備高比強度、軸對稱剛度、孔洞連通性和良好的吸能特性等優點[12-13]。然而,現有極小曲面的研究側重于構型的拓撲以揭露其物理特性[14],難以充分發揮其性能優勢。此外,梯度多孔結構作為一種材料梯度分布的多孔結構,其功能特性呈現出漸進性和局部性的變化。相比于均勻多孔結構,梯度多孔結構在提升結構剛度、抗屈曲能力[15]和吸能特性[16]等方面有明顯的優勢。因此,為確定梯度多孔結構的最優材料分布,需進一步研究多孔結構的優化設計方法。

近年來,采用拓撲優化方法對梯度多孔結構進行優化設計已成為一種趨勢[17-18]。在本質上,拓撲優化和梯度多孔結構都考慮材料屬性的連續性變化[19-20],且在已知載荷和邊界條件下,拓撲優化能夠確定梯度多孔結構在空間中最優材料分布形式。為實現多孔單胞的宏觀梯度分布,Brackett等將多孔結構的體積分數映射到未懲罰的固體各向同性材料懲罰模型法的中間密度上,得到了梯度多孔結構[21]。在此基礎上,還可以桁架單元、六邊形蜂窩單元等為代表性體積單胞[22],設計非均勻壁厚的梯度多孔結構[23]。此外,張衛紅等將均勻化理論和拓撲優化方法相結合提出一種以宏觀結構性能為目標、材料表征體胞構型為變量的梯度多孔結構優化設計方法,實現了多孔單胞構型與宏觀材料分布的并行設計[24]。然而,采用拓撲優化對梯度多孔結構進行優化設計時,多孔單胞構型的優化過程復雜、多孔單胞間的連續性難以保障。因此,有必要對多孔單胞的設計方法、多孔結構的連續性展開研究,豐富多孔結構的多樣性,釋放其工程應用潛力。

為實現多孔單胞構型對多孔結構功能特性的調控,本文提出了隱式曲面梯度多孔結構拓撲優化設計方法。結合數值均勻化法[25]建立隱式曲面梯度多孔結構柔度最小化拓撲優化模型,得到了幾何和功能呈梯度分布的隱式曲面多孔結構,并通過數值案例和三點彎曲實驗驗證了本文所提方法的可行性與有效性。

1 極小曲面多孔結構幾何建模

1.1 三周期極小曲面幾何描述

三周期極小曲面是一種由隱式水平集函數定義的曲面式結構,具有較強的可設計性。因此,本文采用極小曲面作為梯度多孔結構的代表性體積單元。為控制極小曲面的體積分數,在其隱式水平集函數中引入水平參數t,如下式

ΦP=C(X)+C(Y)+C(Z)-t

(1)

ΦG=C(X)S(Y)+C(Y)S(Z)+C(Z)S(X)-t

(2)

ΦI-WP=-C(X)C(Y)-C(Y)C(Z)-C(Z)C(X)-t

(3)

式中:C(·)為余弦函數;S(·)為正弦函數;X=2πx/L,Y=2πy/L,Z=2πz/L,x、y、z為高維物理空間坐標;L為單胞尺寸。式(1)～(3)分別為Primitive(P)型曲面、Gyroid(G)型曲面和I-Wrapped Package(I-WP)型曲面的四維隱式水平集函數,根據水平集函數定義,可得到極小曲面的實體區域。

(4)

式中:Φ是極小曲面隱式水平集函數;Ω是曲面實體區域;?Ω是曲面邊界;D是包含實體區域和曲面邊界的空間。極小曲面多孔結構隱式建模過程如圖1所示。

圖1 極小曲面多孔結構隱式建模過程Fig.1 The implicit modeling of TPMS

在極小曲面四維隱式水平集函數中,t代表水平集函數的一個水平面,通過改變水平參數t可以改變水平面到體心的距離,實現實體區域大小變化,從而控制多孔結構的體積分數。當t極大或極小時,極小曲面會出現斷裂現象,導致結構在歐式空間中不連續,極小曲面多孔結構如圖2所示。

(a)P型曲面(b)G型曲面 (c)I-WP型曲面圖2 極小曲面多孔結構Fig.2 TPMS-based cellular structures

在極小曲面隱式水平集函數中引入罰函數,獲得如圖3所示的極小曲面骨架結構,保證極小曲面多孔結構在小體積分數下具有良好的連續性,表達式如下

(5)

(6)

C(Z)C(X))-(C(2X)+C(2Y)+C(2Z))-t

(7)

(a)P型曲面(b)G型曲面 (c)I-WP型曲面圖3 極小曲面骨架結構Fig.3 TPMS-based skeleton structures

1.2 多孔結構混合參數化建模

拓撲優化方法設計的結構幾何構型復雜,構造幾何參數驅動模型梯度漸變是梯度多孔結構優化設計的基礎。圖4給出了三周期極小曲面參數化模型。為實現極小曲面與拓撲優化的結合,對極小曲面隱式水平集函數線性加權,構造了混合水平集函數

(8)

圖4 三周期極小曲面參數化模型Fig.4 The parameterized TPMS models

根據參與建模的極小曲面類型和構造方法,將所設計的極小曲面多孔結構分為實心多孔結構、空心多孔結構和混合多孔結構。本文采用I-WP型和P型兩種極小曲面進行混合建模。

為建立空心多孔結構參數化模型,對式(8)與同類型極小曲面隱式水平集函數進行差集布爾運算,得到空心多孔結構數學表達式ΦH。差集布爾運算數學表達式為

(9)

為獲得混合多孔結構,對兩種極小曲面的混合水平集函數進行并集布爾運算,得到混合多孔結構的數學表達式Φhyb。并集布爾運算數學表達式為

(10)

2 三維數值均勻化法

三維數值均勻化法廣泛應用于計算周期性多孔微結構的宏觀等效屬性?；跀抵稻鶆蚧ǖ淖兠芏韧負鋬灮椒軌蛏纱罅恐虚g密度,滿足多孔結構的材料分布需求。因此,本文采用三維數值均勻化法計算極小曲面多孔結構的宏觀等效彈性張量。

2.1 極小曲面多孔結構等效彈性屬性分析

基于均勻化理論,極小曲面多孔結構的等效彈性張量可表示為

(11)

為得到擾動位移χij,構建均勻化平衡方程的矩陣形式如下

Kχij=fij

(12)

剛度矩陣為

(13)

虛擬載荷為

(14)

式中:Be為單元應變-位移矩陣;εij為6個單位應變,ε11=(1,0,0,0,0,0)T,ε22=(0,1,0,0,0,0)T,ε33=(0,0,1,0,0,0)T,ε12=(0,0,0,1,0,0)T,ε23=(0,0,0,0,1,0)T,ε13=(0,0,0,0,0,1)T。

將擾動位移代入式(11),得到多孔結構的宏觀等效彈性張量。由于極小曲面多孔結構為體心立方結構,具有高度對稱性,其宏觀等效彈性張量為6×6的對稱矩陣,僅有3個獨立變量,簡化形式如下式

(15)

2.2 基于徑向基函數的參數化等效彈性屬性

徑向基函數由于其插值效率高、收斂性好、插值系統解的唯一性等優點,被廣泛應用于離散數據的插值擬合[26-27]。宏觀等效彈性張量關于混合權重因子的函數曲線如圖5所示。采用徑向基函數建立多孔結構宏觀等效彈性張量關于權重因子α的函數,表達式如下

(a)P型多孔結構

(b)I-WP型實心多孔結構

(c)I-WP型空心多孔結構

(d)I-WP&P型混合多孔結構圖5 宏觀等效彈性張量關于混合權重因子的函數曲線Fig.5 Function curves of the macroscopic equivalent elasticity tensor with respect to the hybrid weight factor

(16)

(17)

其中γ是形狀參數,即水平集網格體積的倒數。

同理,基于高斯徑向基函數建立多孔結構體積分數關于權重因子α的函數V(α)表達式如下

(18)

式中:vi為擴展系數,本文為體積分數。體積分數關于權重因子的函數曲線如圖6所示。

圖6 體積分數關于權重因子的函數曲線Fig.6 Curves of the volume fraction with respect to the weight factor

3 梯度多孔結構高階連續建模

梯度多孔結構的設計過程中,多孔單胞間的連續性至關重要,良好的連續性有助于降低梯度多孔結構在低連通區域的應力集中[28-30]。

極小曲面的隱式水平集函數通過歐拉網格定義,結構邊界為局部水平集函數的零水平面。局部水平集函數可由混合水平集函數進行定義,如下式

(19)

式中:M為混合建模的極小曲面數;X=(x,y,z)為高維物理空間坐標;βj為第j個極小曲面的權重函數,如下式

(20)

其中Ωj為極小曲面j的區域。

(a)權重函數曲線

(b)原始梯度多孔結構圖7 局部水平集權重函數Fig.7 The weight function of local level set

局部水平集權重函數如圖7所示。由圖7b可知,由局部水平集函數直接描述的結構存在幾何突變的特征。為保證不同多孔單胞間的光滑過渡,采用Heaviside函數構造局部插值模型,實現了梯度多孔結構的高階幾何連續。局部插值模型數學表達式如下

H(Φj)=

(21)

式中:ζ為一個正極小值,用于避免數值奇異,Δ是Heaviside近似的半帶寬。高階連續局部插值模型如圖8所示。

(a)局部插值模型曲線

(b)基于局部插值模型的梯度多孔結構圖8 高階連續局部插值模型Fig.8 The high-order continuity local interpolation model

4 梯度多孔結構優化設計

4.1 拓撲優化模型

為實現極小曲面梯度多孔結構剛度拓撲優化,建立了體積約束下混合權重因子α為設計變量、結構柔度最小化的優化模型,如下式

(22)

式中:N為設計域內的單元數;目標函數J為梯度多孔結構的柔度;K、U和F分別為結構的全局剛度矩陣、全局位移向量和全局載荷向量;V(α)為設計域的體積約束;V0為設計域的體積;f為體積分數。

4.2 靈敏度分析

通過獲取目標函數的梯度信息來驅動優化算法有效搜索給定設計域內最優材料分布是拓撲優化中的關鍵一步[19,31]。在拓撲優化中,梯度信息通常被稱為敏度信息。因此,采用基于梯度信息的優化準則算法求解優化模型式(22),需要計算目標函數的梯度信息,目標函數和約束條件對設計變量的一階導數推導如下。

對于結構柔度,在目標函數中代入結構場平衡方程KU=F,有J(α)=FTU,在確定的載荷條件下F為常量,結構柔度關于設計變量的導數為

(23)

由上式可得

(24)

對結構場平衡方程兩邊求導可得

(25)

將式(25)代入式(24)中,結構柔度關于設計變量的導數為

(26)

其中

(27)

將式(27)代入式(26)中,結構柔度關于設計變量的導數為

(28)

根據式(16),有

(29)

根據式(17),高斯徑向基函數關于設計變量的導數如下式

(30)

5 數值案例及實驗分析

本文通過以下數值案例驗證多孔結構剛度拓撲優化方法的有效性,每個算例討論實心多孔結構(方案1)、空心多孔結構(方案2)、混合多孔結構(方案3)3種不同代表性體積單胞的結構柔度最小化問題。假定構成實體材料的彈性模量為1、泊松比為0.3,優化迭代中前后目標函數差值小于0.001或迭代次數達到200次時,優化結束。

5.1 三維懸臂梁結構

5.1.1 數值案例 懸臂梁結構邊界條件如圖9所示,固定336×42×168的懸臂梁結構左端面,結構右端面下邊界處施加F=-1的均布載荷。將宏觀結構離散為16×2×8個八節點六面體單元,初始結構體積分數為0.28,開展拓撲優化設計。

圖9 懸臂梁結構邊界條件Fig.9 Boundary condition of a cantilever beam structure

極小曲面梯度多孔結構如圖10所示,優化的結構呈現明顯的梯度分布,實現了功能的梯度變化,且滿足高階連續。對于方案1,通過優化權重因子α,初始結構柔度為4 635,優化后結構柔度為1 220,結構剛度提升約為74%;對于方案2,通過優化權重因子α,改變了空心多孔結構壁厚,初始結構柔度為2 857,優化后結構柔度為1 085,結構剛度提升約為62%;對于方案3,通過優化權重因子α,較大體積模量的P型多孔結構向較大剪切模量的I-WP型多孔結構轉變,初始結構柔度為1 897,優化后結構柔度為1 268,結構剛度約提升33%。數值計算結果表明,在未優化的初始均勻結構中,多孔結構的剛度性能優劣依次為混合多孔結構、空心多孔結構、實心多孔結構;相比于未優化的均布的極小曲面多孔結構,優化的梯度多孔結構的剛度得到了顯著提升。

(a)方案1

(b)方案2

(c)方案3圖10 極小曲面梯度多孔結構Fig.10 TPMS-based graded cellular structures

5.1.2 魯棒性分析 在結構實際工作期間,載荷往往不是恒定的,在非預期邊界條件下的結構性能是衡量結構魯棒性的重要標準。為驗證多孔結構的魯棒性,在非預期載荷的情況下對比了本文的3種方案、均勻結構和傳統拓撲優化的實體結構的結構柔度變化。非預期載荷的大小與預期載荷相等,方向沿y軸偏轉1.37°,受非預期載荷懸臂梁結構的邊界條件如圖11所示。在邊界條件不變的情況下計算預期與非預期載荷的結構柔度,結果如圖12所示。

圖11 受非預期載荷懸臂梁結構的邊界條件Fig.11 Boundary conditions of the cantilever beam structure with unexpected load

由圖12可知,傳統實體結構的結構柔度由預期載荷下的1 664升高至1 923,結構柔度上升了259。實心多孔結構的結構柔度由1 219升高至1 416,結構柔度上升了180?？招亩嗫捉Y構的結構柔度由1 085升高至1 265,結構柔度上升了197?；旌隙嗫捉Y構的結構柔度由1 268升高至1 469,結構柔度上升了201。均勻結構的結構柔度由4 634升高至5 057,結構柔度上升了423。數值計算結果表明,相比傳統實體結構和均勻多孔結構,本文所設計的梯度多孔結構在非預期載荷下的結構柔度變化較小,可以保證較高的結構剛度,具有良好的魯棒性。

圖12 預期與非預期載荷的結構柔度Fig.12 Structural compliance under unexpected and expected loads

5.2 三維Michell梁結構

5.2.1 數值案例 Michell梁結構邊界條件如圖13所示,固定尺寸為294×42×82的Michel梁結構底部兩端,結構的上端面中部處施加F=-1的均布載荷。將宏觀結構離散為14×2×4個八節點六面體單元,初始結構體積分數為0.28,為簡化計算過程,僅對總設計域的右半部分結構開展拓撲優化設計。

圖13 Michell梁結構邊界條件Fig.13 Boundary condition of the Michell beam

圖14給出了Michell梁梯度多孔結構。方案1中初始結構柔度為1 218,優化后結構柔度為310,結構剛度提升約為74%;方案2中初始結構柔度為797,優化后結構柔度為293,結構剛度提升約為63%;方案3中初始結構柔度為625,優化后結構柔度為426,結構剛度提升約為32%。數值計算結果表明,優化的梯度多孔結構的剛度得到了顯著提升,且不同代表性體積單胞的Michell梁結構剛度優化程度與對應的懸臂梁結構剛度優化程度基本一致。

(a)方案1

(b)方案2

(c)方案3圖14 Michell梁梯度多孔結構Fig.14 Michell beam with graded cellular structures

5.2.2 Michell梁結構的三點彎曲實驗分析 為進一步揭示梯度多孔MBB梁的力學性能,對3種優化的MBB梁結構、均勻I-WP實體多孔結構進行三點彎曲實驗。由于數值案例簡化了結構設計域,在進行實驗前,通過對稱平面將優化結構鏡像為完整MBB梁。為保證實驗結果的準確性,實驗的結構尺寸和體積分數與5.2.1節中數值案例一致。選EOS-P760型3D打印機,采用選擇性激光燒結技術,制造了4種MBB梁結構試樣,不同方案的3D打印試樣如圖15所示。結構材料為PA2200,彈性模量為741 MPa,泊松比為0.3,屈服強度為54 MPa。

(a)均勻多孔結構

(b)方案1

(c)方案2

(d)方案3圖15 不同方案的3D打印試樣Fig.15 3D-printed specimens

在溫室條件下,采用萬能試驗機進行三點彎曲實驗,在保證結構邊界條件一致的情況下,以50 mm/min的動態載荷加載,加載時間為8 s,不同方案的實驗平臺如圖16所示。

(a)均布多孔結構

(b)方案1

(c)方案2

(d)方案3圖16 不同方案的實驗平臺Fig.16 The experimental platform in different schemes

彎曲載荷-位移曲線如圖17所示。由圖17可知,4種結構所受到的載荷與位移呈線性變化。相對于均布的I-WP型均布結構,本文所設計的梯度結構的斜率更大,具有更優的承載特性,從而驗證了本文方法的有效性。由于實驗樣件的制造誤差以及I-WP實心梯度結構與空心梯度結構的理論柔度相差較小,從而造成兩種梯度結構的實驗曲線斜率相近,但實驗結果的變化趨勢符合理論分析的預期。

圖17 彎曲載荷-位移曲線Fig.17 The load-displacement curves

6 結 論

本文提出了一種隱式曲面梯度多孔結構拓撲優化設計方法。通過混合水平集函數實現了極小曲面多孔結構混合參數化建模。結合數值均勻化法建立柔度最小化的拓撲優化模型,基于局部插值模型,實現了高階連續的梯度多孔結構優化設計。數值案例和實驗結果表明,所設計的梯度多孔結構較均布極小曲面多孔結構具有更優的魯棒性和承載特性,且單胞構型的不同會造成梯度多孔結構功能特性差異。所提方法能有效實現多孔單胞構型對結構功能特性的調控,豐富了多孔結構的力學內涵。