考慮全局應(yīng)力約束的大變形柔順機構(gòu)拓?fù)錁?gòu)型設(shè)計

占金青 彭怡平 羅 震 劉 敏

(1.華東交通大學(xué)機電與車輛工程學(xué)院， 南昌 330013； 2.悉尼科技大學(xué)機械與機電工程學(xué)院， 悉尼 NSW2007)

0 引言

柔順機構(gòu)利用自身柔性構(gòu)件引起的變形來實現(xiàn)運動、力或能量傳遞或轉(zhuǎn)換，具有免裝配、無需潤滑和高精度等優(yōu)點[1-3]，因此廣泛應(yīng)用于精密定位、生物工程和微機電系統(tǒng)等領(lǐng)域[4-7]。

柔順機構(gòu)拓?fù)鋬?yōu)化設(shè)計在滿足約束條件下尋求機構(gòu)設(shè)計域內(nèi)材料的最佳分布，使機構(gòu)的某種性能達到最優(yōu)[8-9]。在實際工作中柔順機構(gòu)往往發(fā)生大變形或大轉(zhuǎn)動，因此采用非線性有限元理論進行大變形柔順機構(gòu)拓?fù)鋬?yōu)化設(shè)計非常有必要。PEDERSEN等[10]提出柔順機構(gòu)幾何非線性拓?fù)鋬?yōu)化設(shè)計方法，并且研究了考慮軌跡規(guī)劃的大變形柔順機構(gòu)設(shè)計問題。BRUNS等[11]采用弧長法進行非線性有限元求解，基于幾何非線性和材料非線性進行柔順機構(gòu)拓?fù)鋬?yōu)化設(shè)計。BRUNS等[12]提出考慮復(fù)雜跳躍性能的大變形柔順機構(gòu)拓?fù)鋬?yōu)化設(shè)計方法。DU等[13]采用無網(wǎng)格有限元方法進行大變形柔順機構(gòu)拓?fù)鋬?yōu)化設(shè)計，以避免低密度單元引起的不收斂問題。LUO等[14]采用參數(shù)化水平集方法進行大變形柔順機構(gòu)拓?fù)鋬?yōu)化設(shè)計，獲得具有光滑邊界的機構(gòu)構(gòu)型。LIU等[15]采用改進的超彈材料模型進行大變形柔順機構(gòu)拓?fù)鋬?yōu)化設(shè)計，能夠克服低密度單元的不穩(wěn)定性問題。ZHU等[16]提出基于完全解耦的大變形柔順機構(gòu)拓?fù)鋬?yōu)化設(shè)計方法，可抑制多自由度機構(gòu)的輸入、輸出耦合效應(yīng)。柔順機構(gòu)的大變形容易引起靜強度失效，因此需考慮應(yīng)力約束限制大變形柔順機構(gòu)最大的應(yīng)力。上述研究進行柔順機構(gòu)幾何非線性拓?fù)鋬?yōu)化設(shè)計時并未考慮應(yīng)力約束。

目前，柔順機構(gòu)應(yīng)力約束拓?fù)鋬?yōu)化研究主要基于小變形線彈性有限元理論。DUYSINX等[17]采用局部應(yīng)力約束進行柔順機構(gòu)拓?fù)鋬?yōu)化設(shè)計，但是局部約束數(shù)目太大，導(dǎo)致計算效率降低。LEON等[18]采用全局應(yīng)力約束進行柔順機構(gòu)拓?fù)鋬?yōu)化設(shè)計，有效地提高了計算效率。占金青等[19]采用P范數(shù)方法進行柔順機構(gòu)應(yīng)力約束拓?fù)鋬?yōu)化設(shè)計，能夠有效避免單節(jié)點鉸鏈。CHU等[20]采用水平集方法進行多相材料柔順機構(gòu)應(yīng)力約束拓?fù)鋬?yōu)化設(shè)計，使每相材料結(jié)構(gòu)滿足相應(yīng)的應(yīng)力約束。

為了避免大變形引起靜強度失效，本文提出一種考慮全局應(yīng)力約束的大變形柔順機構(gòu)拓?fù)鋬?yōu)化設(shè)計方法。采用Total-Lagrangian描述方法求解機構(gòu)的大變形響應(yīng)問題；采用基于假設(shè)密度場的單元勢能插值法，以避免幾何非線性拓?fù)鋬?yōu)化出現(xiàn)的數(shù)值不收斂問題；以輸出位移最大化作為優(yōu)化目標(biāo)，采用改進的P范數(shù)將所有單元局部應(yīng)力約束轉(zhuǎn)換為一個全局應(yīng)力約束，構(gòu)建考慮全局應(yīng)力約束的柔順機構(gòu)幾何非線性拓?fù)鋬?yōu)化模型；采用移動漸近算法求解大變形柔順機構(gòu)應(yīng)力約束拓?fù)鋬?yōu)化問題。

1 幾何非線性有限元分析

采用Total-Lagrangian描述方法，以初始構(gòu)型作為參考構(gòu)型，大變形結(jié)構(gòu)的任一點i的構(gòu)型與位移之間關(guān)系表示為

xi=ui+Xi

(1)

式中Xi——初始構(gòu)型

ui——任一點的位移矢量

xi——現(xiàn)時構(gòu)型

對式(1)關(guān)于初始構(gòu)型Xi求導(dǎo)，可求導(dǎo)變形梯度F為

(2)

式中δij——Kronecker符號函數(shù)

變形梯度F左乘其轉(zhuǎn)置矩陣可得右Cauchy-Green變形C為

C=FTF

(3)

當(dāng)僅為剛體轉(zhuǎn)動時，Cij=δij。

變形梯度F的行列式稱為Jacobian行列式J，表示結(jié)構(gòu)體積的變化，其表示為

J=det(F)

(4)

由右Cauchy-Green變形C可得初始構(gòu)型Xi下的應(yīng)變，即為Green-Lagrange應(yīng)變E，計算式為

(5)

式中I——單位矩陣

與Green-Lagrange應(yīng)變E能量共軛的應(yīng)力是第二類Piola-Kirchhoff應(yīng)力S，可表示為

(6)

式中Φ——彈性勢能密度函數(shù)，本文采用St. Venant-Kirchhoff模型

第二類Piola-Kirchhoff應(yīng)力S對Green-Lagrange應(yīng)變E求導(dǎo)，可得彈性矩陣D為

(7)

將機構(gòu)設(shè)計域離散為有限個單元，大變形結(jié)構(gòu)的平衡方程可表示為

(8)

式中R——結(jié)構(gòu)殘余力

Fin——外載荷向量

BL——應(yīng)變矩陣V0——結(jié)構(gòu)體積

采用增量式Newton-Raphson法求解幾何非線性有限元平衡方程

KTΔu=R

(9)

其中

式中 Δu——節(jié)點位移矢量的增量

KT——切線剛度矩陣

結(jié)構(gòu)的節(jié)點位移矢量可更新為

un+1=un+Δu

(10)

第二類Piola-Kirchhoff應(yīng)力S也可表示為

S=JF-1σ(F-1)T

(11)

式中σ——Cauchy應(yīng)力，它為現(xiàn)時構(gòu)型xi下的應(yīng)力，表示大變形結(jié)構(gòu)的真實應(yīng)力

由式(11)可求得Cauchy應(yīng)力σ為

σ=J-1FSFT

(12)

(13)

2 單元勢能插值法

對于大變形柔順機構(gòu)拓?fù)鋬?yōu)化問題，低單元密度區(qū)域的網(wǎng)格容易發(fā)生畸變，從而導(dǎo)致數(shù)值不收斂現(xiàn)象。采用基于假設(shè)密度場的單元勢能插值法[21]克服這種不收斂現(xiàn)象，對實體單元和空洞單元分別采用非線性理論和線性理論求解；任一單元e的勢能函數(shù)Φe表示為

Φe(ue)=(Φ(γeue)-ΦL(γeue)+ΦL(ue))Ee

(14)

式中Ee——單元e彈性模量

ue——單元節(jié)點位移矢量

Φ——單位彈性模量對應(yīng)的勢能密度函數(shù)

ΦL——基于線性理論的勢能密度函數(shù)

γe——插值系數(shù)

對于空洞單元，采用線性理論，即插值系數(shù)γe=0；對于實體單元，采用非線性理論，即插值系數(shù)γe=1。

采用改進的固體各向同性材料懲罰模型[22]描述機構(gòu)的材料分布，單元彈性模量Ee表示為

(15)

式中E1——材料彈性模量

E0——空洞材料彈性模量，取10-9E1

k——懲罰系數(shù)，取3

空洞單元和實體單元插值系數(shù)分別為0和1，但是拓?fù)鋬?yōu)化設(shè)計采用連續(xù)單元密度變量，因此存在介于實體單元和空洞域單元之間的中間密度單元，需要對插值系數(shù)進行插值，使其在取值區(qū)間[0,1]連續(xù)變化。因此，采用閾值映射函數(shù)[11]對插值系數(shù)γe進行插值處理，即

(16)

式中ρ0——界定為實體單元的閾值

β1——控制參數(shù)

3 大變形柔順機構(gòu)全局應(yīng)力約束拓?fù)鋬?yōu)化模型

3.1 拓?fù)鋬?yōu)化模型

為了滿足機構(gòu)的運動需求，以大變形柔順機構(gòu)的輸出位移最大化為目標(biāo)函數(shù)；為了避免大變形引起的靜強度失效，采用單元勢能插值法，建立考慮全局應(yīng)力約束的柔順機構(gòu)幾何非線性拓?fù)鋬?yōu)化數(shù)學(xué)模型為

(17)

式中ρ——單元密度矢量

Uout——柔順機構(gòu)輸出位移

V——優(yōu)化后體積

ve——單元e體積

ρe——單元e密度

f*——允許體積份數(shù)

σmax——機構(gòu)最大von Mises應(yīng)力

σ*——允許應(yīng)力極限值

N——有限單元數(shù)目

3.2 全局應(yīng)力約束

為了避免大變形柔順機構(gòu)應(yīng)力約束拓?fù)鋬?yōu)化出現(xiàn)奇異解現(xiàn)象，采用應(yīng)力松弛方法[23]對所有單元的von Mises應(yīng)力進行懲罰

(18)

式中σe——松弛應(yīng)力

q——松弛因子，取值為1

同時，為了提高計算效率，采用改進的P范數(shù)方法[24]將柔順機構(gòu)所有的單元局部應(yīng)力凝聚化為一個全局的P范數(shù)應(yīng)力，即

(19)

式中p——P范數(shù)參數(shù)

理論上，P范數(shù)參數(shù)p→∞時，P范數(shù)應(yīng)力趨近最大應(yīng)力σmax。然而，p越大，大變形柔順機構(gòu)應(yīng)力約束拓?fù)鋬?yōu)化問題的高度非線性程度越明顯，甚至導(dǎo)致迭代不收斂。p的取值不能過大，否則導(dǎo)致P范數(shù)應(yīng)力和機構(gòu)的最大應(yīng)力存在一定的數(shù)值差異。為了消除這種數(shù)值差異，采用自適應(yīng)縮放方法[25]修正機構(gòu)的P范數(shù)應(yīng)力，即

σmax≈cσPN

(20)

式中c——自適應(yīng)系數(shù)

4 靈敏度分析及過濾技術(shù)

4.1 靈敏度分析

柔順機構(gòu)輸出位移Uout可表示為

Uout=LTu

(21)

式中L——伴隨載荷矢量，除了輸出位移方向?qū)?yīng)的載荷為1，其他自由度方向的載荷均為0

引入伴隨矩陣λ1，機構(gòu)輸出位移Uout可表示為

(22)

由式(22)可求得機構(gòu)輸出位移的靈敏度為

(23)

(24)

因此，機構(gòu)輸出位移Uout靈敏度為

(25)

同理，引入伴隨矩陣λ2，由式(18)～(20)可求得應(yīng)力約束的靈敏度為

(26)

(27)

式(26)可簡寫為

(28)

由式(18)～(20)可求得σmax對位移u的導(dǎo)數(shù)為

(29)

其中

σmax對位移ue的導(dǎo)數(shù)可改寫為

(30)

其中

A=(-J-1C-1BL)TFSFT

式中G——輔助矩陣

通過應(yīng)力約束對單元物質(zhì)密度求偏導(dǎo)可得

(31)

由式(26)和式(32)可得應(yīng)力約束的靈敏度為

(32)

4.2 映射過濾方法

基于變密度法的大變形柔順機構(gòu)拓?fù)鋬?yōu)化設(shè)計結(jié)果容易出現(xiàn)棋盤格、網(wǎng)格依賴性等不穩(wěn)定性現(xiàn)象。采用密度過濾方法[26]避免上述數(shù)值不穩(wěn)定性現(xiàn)象，即

(33)

其中

w(yj)=rmin-‖yj-ye‖

yj——單元j中心坐標(biāo)

vj——單元j體積ρj——單元j密度

Ne——與單元e中心的距離不超過最小過濾半徑rmin的所有單元集合

w(yj)——權(quán)重因子

ye——單元e中心坐標(biāo)

采用密度過濾獲得的拓?fù)錁?gòu)型邊界仍然出現(xiàn)較多的中間單元。采用光滑的Heaviside映射函數(shù)[27]修正過濾的單元密度，使得單元密度向0和1兩端集中，即

(34)

式中η——閾值，取值為0.5

β——控制光滑程度的參數(shù)

對于大變形柔順機構(gòu)應(yīng)力約束拓?fù)鋬?yōu)化問題，采用移動漸近線算法(Method of moving asymptotes，MMA)[28]進行求解。綜上所述，本文提出的考慮全局應(yīng)力約束的大變形柔順機構(gòu)拓?fù)鋬?yōu)化設(shè)計方法流程圖如圖1所示。

5 數(shù)值算例

圖2表示柔性反向器的設(shè)計域、固定邊界、作用載荷和輸出作用。設(shè)計域L×L為300 μm×300 μm，厚度d為7 μm，左上端和左下端為固定邊界，載荷Fin作用在左端中點，其大小為200 mN，輸入彈簧剛度kin和輸出彈簧剛度kout分別為3 000、300 N/m。選用多晶硅材料，其彈性模量E1為180 GPa，泊松比μ為0.3，強度極限范圍一般為1～10 GPa[29]。允許體積份數(shù)f*為0.25，控制參數(shù)β初始值取為1，每間隔50迭代步增大1倍，最大值為14，P范數(shù)參數(shù)p取值為8，過濾半徑rmin為4.5倍單元尺寸。輸入、輸出作用區(qū)域的應(yīng)力不作限制，如圖2黑色區(qū)域所示。

首先，采用無應(yīng)力約束的柔順機構(gòu)幾何非線性拓?fù)鋬?yōu)化模型進行柔性反向器設(shè)計，拓?fù)鋬?yōu)化結(jié)果如圖3所示，反向器構(gòu)型出現(xiàn)集中式柔度的單節(jié)點鉸鏈，最大應(yīng)力為5.70 GPa，應(yīng)力分布極為不均勻，輸出位移為-6.34 μm。允許應(yīng)力約束值取為1.50 GPa，進行考慮全局應(yīng)力約束的大變形反向器幾何非線性拓?fù)鋬?yōu)化設(shè)計，拓?fù)鋬?yōu)化結(jié)果如圖4所示，反向器構(gòu)型中的鉸鏈區(qū)域被擴大，使得機構(gòu)的柔度分布更加均勻，從而有效地抑制了集中式柔度的單節(jié)點鉸鏈出現(xiàn)，且機構(gòu)的應(yīng)力分布更加均勻，但應(yīng)力約束導(dǎo)致輸出位移有所減小。

圖5為考慮全局應(yīng)力約束的大變形反向拓?fù)鋬?yōu)化迭代過程曲線，優(yōu)化初始階段機構(gòu)的最大應(yīng)力急劇下降，隨后又急劇增大，這是由材料突然被大量移除而導(dǎo)致；由于參數(shù)β值每間隔50迭代步發(fā)生變化，引起最大應(yīng)力間隔50迭代步會發(fā)生振蕩；最終很好地收斂于允許的應(yīng)力約束值1.50 GPa。

考慮有、無全局應(yīng)力約束拓?fù)鋬?yōu)化獲得的大變形柔性反向器的變形圖如圖6所示。與無應(yīng)力約束優(yōu)化結(jié)果相比，應(yīng)力約束(σ*為1.50 GPa)拓?fù)鋬?yōu)化獲得的大變形反向器構(gòu)型中的鉸鏈區(qū)域被擴大，使得發(fā)生變形的區(qū)域相對更大，而不是依靠集中式柔度的單節(jié)點鉸鏈實現(xiàn)大變形。

允許應(yīng)力約束值σ*分別為3.00、2.75、2.50、2.25、2.00、1.50 GPa進行反向器拓?fù)鋬?yōu)化設(shè)計，獲得的優(yōu)化結(jié)果如圖7所示。隨著允許的應(yīng)力約束值減小，反向器構(gòu)型中的鉸鏈區(qū)域逐漸被擴大，使得機構(gòu)的柔度分布更加均勻，反向器構(gòu)型逐漸由集中式柔度分布轉(zhuǎn)變?yōu)榉植际饺岫确植?，鉸鏈區(qū)域的應(yīng)力分布更加均勻，最大應(yīng)力均很好地滿足應(yīng)力約束，表明應(yīng)力約束的大變形柔順機構(gòu)拓?fù)鋬?yōu)化方法能夠有效地抑制集中式柔度的單節(jié)點鉸鏈；但是機構(gòu)的輸出位移逐漸減小，如表1所示。

表1 不同允許應(yīng)力獲得的反向器拓?fù)鋬?yōu)化結(jié)果Tab.1 Results of inverters obtained by topology optimization using different allowable stress limits

6 結(jié)論

(1)為了避免大變形引起靜強度失效，采用單元勢能插值法以避免幾何非線性拓?fù)鋬?yōu)化出現(xiàn)的數(shù)值不收斂問題，利用改進的P范數(shù)引入全局的應(yīng)力約束，提出一種考慮全局應(yīng)力約束的大變形柔順機構(gòu)拓?fù)鋬?yōu)化設(shè)計方法。

(2)與無應(yīng)力約束的大變形柔順機構(gòu)拓?fù)鋬?yōu)化結(jié)果相比，采用全局應(yīng)力約束的大變形柔順機構(gòu)拓?fù)鋬?yōu)化方法獲得的機構(gòu)構(gòu)型鉸鏈區(qū)域被擴大，使柔度分布更加均勻，抑制了集中式柔度單節(jié)點鉸鏈的出現(xiàn)，應(yīng)力分布更加均勻。

(3)隨著允許應(yīng)力約束值的減小，大變形柔順機構(gòu)構(gòu)型中的鉸鏈區(qū)域逐漸擴大，使機構(gòu)的柔度分布更加均勻，鉸鏈區(qū)域的應(yīng)力分布更加均勻，并且機構(gòu)的最大應(yīng)力均很好地滿足應(yīng)力約束，但是輸出位移逐漸減小。