含自重載荷的功能梯度材料結構時域動力學拓撲優化設計

文桂林 陳高錫 王洪鑫 薛 亮 魏 鵬 劉 杰，

1.燕山大學機械工程學院，秦皇島，066004

2.廣州大學機械與電氣工程學院，廣州，510006

3.湖南大學機械與運載工程學院，長沙，410082

4.華南理工大學土木與交通學院，廣州，510641

0 引言

具有良好力學性能的輕質結構設計近年來得到學術界和工程界越來越多的關注[1-3]，其中，功能梯度材料(functionally gradient materials，FGM)的力學性能可沿著一個或多個方向逐漸變化，具有相對于傳統均一材料更多的設計自由度和更優異的力學性能，近年來被廣泛應用于機械工程和航空航天等領域[1]。然而，傳統經驗設計方法很難最大程度地發揮FGM結構的最優機械性能。

結構拓撲優化方法通過尋求材料的最優分布可以設計出具有優良機械性能的FGM結構。NELLI等[4]首次提出了利用連續體結構拓撲優化來設計FGM結構，根據固體各向同性材料懲罰插值(solid isotropic material with penalization，SIMP)材料模型，提出一種新的材料模型FGM-SIMP來表示一個完全連續的FGM材料。在水平集方法框架下，XIA等[5]提出一種功能梯度結構拓撲優化設計新方法，實現了材料屬性和機械性能的同時優化。采用凸規劃求解策略以及周長控制方法，LI等[6]對具有拉脹特性的新型功能梯度蜂窩復合材料進行了拓撲優化設計。邱克鵬等[7]利用SIMP法實現了功能梯度MBB梁和功能梯度夾層結構夾芯的拓撲構型設計，并揭示了材料性能和材料插值模型對結構優化中材料分布的影響規律。李信卿等[8]對周期性功能梯度結構進行拓撲優化設計，獲得了具有較好散熱性能的FGM結構。雖然上述方法可以很好地提高FGM結構的機械性能，但都忽略了自重載荷和動力學特性，而面向重型機械裝置、船舶和航空航天裝備領域中FGM結構的設計時，自重載荷和動力學特性往往無法忽視。

已有研究對均一化材料組成的機械結構進行了考慮自重載荷的拓撲優化設計。CHEN等[9]提出了一種基于設計相關載荷的線彈性結構拓撲優化方法。BRUYNEEL等[10]研究了體積約束下柔順度最小化結構拓撲優化問題，發現當結構自重占主導地位時，會出現柔順度非單調性行為、最優結果無約束行為等，針對這些缺陷，張暉等[11]利用材料屬性的有理近似(rational approximation of material properties，RAMP)材料插值模型和平均敏度過濾技術很好地解決了該問題。還有學者利用雙線漸進結構優化方法和基于非均勻有理B樣條基函數插值的拓撲描述函數方法，成功解決了考慮自重載荷的拓撲優化設計問題[12-13]。然而，上述研究大多針對均一材料的結構拓撲優化設計，而FGM的材料密度不均勻分布時，自重載荷對拓撲優化設計結果影響更大。

此外，考慮動力學特性的拓撲優化設計方法在機械工程領域得到了廣泛的關注。OLIVER等[14]提出了一種增廣拉格朗日法，以解決一般動態載荷作用下具有應力約束的結構拓撲優化設計問題?；谝环N新的凝聚函數策略和時域求解結構動態響應的思路，ZHAO等[15]提出了一種新的拓撲優化方法，可以有效降低結構整個振動階段最大動態響應?；谝环N逐次迭代分析策略，KANG等[16]研究了考慮自然頻率的動態大規模拓撲優化問題。綜上研究發現，關于FGM結構的動力學拓撲優化設計研究較少。

基于此，本文提出了一種含自重載荷的功能梯度材料結構時域動力學拓撲優化方法。給出了針對FGM結構的自重載荷分配策略，研究了有/無自重載荷作用下，均一材料和功能梯度材料結構的優化構型差異性，并利用數值仿真和實驗方法研究了均一材料組成結構在考慮自重載荷和無自重載荷作用下優化結構的動力學性能。

1 功能梯度材料結構的材料模型

在SIMP材料插值模型框架下，建立相對單元密度和功能梯度材料彈性模量、功能梯度材料密度之間的關系[4，7]：

(1)

式中，EH、ρH分別為FGM的彈性模量和密度；E0、ρ0分別為初始設定的彈性模量和密度；a、b為FGM材料梯度變化的相關設計參數；x、y為設定在FGM區域的坐標方向值；α為小于1的正實數，通常令α=15/16[15]；η為單元密度即設計變量；p為懲罰系數，通常令p=3，該懲罰系數可以有效地消除中間密度單元并推動拓撲優化趨向離散0-1解。

當a和b等于0時，FGM沒有發生梯度變化，即相當于傳統的均一材料，為功能梯度材料的特例。

如圖1所示，FGM結構設計域離散后，單元彈性模量和單元密度分別為該單元中每個節點處彈性模量和密度值的平均值，即

圖1 二維功能梯度材料結構的材料分配策略

(2)

2 功能梯度材料結構自重載荷的分配策略

如圖2所示，考慮4節點的四邊形有限元網格和沿垂直y方向施加重力載荷，每個單元質量的四分之一分布在4個節點上，相鄰單元質量在共同節點上疊加。由于功能梯度材料結構的材料屬性呈梯度變化，所以功能梯度材料結構單元的自重載荷也呈梯度變化。因此，功能梯度材料結構受自重作用的單元載荷矢量為

圖2 二維功能梯度材料結構的自重載荷分配策略

(3)

式中，g為重力加速度；Vi為單元體積。

自重載荷相對于一個單元的變化可表示為

(4)

3 功能梯度材料結構動力學拓撲優化

3.1 拓撲優化模型

在上述材料插值模型和自重載荷分配策略的基礎上，以動柔度最小為優化目標，以結構整體的體積為約束，建立FGM結構拓撲優化列式：

(5)

式中，MH、ZH、KH分別為全局質量、阻尼和剛度矩陣；N為單元數量；Ve為單元e的體積；Vmax為期望的總材料體積；為了避免有限元方程在求解過程剛度矩陣的奇異性，單元最小密度不取零，而是小的正數，ηmin=0.001；F(t)包括了自重載荷和外載荷；u(t)為在t時刻的結構位移響應；tf為外載荷加載時長；η為設計變量向量。

目標函數表達為

(6)

式中，C為結構柔度。

3.2 靈敏度分析

(7)

外載荷與設計變量無關，將?f/?ηe表示為自重載荷對設計變量求導，進一步可得

(8)

式(7)中引入伴隨變量λ的微分方程為

(9)

(10)

s∈[0，tf]

Λ(s)=λ(tf-t)

求解出伴隨變量，式(8)可進一步寫為

(11)

設瑞利阻尼Z=αrM+βrK，且系數αr和βr相互獨立，則

(12)

(13)

(14)

體積約束對設計變量的導數為

(15)

求得目標函數和體積約束對設計變量的靈敏度信息后，利用移動漸近線方法(the method of moving asymptotes，MMA)[17]可以對式(5)中的優化問題進行求解。

4 算例

算例中使用的參數如下：彈性模量E0=200 GPa，泊松比υ=0.3，密度ρ=7800 kg/m3，重力加速度g=9.81 m/s2，結構自重G=體積×密度×重力加速度×體積分數。多項插值模型中的懲罰系數p=3，α=15/16。動態載荷的單元時間步長Δt=0.005 s。瑞利阻尼系數αr=10、βr=10-5，動柔度單位為N·m。為了便于描述，將無自重載荷和考慮自重載荷的動力學拓撲優化設計分別簡寫為W/oSL(without self-weight load)和WSL(with self-weight load)。

4.1 均一材料二維結構算例

圖3 二維兩端固定梁結構和外載荷

圖4所示為二維兩端固定梁結構施加靜態和動態載荷對優化結果的影響。其中，圖4a所示是考慮自重的靜力學拓撲優化結果；圖4b～圖4d所示分別為使外載荷激勵時間tf為0.1 s、0.03 s和0.01 s時的WSL結果。表1所示為二維兩端固定梁結構在加載時長tf=0.01 s時，不同外載荷(以與自重的比值量化)的優化結果，以及WSL和W/oSL的目標函數比值γ(γ=JWSL/JW/oSL)。

(a)靜力學，J=7.602 N·m (b)tf=0.1 s，J=2.879 N·m

表1 加載時長tf=0.01 s，不同外載荷對優化結果的影響

由圖4可知，動力學優化結果的目標函數值均小于靜力學優化結果，且當加載時長較大時(如tf=0.1 s)，動力學優化構型與靜力學優化構型近似，這與預期的一致；當加載時長較小時(如tf=0.03 s)，動力學優化構型明顯異于靜力學優化構型，加載時間越短(角頻率越大)，構型差異越明顯。由表1可知，受自重載荷的影響，優化結構的材料從固定端到外載荷施加點呈遞減分布，這是由于重力使簡支梁中間位置動態撓度變大，需要更多的材料分布到固定端以增加整體結構的動剛度。此外，當外載荷占自重比值逐漸變小時，動力學拓撲優化問題由自重載荷和外載荷共同作用變為自重載荷占主導，此時會有更多的材料向兩側固定端分布，且目標函數由2.291 N·m逐漸減小為0.339 N·m。此外，相同外載荷作用下，以WSL和W/oSL目標函數比值γ作為自重載荷對優化結構的影響程度指標，當外載荷占自重比值逐漸變小時，目標函數比值逐漸增大，與預期的結果一致。圖5是在加載時長tf=0.01 s、外載荷F0=100%G時WSL和W/oSL的目標函數和體積約束的迭代過程圖，可以發現，這兩種情況目標函數均收斂良好，體積分數均收斂至50%。

(a)WSL

4.2 功能梯度材料二維結構算例

初始設計域和邊界條件均與4.1節一致，進行功能梯度材料算例分析。表2所示為兩端固定二維梁結構在加載時間tf=0.01 s和外載荷F0=100%G時，不同類型功能梯度材料的WSL和W/oSL結果。

表2 不同二維功能梯度材料結構的優化結果對比

(1)當a=0、b=0時，二維功能梯度材料是均一材料。WSL和W/oSL的優化目標函數值分別為2.291 N·m和1.653 N·m，兩者目標函數比值為1.39。

(2)當a=0、b=0.06時，材料密度從下往上遞增且每層相同，FGM梯度方向與載荷施加方向一致。受FGM梯度方向的影響，優化最終構型的材料分布與均一材料有明顯不同，即有大量材料移動到結構下半部分的密度較低區域，以提供更大的結構動剛度。此外，FGM結構WSL和W/oSL的目標函數值分別為1.914 N·m和0.562 N·m，兩者目標函數比值為3.41，考慮自重情況下，有更多的材料向兩端約束位置分布，以抵抗自重載荷對FGM材料變形的影響。

(3)a=0、b=-0.06時，材料密度從上往下遞增且每層相同。材料梯度方向正好與a=0、b=0.06時的情況相反，FGM結構的最優材料分布規律也與之相反，即大量材料移動到結構上部分密度較低區域，以提高結構整體動剛度。WSL和W/oSL的目標函數值相差較小，分別為4.816 N·m和4.155 N·m，但是最終優化構型上有較明顯的差異，即為了抵抗自重載荷的影響，前者的材料更多地向兩側約束處移動。

(4)當a=0.02、b=0時，材料性能從左往右遞增且每列相同，梯度方向與載荷施加方向垂直，此時FGM結構產生撓度的方向和材料梯度方向垂直。對比WSL和W/oSL，重力載荷會增大FGM結構在垂直梯度方向的撓度，需要更多的材料分布到結構左側低密度區域以增加整體結構的動剛度；此外，由于FGM材料不對稱，最終的拓撲優化構型存在明顯的不對稱，且WSL的目標函數是W/oSL的8倍左右，說明了材料梯度方向對FGM結構動力學拓撲優化結果具有很大影響。

(5)當a=0.01、b=0.01時，材料性能從左往右遞增，梯度方向與載荷施加方向相交45°。此時，FGM結構變形方向與材料梯度方向相交45°，優化結構的材料部分向FGM結構左下角位置分布。同時，材料梯度分布不對稱導致最終的優化構型材料分布不對稱。

4.3 功能梯度材料三維結構算例

基于二維的FGM材料模型，三維FGM材料模型在z方向(厚度方向)梯度變化的相關設計參數設為q，則三維FGM的彈性模量和密度為

(16)

一個兩端固定三維梁結構，在其上表面中心位置施加矩形載荷F0=G，其幾何尺寸、邊界條件和外載荷如圖6所示。初始設計區域的長度、寬度和厚度分別為15 m、3 m和1 m，網格單元個數為150×30×10，最終設計體積為初始設計區域的40%，過濾半徑rmin=2.5，結構自重G=1 375 920 N。為了簡化問題且不失一般性，這里假設FGM材料模型在厚度方向不發生變化，即q=0。不過，對本文方法進行適當修改，可以很容易解決三個方向都有梯度變化的情況。

圖6 三維兩端固定梁結構和外載荷

圖7所示為自重載荷作用下，三維FGM結構在加載時間tf=0.1 s、外載荷F0=100%G情況下，不同類型的功能梯度材料的優化結果?？梢钥闯?，三維優化結構的邊界清晰，優化結果的規律與二維情況類似，表明該方法可以有效應用到三維功能梯度材料結構優化設計中，并進一步驗證所提方法的有效性和良好的工程應用潛力。

(a)WSL，F0=100%G，J=5.224 N·m

5 實驗與數值仿真研究

5.1 二維懸臂梁結構優化設計

圖8 懸臂梁結構初始設計域和邊界條件

(a)W/oSL (b)WSL

5.2 懸臂梁結構動力學數值仿真分析

對上述優化結構進行CAD幾何重建模，然后利用ABAQUS穩態動力學模塊進行動力學仿真分析。采用4節點殼單元(S4)對仿真結構進行網格劃分，全局單元網格尺寸為1 mm。圖10a和圖10b分別表示W/oSL和WSL前4階模態振型圖。注意到，懸臂梁優化結構在oxz平面的彎曲模態時，結構末端點在z方向的撓度最大，因此著重研究在oxz平面的彎曲模態。由圖10a可知，W/oSL優化結構的前4階固有頻率分別為129.41，504.33，633.39，715.39 Hz，其oxz平面彎曲模態是第3階模態，固有頻率為633 Hz；由圖10b可知，WSL優化結構的前4階固有頻率分別為152.94，558.25，728.56，768.95 Hz，其oxz平面彎曲模態是第4階模態，固有頻率為769 Hz?？梢?，WSL優化結構的每階固有頻率均有所提高，特別地，在工程中提高結構的第1階固有頻率可以在一定程度上避免結構發生共振。

(a)W/oSL

5.3 動力學測試實驗設置

動力學測試實驗系統如圖11所示，包括數據采集分析儀(型號AVANT-MI-6008)、功率放大器(型號HEA-200C)、激振器(型號JZK-10)、力傳感器(型號PCB-208C02)和加速度傳感器(型號KISTLER-8776A50M3)；實驗平臺采用歐標3030鋁型材進行搭建；實驗樣品采用三維打印機(型號JGMAKER-A6)制備，所用材料為PLA，打印層分辨力為0.1 mm，打印方向平行于結構的長度(圖12)。

圖11 實驗裝置

(a)W/oSL優化結構

5.4 實驗與數值仿真結果對比分析

由上述數值仿真結果可知，W/oSL優化結構的oxz平面彎曲模態為第3階模態，固有頻率為633 Hz；WSL優化結構的oxz平面彎曲模態為第4階模態，固有頻率為769 Hz。為了激發懸臂梁結構在oxz平面的彎曲共振，對W/oSL和WSL優化結構分別施加頻率為633 Hz和769 Hz、幅值均為0.42 N的正弦激勵。圖13和圖14所示為W/oSL和WSL優化結構的實驗和數值仿真對比結果。圖15從實驗上對比了WSL和W/oSL優化結構的位移響應。

(a)加速度-時間曲線

(a)加速度-時間曲線

圖15 WSL和W/oSL優化結構的位移-時間曲線實驗結果對比

可以發現，W/oSL優化結構的實驗和數值仿真的加速度幅值均為7.1g左右，而WSL優化結構均為5.3g左右。此外，兩種方法得到的位移幅值都十分相近(W/oSL為0.0039 mm和0.0045 mm；WSL為0.0021 mm和0.0022 mm)，表明實驗結果和仿真結果具有良好的一致性。通過實驗對比WSL和W/oSL優化結構的位移可以發現，WSL比W/oSL的位移幅值減小了近46%，說明本文所提方法得到的WSL結構具有更大的動剛度。這些結果均證實了所提方法的有效性。

6 結論

本文提出了一種含自重載荷的功能梯度材料結構時域動力學拓撲優化方法，系統研究了自重載荷和材料梯度方向對動力學拓撲優化設計結果的影響，并用實驗和數值方法驗證了所提方法的有效性。研究發現：①自重載荷和材料梯度分布形式對FGM結構的拓撲優化構型和動剛度具有較大的影響，且為了抵抗自重載荷對FGM結構變形的影響，拓撲優化構型的材料分布會向材料低密度區域移動；②考慮自重載荷的拓撲優化設計可以提高結構固有頻率和結構動剛度，在實際工程服役中具有更高的可靠度。需要指出的是，實際工程機械結構會受到不同的載荷激勵，針對不同的具體載荷激勵特征進行拓撲優化設計和分析十分有必要，筆者將在后續工作中進一步解決該問題。