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

（合肥工業(yè)大學(xué)機(jī)械與汽車(chē)工程學(xué)院，安徽 合肥 230009）

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

對(duì)材料的空間分布進(jìn)行優(yōu)化則以拓?fù)鋬?yōu)化為代表，主要采用均勻化方法、水平集方法和漸進(jìn)結(jié)構(gòu)優(yōu)化方法等[10-11]。Guest等[12]將節(jié)點(diǎn)變量引入拓?fù)鋬?yōu)化，但采用鄰域映射法將節(jié)點(diǎn)解轉(zhuǎn)換為單元解進(jìn)行后續(xù)計(jì)算。優(yōu)化異質(zhì)材料實(shí)體本質(zhì)上就是確定組成材料在空間中的體積分?jǐn)?shù)的分布。Huang等[13]基于復(fù)合材料微觀結(jié)構(gòu)分析空間點(diǎn)的有效屬性并通過(guò)優(yōu)化微觀結(jié)構(gòu)單元的幾何尺寸和偏置角度來(lái)優(yōu)化宏觀材料的分布。Carbonari等[14]設(shè)計(jì)了壓電梯度材料致動(dòng)器、分類(lèi)進(jìn)行靈敏度分析并運(yùn)用均勻化方法對(duì)材料分布進(jìn)行優(yōu)化。本文提出了基于節(jié)點(diǎn)變量和梯度有限元的材料優(yōu)化模型、運(yùn)用結(jié)構(gòu)優(yōu)化移動(dòng)漸近線算法實(shí)現(xiàn)異質(zhì)材料實(shí)體功能最優(yōu)的完整設(shè)計(jì)流程。

1 梯度有限元方法

建立計(jì)算模型首先需要建立復(fù)合材料的有效屬性與組分材料體積分?jǐn)?shù)的映射關(guān)系。異質(zhì)材料實(shí)體的每一個(gè)點(diǎn)都是各向同性的復(fù)合材料，本文采用Hashin-Shtrikman下邊界作為插值模型來(lái)擬合由快速成型（RP）技術(shù)制備出來(lái)的梯度材料的彈性模量[15-16]。擬合的復(fù)合材料的有效楊氏模量算式為

式中， x ∈ [0 ，1]為組分材料的體積分?jǐn)?shù)，E1≥E2分別為兩相組分材料的楊氏模量，該公式適用于二維彈性問(wèn)題。鋼、鋁復(fù)合材料有效楊氏模量的幾種插值函數(shù)與體積分?jǐn)?shù)變量關(guān)系，如圖1所示。Voigt邊界和Reuss邊界適用于每一個(gè)點(diǎn)都是各向異性材料的結(jié)構(gòu)體。

以一階拉格朗日四邊形單元（Q4）對(duì)二維實(shí)體空間劃分網(wǎng)格，每個(gè)節(jié)點(diǎn)處的體積分?jǐn)?shù)向量作為優(yōu)化變量，即 x =[x1， …，xnn]，式中nn為網(wǎng)格總節(jié)點(diǎn)數(shù)。其函數(shù)在空間 C0連續(xù)，單元內(nèi)的體積分?jǐn)?shù)分布由雙線性插值決定，保證了函數(shù)的協(xié)調(diào)性和連續(xù)性，這種有限元即稱為梯度有限元。節(jié)點(diǎn)處的楊氏模量ki由節(jié)點(diǎn)處的體積分?jǐn)?shù)xi代入式（1）進(jìn)行計(jì)算，則網(wǎng)格第e個(gè)單元內(nèi)的楊氏模量分布函數(shù)為

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

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

式中，ne為網(wǎng)格總單元數(shù)， Le為第e個(gè)單元的離散關(guān)系矩陣， Ke為第e個(gè)單元的剛度矩陣，Be為形函數(shù)第e個(gè)單元的空間偏導(dǎo)數(shù)矩陣，將式（3）代入即可求解。此時(shí)，根據(jù)Ku=f可求解有限元模型得出位移向量u并進(jìn)行后處理。

拓?fù)鋬?yōu)化通常采用單元變量進(jìn)行分析和優(yōu)化，即每個(gè)單元內(nèi)材料是均質(zhì)的。圖2比較了基于單元變量的均質(zhì)四邊形單元(a)和基于節(jié)點(diǎn)體積分?jǐn)?shù)變量的梯度四邊形單元(b)的可視化效果。

圖2 均質(zhì)單元和梯度單元的可視化

2 材料分布優(yōu)化模型

材料分布的優(yōu)化問(wèn)題可抽象為以下數(shù)學(xué)模型

本文求解計(jì)算基于節(jié)點(diǎn)變量、梯度單元的有限元模型，上述優(yōu)化問(wèn)題以及設(shè)計(jì)異質(zhì)材料實(shí)體的過(guò)程如圖3所示。

2.1 靈敏度分析

移動(dòng)漸近線算法是一階方法，需要目標(biāo)函數(shù)和約束函數(shù)的一階偏導(dǎo)數(shù)即靈敏度。本文用解析方法求解靈敏度。由式（1）計(jì)算得復(fù)合材料有效楊氏模量函數(shù)對(duì)節(jié)點(diǎn)變量體積分?jǐn)?shù)xj的靈敏度為

①北京市委組織部把生態(tài)清潔小流域建設(shè)、農(nóng)村治污達(dá)標(biāo)率作為各區(qū)縣委、政府領(lǐng)導(dǎo)班子屆中考核的指標(biāo)，明確了地方政府責(zé)任。政府出臺(tái)《關(guān)于推進(jìn)山區(qū)小流域綜合治理和關(guān)停廢棄礦山生態(tài)修復(fù)的意見(jiàn)》（京政辦發(fā)〔2006〕66 號(hào)），市發(fā)改委、市財(cái)政局每年從基本建設(shè)資金、水資源費(fèi)和土地出讓金中安排一定比例資金用于生態(tài)清潔小流域建設(shè)；投資標(biāo)準(zhǔn)從25萬(wàn)元/km2提高到 50萬(wàn)元/km2；制定廢棄礦山治理規(guī)劃，實(shí)現(xiàn)全市現(xiàn)有3 667 hm2廢棄礦山全部生態(tài)修復(fù)，裸露礦山重新披綠；2014年標(biāo)準(zhǔn)又提高到 65萬(wàn)元/km2。

圖3 異質(zhì)材料實(shí)體設(shè)計(jì)優(yōu)化流程圖

由式（4）計(jì)算得有限元模型總體剛度矩陣對(duì)節(jié)點(diǎn)變量體積分?jǐn)?shù)xj的靈敏度為

式中，nae為包含全局第j個(gè)節(jié)點(diǎn)的單元數(shù)，將式（6）代入即可求解。對(duì)平衡方程K(x)u(x ) =F求偏導(dǎo)數(shù)得

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

由二維Q4單元?jiǎng)澐志W(wǎng)格的異質(zhì)材料實(shí)體的重量寫(xiě)作

式中， Ae為第e個(gè)單元的面積，t為厚度。由此可計(jì)算重量對(duì)節(jié)點(diǎn)變量體積分?jǐn)?shù)xj的靈敏度為

2.2 移動(dòng)漸近線算法

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

移動(dòng)漸近線算法引入如下中間變量

式中Lj及Uj即所謂的移動(dòng)漸近線，隨著迭代而變化，且對(duì)第k次迭代來(lái)說(shuō)始終滿足算法對(duì)目標(biāo)函數(shù)和約束函數(shù)gi,i = 0,… ,l 在第k次迭代的迭代解 xk處的近似為

式中， r ik為殘值，且有

則優(yōu)化問(wèn)題（5）在第k次迭代處的移動(dòng)漸近線近似可寫(xiě)作

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

3 梯度控制算法

在異質(zhì)材料實(shí)體設(shè)計(jì)方法中，材料的體積分?jǐn)?shù)和物理屬性通過(guò)Q4單元的雙線性形函數(shù)實(shí)現(xiàn)了在二維設(shè)計(jì)空間中的0C 連續(xù)，避免了材料空間分布突變?cè)斐蓱?yīng)力集中[1]，為了提高該設(shè)計(jì)優(yōu)化方法的魯棒性和自主性，有必要使設(shè)計(jì)者可以直接控制異質(zhì)材料在空間中的梯度大小。Petersson等[19]將基于均質(zhì)單元和單元變量的鄰域單元梯度約束引入拓?fù)鋬?yōu)化，提供了一種直接控制局部材料梯度的方法，并且因此減少了數(shù)值計(jì)算中的棋盤(pán)模式和網(wǎng)格依賴現(xiàn)象，但這種算法大大增加了運(yùn)算量。本文將該局部材料梯度控制方法應(yīng)用于梯度單元和節(jié)點(diǎn)變量，并演化成一種新算法提高計(jì)算性能。

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

圖4 節(jié)點(diǎn)鄰域梯度控制示意圖

式中，xi和xj分別為節(jié)點(diǎn)i和節(jié)點(diǎn) j處的設(shè)計(jì)變量， 此式引入了個(gè)線性約束函數(shù)，直接計(jì)算不具備經(jīng)濟(jì)性，因此提出實(shí)施如下設(shè)計(jì)變量自適應(yīng)下界來(lái)實(shí)現(xiàn)式（16）的約束

4 優(yōu)化設(shè)計(jì)實(shí)例

本文以金屬夾鉗為例來(lái)驗(yàn)證說(shuō)明前述的異質(zhì)材料實(shí)體設(shè)計(jì)優(yōu)化方法。以鋼、鋁為兩相組成材料設(shè)計(jì)如圖5所示的夾鉗，設(shè)厚度t=1mm。設(shè)定鋼和鋁的楊氏模量分別為200GPa和68GPa，密度分別為7.85g/cm3和2.64g/cm3，兩者泊松比差別細(xì)微，簡(jiǎn)便起見(jiàn)均設(shè)為 0.3。將該夾鉗劃分為Q4單元的網(wǎng)格，以重量為目標(biāo)函數(shù)，A節(jié)點(diǎn)和B節(jié)點(diǎn)的Y方向位移為約束函數(shù)建立優(yōu)化數(shù)學(xué)模型如下

式中，W為夾鉗重量，uyA和uyB分別為A點(diǎn)和B點(diǎn)的Y方向位移，x為設(shè)計(jì)變量，設(shè)為各節(jié)點(diǎn)處鋼的體積分?jǐn)?shù)向量，初始化為，按前文所述方法進(jìn)行靈敏度分析、移動(dòng)漸近線優(yōu)化，得到的最優(yōu)材料分布如圖6(a)所示，與圖6(b)所示的由商業(yè)軟件 Optistruct計(jì)算得到的相同結(jié)構(gòu)的載荷路徑拓?fù)鋱D相比，可看到強(qiáng)度較高的鋼的分布與載荷路徑拓?fù)溆挟惽ぶ帯?/p>

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

圖6 夾鉗的梯度材料最優(yōu)分布與載荷路徑對(duì)比

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

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

5 結(jié) 論

本文在簡(jiǎn)明的梯度有限元概念基礎(chǔ)上提出的啟發(fā)式異質(zhì)材料實(shí)體設(shè)計(jì)優(yōu)化方法，流程清晰，運(yùn)算量合理，可以在二維空間得出實(shí)體材料的最優(yōu)分布，從而使設(shè)計(jì)目標(biāo)和設(shè)計(jì)約束達(dá)到最優(yōu)的均衡。這種設(shè)計(jì)方法避免了基于單元變量的材料分布優(yōu)化所導(dǎo)致的材料突變，具有材料梯度可控且控制方法運(yùn)算量小的特點(diǎn)。因?yàn)槔昧擞邢拊陨淼男魏瘮?shù)對(duì)材料分布進(jìn)行空間插值從而使得計(jì)算過(guò)程大大簡(jiǎn)化，并且同時(shí)達(dá)到了材料分布的C0連續(xù)。這種設(shè)計(jì)方法可以向三維空間擴(kuò)展并廣泛應(yīng)用到梯度功能材料的分布設(shè)計(jì)中去。該方法的缺點(diǎn)是只能設(shè)計(jì)兩相材料實(shí)體，如何處理兩相以上的材料的分布優(yōu)化成為了下一步工作的重點(diǎn)。

[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零件異質(zhì)材料建模方法[J]. 機(jī)械工程學(xué)報(bào),2004,40(5): 111-117.

[7] 寇欣宇,王以忠,彭一準(zhǔn). 基于非流形幾何與特征樹(shù)的異質(zhì)材料實(shí)體可視化方法[J]. 計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào),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.