梯度空心多孔結構優化設計方法

田啟華 舒正濤 付君健 杜義賢 周祥曼 田 磊

1.三峽大學機械與動力學院，宜昌，443002

2.水電機械設備設計與維護湖北省重點實驗室，宜昌，443002

0 引言

空心多孔結構是基于實心多孔結構改進的新形式，以實心多孔結構為基礎，去除內部材料，保留一定外殼厚度，實現結構的空心化[1]?？招亩嗫捉Y構具有承載、吸能、吸聲以及隔振等優良的性能[2-5]，廣泛應用于汽車、船舶和航空航天等領域。然而傳統的建模方法難以實現空心結構復雜的幾何特征描述，且均勻的周期性結構無法充分發揮其力學性能，在實際工程應用中存在一定的局限性，因此，需要進一步探索和發展空心多孔結構參數化建模技術，實現結構功能特性的調控。

先進制造技術的發展為空心多孔結構的應用提供了強有力的支撐。目前，可采用焊接成形[6-7]、熱膨脹模塑成形[8]和增材制造[9-10]等技術，制備出菱形、金字塔型、四面體型及三周期極小曲面等空心多孔結構?？招亩嗫捉Y構壓縮和剪切工況是工程應用中最常見的應用場景，通過研究不同材料的單胞構型發現，空心多孔結構的壓縮和剪切強度優于實心多孔結構，特別是面外壓縮強度，是相同密度下實心多孔結構的2倍以上[6-7]。有限元數值仿真也表明，相對密度為0.9%～5.8%的空心金字塔型多孔結構的壓縮和剪切強度可達相同密度實心多孔結構的3～5倍[11]。盡管在空心多孔結構的制備上具有成熟的制造技術，但目前涉及空心多孔結構的優化設計方法較少[12-13]。大多研究都采用形狀單一的空心圓管或者方管進行設計，通過調整空心桿件的尺寸、偏角或壁厚等幾何參數，實現空心多孔結構力學性能的調控[1，4，14]。然而，空心多孔結構的力學性能與其幾何構型存在復雜的耦合關系，進一步研究參數化建模技術，實現空心多孔結構參數對幾何構型和力學性能的調控是一個亟待解決的關鍵問題。

此外，空心多孔結構的分布形式對力學性能也有重要的影響。相對密度均勻的周期性多孔結構無法充分發揮空心結構單胞的性能優勢，通過空心多孔結構的梯度設計可進一步提高空心結構力學性能。功能梯度多孔結構的性能依賴于周期性細/微觀結構的拓撲構型，無需改變結構的材料組分即可實現結構特定的功能設計[15]。一些學者提出了關于極小曲面的功能梯度多孔結構建模方法，如B樣條插值曲線和多項式插值函數[16-17]等，實現了多孔結構的梯度漸變。梯度多孔結構的設計實際上是考慮設計域內材料屬性的連續變化，本質上與拓撲優化方法類似[18-19]，故可利用拓撲優化方法開展梯度多孔結構設計[20-21]，如密度映射方法[22-23]、材料/結構一體化優化設計方法[24-25]等。上述幾類方法豐富了梯度結構的設計理論和工程應用，但仍存在一些亟待解決的問題，如多孔結構整體體積分數不易準確控制、結構邊界連續性難以保證。需要說明的是，上述研究大多側重于梯度結構吸能特性的提升，結構的承載性能還存在一定的不足[26]。因此，引入性能更優的空心多孔結構單胞，研究穩定、高效的梯度空心多孔結構優化設計方法，將為空心多孔結構的創新設計和承載性能的提升提供一個新的思路。

為實現空心多孔結構的參數化建模、提高多孔結構的承載性能，本文提出了一種梯度空心多孔結構優化設計方法。對I-WP型三周期極小曲面進行布爾運算，實現了空心多孔結構單胞的隱式建模；通過構造混合水平集函數，建立了梯度多孔結構的參數化模型；利用水平集演化出空心多孔結構在空間中最優材料分布形式，得到了具有良好連續性和優良承載性能的梯度空心多孔結構，并通過數值案例和實驗驗證了本文所提方法的可行性與有效性。

1 I-WP型空心多孔結構參數化建模

水平集方法的基本思想是將曲線(或曲面)隱式地表達為一個高維函數Φ(X)的等值面，通過追蹤高維函數Φ(X)的運動來達到描述曲線(或曲面)變化的目的[27]。其數學表達式為[28]

?Ω(X)={X：Φ(X)=0，X∈D}

(1)

式中，X為高維空間的物理坐標；D為包含曲線?Ω(或曲面的)的空間；Ω為幾何結構的實體區域；?Ω為實體區域與空洞區域的邊界。?Ω=ΓD∪ΓN∪Γf，包含Dirichlet邊界ΓD、Neumann邊界ΓN和無牽引力邊Γf。

水平集函數隱式地描述了結構形狀Ω，定義如下：

(2)

圖1所示為水平集函數的隱式描述方式，對于三維的結構邊界，需要四維的水平集函數來描述。

圖1 四維水平集函數及零等值面

三周期極小曲面結構具有低密度、高比強度、設計性強的結構優勢以及良好的制造優勢，如具備自支撐能力、良好的排屑性能[29]，是梯度多孔結構的主要代表性體積單胞之一。本文將采用I-WP型極小曲面作為宏觀結構的代表性體積單胞，實現梯度空心多孔結構的參數化建模。

極小曲面的數學表達式是一種三維隱式水平集函數，在I-WP型極小曲面函數中引入水平參數r可控制極小曲面體積分數的變化[30]，其表達式為

(3)

式中，ФI(x，y，z，r)為極小曲面數學表達式；x，y，z為高維空間的物理坐標；L為極小曲面單胞的邊長。

為了設計I-WP型空心多孔結構，定義兩個I-WP型極小曲面，水平參數分別定義為r1和r2，得到兩個體積分數不同的極小曲面多孔結構。將水平參數為r2的極小曲面多孔結構作為基礎結構，通過差集布爾運算，得到I-WP型空心多孔結構的數學表達式[31]。差集布爾運算的數學表達式為

Φ=min(Φ(x，y，z，r2)，-Φ(x，y，z，r1))

(4)

通過調整水平參數r1可以得到不同壁厚的空心單胞，布爾運算過程如圖2所示。

圖2 I-WP型空心多孔結構布爾運算過程

極小曲面是隱式函數，其隱式特性使其具有布爾運算簡便、結構邊界清晰、拓撲可變能力強等特點，與CAD/CAE等軟件的結合也較為方便。將水平集函數的零等值面轉化為模型三角面片的幾何信息即每個三角面片的法矢量和頂點坐標，由此可得到二進制STL文件，進而對模型進行后處理。模型后處理過程如圖3所示，首先對原始STL模型表面三角面片上存在的孔洞、縫隙、重疊或相交等缺陷進行修復，使其滿足模型實體化要求；再將模型表面三角面片縫合，形成實體。實體化是模型在CAD軟件中進行布爾運算、打孔、修剪等操作的前提，實體化后的模型也便于與CAE軟件結合，進行力學仿真等操作。

圖3 I-WP型空心多孔結構模型處理

2 空心多孔結構彈性性能分析

為說明I-WP型空心多孔結構的性能優勢，采用數值均勻化法[32]計算I-WP型空心多孔結構和實心多孔結構的等效彈性矩陣，揭示空心多孔結構的宏觀等效力學特性。

多孔結構等效彈性張量CH的表達式為

(5)

式中，|V|為單胞的體積；Ve為單元e的體積；I為六階單位矩陣；Be為單元應變-位移矩陣；χe為單元位移向量；De為單元本構矩陣。

基于I-WP型空心多孔結構的高度對稱性，其等效彈性矩陣可以簡化為

(6)

使用等效彈性矩陣中的分量來計算多孔結構的等效彈性模量EH、等效剪切模量GH和等效泊松比υH：

EH=(C112+C11C12-2C122)/(C11+C12)

(7)

GH=C44

(8)

(9)

剛度θ作為多孔結構彈性性能的度量指標[33]：

(10)

(11)

(12)

(13)

(14)

式中，EHSU、GHSU和νHSU分別為引入HS理論極限的體積模量、剪切模量和泊松比；ρ為多孔結構的相對密度；Gb、Kb分別為材料的剪切模量和體積模量。

式(10)～式(14)中，取材料彈性模量(量綱一化)為1，泊松比為0.3，則Gb=1/2.6，Kb=1/1.2。

采用I-WP型空心/實心單胞的梯度多孔結構剛度對比數據如圖4所示，數據結果表明，I-WP型空心結構在任意體積分數下的剛度均高于實心多孔結構，但隨著體積分數的增大，實心單胞剛度逐漸逼近空心單胞。

圖4 不同體積分數的結構剛度對比

為進一步說明I-WP型空心多孔結構單胞的力學性能優勢，分別對體積分數為0.2～0.7的I-WP型空心/實心單胞進行有限元分析。材料設置為尼龍PA2200，與后續的仿真和實驗材料保持一致。具體材料屬性如下：彈性模量741 MPa，屈服強度54 MPa，泊松比0.3，密度1020 kg/m3。多孔結構單胞的邊界條件如圖5a所示，單胞下平面為固定約束，上平面施加豎直向下的均布載荷，F=500 N。

在相同載荷下，兩種結構產生的最大位移與最大應力對比如圖5b、圖5c所示。單胞位移與應力變化趨勢較為相似，隨著體積分數的增大，空心、實心單胞的性能差距逐漸縮小，但空心單胞產生的位移與應力始終比實心單胞小，說明了空心多孔結構單胞具有更優良的力學性能。

(a)邊界條件

3 空心多孔結構混合水平集描述

3.1 混合水平集函數構造

三周期極小曲面多孔結構采用參數化數學函數進行幾何描述，具有較強的可設計性[16]：①通過調整極小曲面單胞的邊長L或水平控制參數r，可控制極小曲面構型和體積分數的變化；②對不同類型的三周期極小曲面進行混合設計，可得到具有新型結構特征的多孔結構?；谝陨咸攸c，可以引入凸優化理論中的仿射概念[34]：對于集合Q?R，如果通過集合Q中任意兩個不同點之間的直線仍然在集合Q中，則稱集合Q為仿射集。

利用相似的概念，對兩種不同的空心極小曲面結構的混合設計過程進行定義，其數學表達式為

ΦH=sΦ1+(1-s)Φ2

(15)

式中，ΦH為具有新型結構特征的空心極小曲面函數；Φ1和Φ2為兩種預定義的空心極小曲面水平集函數；s為設計參數，即兩種預定義空心極小曲面結構的相對權重。

基于式(15)，預定義兩種空心多孔結構單胞的水平集函數，并引入權重系數，構造梯度空心多孔結構的混合水平集函數，實現多孔結構的參數化建模。采用隱式水平集函數進行結構幾何描述：

(16)

式中，w(t)為權重系數矢量，是關于宏觀時間t的變量，每一元素值取值范圍為[0,1]；x為設計域D的內空間物理坐標；φ1(x)和φ2(x)分別為預定義的兩種空心多孔結構單胞的水平集函數。

通過調整權重系數w(t)，可以靈活地控制空心多孔結構單胞的壁厚，同時改變了單胞的體積分數。為進一步說明權重系數對空心多孔結構單胞壁厚和體積分數的影響，預定義兩個量綱一尺寸為50×50×50的空心多孔結構單胞，根據式(3)、式(4)，定義可變水平參數r1為13.3、-4.69，固定參數r2=-11，得到兩個空心多孔結構單胞φ1(x)和φ2(x)，其體積分數為0.7、0.15，單胞壁厚分別為10.44、1.94(量綱一單位)。通過調整權重系數發現，權重系數對單胞壁厚的影響程度主要分為三個階段，即權重系數為0～0.2、0.3～0.8、0.9～1。權重系數與單胞壁厚、體積分數的詳細關系如表1所示。I-WP型極小曲面空心單胞幾何特征較為復雜，因此，權重系數和單胞壁厚并不是嚴格的線性關系。

表1 權重系數與單胞參數關系

3.2 梯度空心多孔結構C0連續性建模

在梯度空心多孔結構設計過程中，良好的連續性是多孔結構的力學性能與后期制造性能的基本保障[35-36]。為了保證多孔結構邊界的連續性過渡，本文基于混合水平集方法構造全局更新-局部插值的設計變量更新策略，實現梯度空心多孔結構C0連續性建模。

將設計域離散為M個單胞，具有m個節點。任意一個單胞C可以進一步離散為若干八節點六面體單元，如圖6所示。

圖6 設計域離散與設計變量定義

全局設計變量可表示為[35]

w(t)=[w1(t)w2(t) …wm(t)]Τ

(17)

式中，w(t)為定義在單胞節點上的宏觀設計變量，即全局權重系數；m為設計域單胞節點數，即全局設計變量的數目。

對于編號為C的三維多孔結構單胞，8個節點上的權重系數是全局權重系數向量的子向量：

(18)

(19)

式中，省略位置的數值均為0。

(20)

對于八節點六面體線性單元，其形函數Ni定義為

(21)

其中，i(i=1，2，…，8)為單元節點序號，ξ0=ξiξ，η0=ηiη，ζ0=ζiζ，ξi=±1，ηi=±1，ζi=±1，ξ、η、ζ為局部坐標系下單元內部點坐標。

單胞C上所有節點的權重系數wC(t)可表示為

(22)

在不添加額外的幾何約束或對優化結構進行后處理的情況下，結構邊界會出現幾何突變，相鄰單胞邊界的連續性無法得到保證，如圖7a所示。由圖6所示的設計變量定義可知，相鄰單胞在公共節點上定義了相同的權重系數，對全局權重系數進行插值時，能夠自然地保證線性插值函數的連續性，進而實現了梯度多孔結構C0連續，而不需要添加額外的幾何約束。圖7b所示為引入插值策略后的梯度多孔結構，可以發現，相鄰空心結構單胞邊界具有良好幾何連續性，實現了單胞壁厚的連續漸變。

(a)原始梯度空心多孔結構

4 拓撲優化數學模型與靈敏度分析

4.1 拓撲優化數學模型

以結構剛度性能為目標函數，多孔結構體積分數為約束條件，建立基于混合水平集方法的梯度空心多孔結構拓撲優化數學模型：

(23)

彈性平衡方程的弱形式中，能量雙線性形式aФ(u，v)和載荷線性形式lФ(v)分別為

(24)

(25)

其中，f為應用在邊界?Ω上的牽引力，p為體積力，|Ф|為水平集函數梯度的模。δ(Ф)為Dirac函數，是Heaviside函數的導數形式：

(26)

這里，γ=0.001，Δ為Heaviside函數近似的半帶寬。

4.2 靈敏度分析

本文采用形狀導數敏度[15，37]分析方法推導目標函數、能量雙線性形式、載荷線性形式關于宏觀時間變量t的微分形式，進一步化簡為

(27)

(28)

式中，vn為法向速度；n為法向單位向量。

以全局權重系數w(t)作為設計變量，對Hamilton-Jacobi偏微分方程進行時空解耦，從而將水平集函數參數化。將混合水平集函數代入Hamilton-Jacobi偏微分方程，得到僅關于時間t的常微分方程，其數學表達式為

(29)

可得水平集函數邊界演化的法向速度：

(30)

將式(30)代入式(27)得

(31)

采用鏈式法則對目標函數J(Φ)直接求其關于時間變量t的導數：

(32)

比較式(31)和式(32)可以得到目標函數J(Φ)關于設計變量w(t)的靈敏度：

(33)

同理得到體積約束g(Φ)關于設計變量w(t)的靈敏度：

(34)

5 數值案例

梯度空心多孔結構設計流程如圖 8所示。首先對I-WP型實心單胞進行布爾運算，獲得空心單胞，再利用混合水平集函數實現多孔結構參數化建模，最后求解拓撲優化模型，得到幾何和功能在空間呈梯度分布的空心多孔結構。

圖8 梯度空心多孔結構設計流程

本文將采用OC(optimality criteria)作為優化求解算法，對數學模型進行求解。OC算法具有簡單直接、迭代收斂快的特點，對單約束條件下的大規模結構優化設計問題求解尤為高效[38]。在優化算法中，給定收斂準則如下：

(35)

式中，J為目標函數值；n為當前迭代步數；nmax為最大迭代步數。

為了驗證本文提出方法的可行性與有效性，分別基于I-WP型空心和實心多孔結構單胞，給出三維懸臂梁梯度多孔結構拓撲優化算例。在本節的優化過程中，本文給出的量均為量綱一單位，在有限元分析中，按照實際模型尺寸以及負載為參數給出了相應的剛度。

5.1 采用空心多孔結構單胞的懸臂梁結構優化

圖9所示為三維懸臂梁結構設計域，設計域長寬高分別為L=120，W=20，H=50。設計域的右端面為固定約束，上平面距離左邊界L/12處施加豎直向下的均布載荷F=-10。設計域離散為12×2×5個單胞，每個單胞進一步離散為50×50×50個八節點六面體單元。設計域體積分數約束為0.45。選取I-WP型空心多孔結構作為預定義多孔結構單胞，如圖10所示，其體積分數分別為0.2、0.6，單胞外形尺寸一致，壁厚尺寸不同。

圖9 三維懸臂梁設計域

圖10 預定義空心多孔結構單胞

圖11所示為采用I-WP型空心單胞的懸臂梁拓撲優化迭代曲線，隨著優化的進行，結構柔度逐漸減小，經過61次迭代后，結構柔度從15 064.27收斂至8551.08，體積分數收斂至0.45。根據優化迭代曲線可知，提出的梯度空心多孔結構混合水平集優化算法運行穩定，收斂性較好。

圖11 采用I-WP型空心單胞的懸臂梁結構優化迭代曲線

懸臂梁結構最優構型如圖12所示，結構輪廓清晰，在主要傳遞力路徑上，空心單胞壁厚較厚，其他區域空心單胞體積分數較小，保留最小壁厚尺寸。優化后的結構在幾何上和密度上呈現出明顯的梯度分布，單胞自身根據外部載荷和邊界條件在空間中表現出不同的方向梯度。從這個意義上說，該方法不僅獲得了結構最佳材料空間分布，還實現了最優梯度單胞構型設計。此外，在優化過程中，由于相鄰單胞公共節點上的權重系數是相同的，故保證了線性插值函數的連續性，從而實現了單胞幾何邊界的連續過渡，也使得優化后的梯度空心多孔結構具有較好的力學性能。

圖12 采用空心單胞的懸臂梁最優拓撲構型

設計域中單胞的梯度變化是通過單胞節點上的權重系數控制的，在結構主要的傳遞力路徑上，權重系數分布的數值較大，生成壁厚較厚的單胞。對于設計域中承載作用貢獻小的區域，權重系數分布的數值較小，從而填充薄壁單胞。在傳統的拓撲優化中，優化后的結構只保留了受力材料，當傳遞力路徑位置上的材料被破壞時，會嚴重影響結構整體的力學性能[16]。而在本文所提出的方法中，通過控制預設單胞的體積分數，使得設計域內分布一系列壁厚不同的空心梯度單胞，極大地提高了結構穩定性和整體剛度性能。

5.2 采用實心多孔結構單胞的懸臂梁結構優化

為了說明優化方法的有效性，采用實心多孔結構單胞對懸臂梁結構進行優化設計。選取圖13所示的I-WP型實心結構作為預定義多孔結構單胞，其體積分數分別為0.2、0.6。圖14所示為采用實心單胞的懸臂梁結構優化迭代曲線，經過61次迭代后，目標函數和體積約束曲線收斂。初始均勻多孔結構的結構柔度為53 731.06，優化后的結構柔度為10 954.80。

圖13 預定義實心多孔結構單胞

圖14 采用I-WP型實心單胞的懸臂梁結構優化迭代曲線

由于預定義的I-WP型實心多孔結構單胞最小體積分數為0.2，故多孔結構被填充到整個設計域內。在傳遞力路徑上的單胞多孔結構單胞的體積分數較大，其他區域則分布較少的材料，結構在幾何和密度上也呈現出明顯的梯度分布。由圖15發現，體積分數不同的單胞邊界之間過渡連續，可有效避免多孔結構承載時出現的應力集中現象，也保證了后期的可制造性。

圖15 采用實心單胞的懸臂梁最優拓撲構型

在上述兩個算例中，分別采用了I-WP型空心和實心單胞對懸臂梁結構進行優化。根據優化結果可知，采用本文提出的混合水平集拓撲優化算法，均可不同程度地提高結構的承載性能，得到幾何上和密度上呈梯度分布的多孔結構。根據圖11和圖14繪制了兩種優化后的結構柔度對比圖，如圖16所示。在結構體積分數設置為0.45條件下，優化后的梯度結構柔度分別為8551.08、10 954.80。優化結果表明，在相同體積分數下，梯度空心多孔結構的承載能力可提高21.94%，說明本文方法具有顯著的優化效果。

圖16 結構柔度對比

5.3 優化結構的有限元分析

為了對比分析優化后的結構在線彈性范圍內的靜承載性能，分別對兩種結構進行了有限元分析。定義結構的材料為尼龍PA2200，與后續實驗樣件材料保持一致，在本文第2節中給出了材料屬性的具體數值。邊界條件按照圖9設置，繪制了結構載荷與最大位移曲線，如圖17所示。由有限元仿真結果可知，在相同的靜載荷下，梯度空心多孔結構的變形位移明顯小于相同體積分數的實心多孔結構，表明空心多孔結構的承載性能更優。在線彈性范圍內，結構變形量與負載成線性關系，根據載荷-位移曲線可知，采用I-WP型空心和實心單胞的梯度多孔結構整體剛度分別為63.86 N/mm、51.95 N/mm，結構剛度提高了約22.92%。

圖17 載荷-最大位移仿真曲線

6 懸臂梁結構實驗分析

為進一步驗證采用空心或實心單胞的梯度多孔結構承載性能，對優化結構進行實驗分析，用于實驗的模型結構尺寸和體積分數與5.1節中給定的參數保持一致，L=120 mm，W=20 mm，H=50 mm。同時，為了便于夾持和固定樣件，在模型右側設置尺寸為15 mm×20 mm×50 mm的實體區域。選用EOS-P760型3D打印機，采用選擇性激光燒結(selective laser sintering，SLS)技術制造了懸臂梁模型樣件，結構材料為尼龍PA2200，如圖18所示。

(a)采用空心單胞的懸臂梁樣件

根據相關實驗標準，在室溫條件下，采用WDW-100E型電子式萬能試驗機進行懸臂梁實驗。圖19所示為樣件實驗裝置，包括夾持樣件的臺鉗、固定臺、萬能試驗機等。在盡可能保證結構邊界條件一致的情況下，萬能試驗機壓頭以20 mm/min的動態載荷垂直加載，模擬準靜態壓縮條件，加載時間為30 s。利用萬能試驗機數據采集系統記錄加載過程中的位移和相應的載荷。

圖19 實驗裝置

萬能試驗機壓頭向下移動行程約為10 mm，取位移區間0～8 mm為有效行程，樣件的實驗荷載-位移曲線如圖20所示。通過線性回歸計算了樣件的整體剛度，采用I-WP型空心和實心單胞的梯度結構樣件的整體剛度分別為65.89 N/mm、47.92 N/mm。得到的剛度值表明，采用空心單胞的梯度結構剛度性能明顯優于采用實心單胞的梯度結構。此外，結構的強度還可以通過在相同變形下承受的載荷大小來評估。如圖20所示，當萬能試驗機壓頭向下位移達到8 mm時，空心單胞梯度結構承受的載荷為487.4 N，實心單胞梯度結構承受的載荷為383.0 N，說明采用空心單胞的梯度結構具有較好的承載性能。

圖20 兩種樣件的負載-位移曲線

本文分別對兩種結構進行了有限元分析與實驗研究，仿真與實驗的結果具有較好的一致性，其誤差均在5%左右。有限元分析結果與實驗結果均表明，相較于采用實心單胞的梯度多孔結構，本文所設計的梯度空心多孔結構的剛度和強度有顯著提高。

7 結論

本文提出了一種梯度空心多孔結構優化設計方法，對I-WP型三周期極小曲面進行布爾運算實現了空心多孔結構單胞的參數化設計；引入凸優化理論的仿射概念，構造混合水平集函數，實現了梯度空心多孔結構的參數化建模。對優化結構進行有限元仿真分析，并采用選擇性激光燒結技術制造了模型樣件，開展了力學實驗。數值案例和實驗結果表明，所提方法能夠有效地解決空心多孔結構復雜的幾何特征建模問題，顯著提高多孔結構的承載性能。在實際工程問題中，結構需要面對復雜的物理場問題，可在本文方法中考慮結構的吸能特性、屈曲穩定性或散熱特性等，實現面向多功能、多物理場的梯度空心多孔結構設計。