基于梯度有限元的異質材料實體優化設計

吳振華， 趙 韓， 董玉德， 柳吉慶

（合肥工業大學機械與汽車工程學院，安徽 合肥 230009）

異質材料實體（HO）是指由兩種或兩種以上材料組成的實體，若材料組分和屬性在空間光滑變化則也稱為梯度功能材料（FGM），其作為一種新型復合材料在化工、電子、機械、生物醫學工程等領域得到日益廣泛的應用[1-2]。近年來，隨著超聲波固結（UC）、激光近凈成形（LENS）、直接金屬沉積（DMD）等快速成型工藝技術的發展，允許材料組分以區、層、點為單位變化，異質材料實體的制備更具有靈活性和經濟性[3-4]。相應的異質材料實體的計算機圖形學表示、設計和分析技術也日趨成熟[5]，吳曉軍等[6]提出了一種基于CAD零件離散體素模型的異質材料建模方法，寇欣宇等[7]提出了一種基于非流形幾何與特征樹的異質材料實體可視化方法，Yang等[8]提出了基于 B樣條曲線和圖形放樣的異質材料實體設計和有限元分析流程，Kim等[9]提出了梯度有限元的概念和理論模型，用以表示各向同性或正交異性的非均質梯度功能材料。

對材料的空間分布進行優化則以拓撲優化為代表，主要采用均勻化方法、水平集方法和漸進結構優化方法等[10-11]。Guest等[12]將節點變量引入拓撲優化，但采用鄰域映射法將節點解轉換為單元解進行后續計算。優化異質材料實體本質上就是確定組成材料在空間中的體積分數的分布。Huang等[13]基于復合材料微觀結構分析空間點的有效屬性并通過優化微觀結構單元的幾何尺寸和偏置角度來優化宏觀材料的分布。Carbonari等[14]設計了壓電梯度材料致動器、分類進行靈敏度分析并運用均勻化方法對材料分布進行優化。本文提出了基于節點變量和梯度有限元的材料優化模型、運用結構優化移動漸近線算法實現異質材料實體功能最優的完整設計流程。

1 梯度有限元方法

建立計算模型首先需要建立復合材料的有效屬性與組分材料體積分數的映射關系。異質材料實體的每一個點都是各向同性的復合材料，本文采用Hashin-Shtrikman下邊界作為插值模型來擬合由快速成型（RP）技術制備出來的梯度材料的彈性模量[15-16]。擬合的復合材料的有效楊氏模量算式為

式中， x ∈ [0 ，1]為組分材料的體積分數，E1≥E2分別為兩相組分材料的楊氏模量，該公式適用于二維彈性問題。鋼、鋁復合材料有效楊氏模量的幾種插值函數與體積分數變量關系，如圖1所示。Voigt邊界和Reuss邊界適用于每一個點都是各向異性材料的結構體。

以一階拉格朗日四邊形單元（Q4）對二維實體空間劃分網格，每個節點處的體積分數向量作為優化變量，即 x =[x1， …，xnn]，式中nn為網格總節點數。其函數在空間 C0連續，單元內的體積分數分布由雙線性插值決定，保證了函數的協調性和連續性，這種有限元即稱為梯度有限元。節點處的楊氏模量ki由節點處的體積分數xi代入式（1）進行計算，則網格第e個單元內的楊氏模量分布函數為

圖1 鋼、鋁復合材料有效楊氏模量的擬合曲線

式中，ν即泊松比，將式（2）代入即可求解。平面應力問題的梯度有限元的總體剛度矩陣為

式中，ne為網格總單元數， Le為第e個單元的離散關系矩陣， Ke為第e個單元的剛度矩陣，Be為形函數第e個單元的空間偏導數矩陣，將式（3）代入即可求解。此時，根據Ku=f可求解有限元模型得出位移向量u并進行后處理。

拓撲優化通常采用單元變量進行分析和優化，即每個單元內材料是均質的。圖2比較了基于單元變量的均質四邊形單元(a)和基于節點體積分數變量的梯度四邊形單元(b)的可視化效果。

圖2 均質單元和梯度單元的可視化

2 材料分布優化模型

材料分布的優化問題可抽象為以下數學模型

本文求解計算基于節點變量、梯度單元的有限元模型，上述優化問題以及設計異質材料實體的過程如圖3所示。

2.1 靈敏度分析

移動漸近線算法是一階方法，需要目標函數和約束函數的一階偏導數即靈敏度。本文用解析方法求解靈敏度。由式（1）計算得復合材料有效楊氏模量函數對節點變量體積分數xj的靈敏度為

①北京市委組織部把生態清潔小流域建設、農村治污達標率作為各區縣委、政府領導班子屆中考核的指標，明確了地方政府責任。政府出臺《關于推進山區小流域綜合治理和關停廢棄礦山生態修復的意見》（京政辦發〔2006〕66 號），市發改委、市財政局每年從基本建設資金、水資源費和土地出讓金中安排一定比例資金用于生態清潔小流域建設；投資標準從25萬元/km2提高到 50萬元/km2；制定廢棄礦山治理規劃，實現全市現有3 667 hm2廢棄礦山全部生態修復，裸露礦山重新披綠；2014年標準又提高到 65萬元/km2。

圖3 異質材料實體設計優化流程圖

由式（4）計算得有限元模型總體剛度矩陣對節點變量體積分數xj的靈敏度為

式中，nae為包含全局第j個節點的單元數，將式（6）代入即可求解。對平衡方程K(x)u(x ) =F求偏導數得

式中，k為優化迭代次序號，此式具有平衡方程的形式，右邊項視為假載[17]（pseudo-load），將式（7）代入并使用求解有限元模型的同一求解器即可解得

由二維Q4單元劃分網格的異質材料實體的重量寫作

式中， Ae為第e個單元的面積，t為厚度。由此可計算重量對節點變量體積分數xj的靈敏度為

2.2 移動漸近線算法

Svanberg[18]提出移動漸近線（MMA）算法后該算法因能適應結構優化問題的特性而得到廣泛應用，其在廣義結構優化問題上近似精確度高于序列線性規劃（SLP），運算量低于序列二次規劃（SQP），收斂速度快于凸線性化方法（CONLIN）[10]。本文基于梯度單元節點變量的異質材料分布優化問題設計變量數目龐大，適于以MMA算法進行優化。

移動漸近線算法引入如下中間變量

式中Lj及Uj即所謂的移動漸近線，隨著迭代而變化，且對第k次迭代來說始終滿足算法對目標函數和約束函數gi,i = 0,… ,l 在第k次迭代的迭代解 xk處的近似為

式中， r ik為殘值，且有

則優化問題（5）在第k次迭代處的移動漸近線近似可寫作

式中， μ ∈ (0,1)，則原優化問題（5）在第k次迭代處被近似為一個具有顯式、可分特性的簡單的凸優化問題（14），利用拉格朗日對偶可得出在第k次迭代處的最優解并代入下一次迭代，直至滿足優化終止條件。

3 梯度控制算法

在異質材料實體設計方法中，材料的體積分數和物理屬性通過Q4單元的雙線性形函數實現了在二維設計空間中的0C 連續，避免了材料空間分布突變造成應力集中[1]，為了提高該設計優化方法的魯棒性和自主性，有必要使設計者可以直接控制異質材料在空間中的梯度大小。Petersson等[19]將基于均質單元和單元變量的鄰域單元梯度約束引入拓撲優化，提供了一種直接控制局部材料梯度的方法，并且因此減少了數值計算中的棋盤模式和網格依賴現象，但這種算法大大增加了運算量。本文將該局部材料梯度控制方法應用于梯度單元和節點變量，并演化成一種新算法提高計算性能。

圖4為梯度控制算法的網格示意圖，Ωi為節點i的鄰域，鄰域半徑為r，d(i,j)為節點i與節點 j的距離， Ωij∈ ，梯度約束可表示為

圖4 節點鄰域梯度控制示意圖

式中，xi和xj分別為節點i和節點 j處的設計變量， 此式引入了個線性約束函數，直接計算不具備經濟性，因此提出實施如下設計變量自適應下界來實現式（16）的約束

4 優化設計實例

本文以金屬夾鉗為例來驗證說明前述的異質材料實體設計優化方法。以鋼、鋁為兩相組成材料設計如圖5所示的夾鉗，設厚度t=1mm。設定鋼和鋁的楊氏模量分別為200GPa和68GPa，密度分別為7.85g/cm3和2.64g/cm3，兩者泊松比差別細微，簡便起見均設為 0.3。將該夾鉗劃分為Q4單元的網格，以重量為目標函數，A節點和B節點的Y方向位移為約束函數建立優化數學模型如下

式中，W為夾鉗重量，uyA和uyB分別為A點和B點的Y方向位移，x為設計變量，設為各節點處鋼的體積分數向量，初始化為，按前文所述方法進行靈敏度分析、移動漸近線優化，得到的最優材料分布如圖6(a)所示，與圖6(b)所示的由商業軟件 Optistruct計算得到的相同結構的載荷路徑拓撲圖相比，可看到強度較高的鋼的分布與載荷路徑拓撲有異曲同工之處。

圖5 夾鉗的尺寸及載荷工況

圖6 夾鉗的梯度材料最優分布與載荷路徑對比

表1對比了不同材料組成的夾鉗在載荷工況下的重量和位移，鋼、鋁最優分布的異質材料夾鉗在滿足位移約束的前提下重量達到了極值，而位移值又比同樣材料比例組成的均質合金的位移值要小。

表1 不同材料金屬夾鉗重量和位移對比

5 結 論

本文在簡明的梯度有限元概念基礎上提出的啟發式異質材料實體設計優化方法，流程清晰，運算量合理，可以在二維空間得出實體材料的最優分布，從而使設計目標和設計約束達到最優的均衡。這種設計方法避免了基于單元變量的材料分布優化所導致的材料突變，具有材料梯度可控且控制方法運算量小的特點。因為利用了有限元自身的形函數對材料分布進行空間插值從而使得計算過程大大簡化，并且同時達到了材料分布的C0連續。這種設計方法可以向三維空間擴展并廣泛應用到梯度功能材料的分布設計中去。該方法的缺點是只能設計兩相材料實體，如何處理兩相以上的材料的分布優化成為了下一步工作的重點。

[1] Miyamoto Y. Functionally graded materials: design,processing,and applications [M]. Chapman & Hall,1999: 1-27.

[2] Pompe W,Worch H,Epple M,et al. Functionally graded materials for biomedical applications [J].Materials Science and Engineering A,2003,362(1-2):40-60.

[3] Hu Yuna,Blouin V,Fadel G M. Design for manufacturing of 3D heterogeneous objects with processing time consideration [J]. Journal of Mechanical Design,2008,130(3): 371-379.

[4] Zhan Y B,Lin D J,An Q. Slicing method for reverse engineering based on image mosaic [J]. Computer Aided Drafting,Design and Manufacturing,2008,(02):33-38.

[5] Hung SKFC,Deng Z M. Thermal design and optimization of heat pipe radiator for satellites [J].Computer Aided Drafting,Design and Manufacturing,2010,(01): 56-64.

[6] 吳曉軍,劉偉軍,王天然. 三維 CAD零件異質材料建模方法[J]. 機械工程學報,2004,40(5): 111-117.

[7] 寇欣宇,王以忠,彭一準. 基于非流形幾何與特征樹的異質材料實體可視化方法[J]. 計算機輔助設計與圖形學學報,2008,20(4): 532-539.

[8] Yang Pinghai,Qian Xiaoping. A B-spline-based approach to heterogeneous objects design and analysis [J]. Computer-Aided Design,2007,39(2):95-111.

[9] Kim J H,Paulino G H. Isoparametric graded finite elements for nonhomogeneous isotropic and orthotropic materials [J]. Journal of Applied Mechanics,2002,69(4): 502-514.

[10] Bends?e M P,Sigmund O. Topology optimization:theory,methods,and applications [M]. Berlin:Springer Verlag,2003: 10-30.

[11] Rozvany G. A critical review of established methods of structural topology optimization [J]. Structural and Multidisciplinary Optimization, 2009, 37(3):217-237.

[12] Guest J K,Prévost J H,Belytschko T. Achieving minimum length scale in topology optimization using nodal design variables and projection functions [J].International Journal for Numerical Methods in Engineering,2004,61(2): 238-254.

[13] Huang Jinhua,Fadel G M,Blouin V Y,et al.Bi-objective optimization design of functionally gradient materials [J]. Materials & Design,2002,23(7): 657-666.

[14] Carbonari R C,Silva E,Paulino G H,et al. Topology optimization design of functionally graded bimorph-type piezoelectric actuators [J]. Smart Materials and Structures,2007,16: 2605-2620.

[15] Walpole L J. On bounds for the overall elastic moduli of inhomogeneous systems—I [J]. Journal of the Mechanics and Physics of Solids,1966,14(3):151-162.

[16] Bends?e M P,Sigmund O. Material interpolation schemes in topology optimization [J]. Archive of Applied Mechanics (Ingenieur Archiv),1999,69(9):635-654.

[17] Christensen P W,Klarbring A. An introduction to structural optimization[M]. Springer Verlag,2008:93-99.

[18] Svanberg K. The method of moving asymptotes: a new method for structural optimization [J].International Journal for Numerical Methods in Engineering,1987,24(2): 359-373.

[19] Petersson J,Sigmund O. Slope constrained topology optimization [J]. International Journal for Numerical Methods in Engineering,1998,41(8): 1417-1434.