李永欣， 賀星星， 李凱倫， 楊權(quán)威

(1.燕山大學(xué) 機械工程學(xué)院，河北省重型智能制造裝備技術(shù)創(chuàng)新中心，秦皇島 066044；2.天津恒天新能源汽車研究院有限公司，天津 300450)

1 引 言

隨著汽車、船舶和航空航天等領(lǐng)域?qū)Y(jié)構(gòu)輕量化和結(jié)構(gòu)性能需求的提高，通過拓?fù)鋬?yōu)化在概念設(shè)計階段獲得顛覆性設(shè)計獲得學(xué)者和設(shè)計人員的青睞。目前，常用的拓?fù)鋬?yōu)化方法有均勻化方法、變密度法、水平集法、獨立連續(xù)映射法和漸進結(jié)構(gòu)優(yōu)化法，其中，變密度法的使用范圍最為廣泛，且已經(jīng)引入到眾多商業(yè)軟件中。但是變密度法對于動態(tài)響應(yīng)尤其是非線性動態(tài)處理上適用性較差。

動態(tài)結(jié)構(gòu)優(yōu)化方法主要有近似代理模型法[1]、慣性釋放法IRM(Inertia Relief Method)[2]、混合元胞自動機法HCA(Hybrid Cellar Automata)[3]和等效靜態(tài)載荷法ESLM(Equivalent Static Loads Method)。等效靜態(tài)載荷法通過將動態(tài)優(yōu)化過程轉(zhuǎn)化為靜態(tài)優(yōu)化過程，提升了計算效率，并通過多重循環(huán)提高了優(yōu)化的穩(wěn)定性，且優(yōu)化目標(biāo)更為多樣，更符合工程實際需求，ESLM法已經(jīng)應(yīng)用于ESLDyna，OPTISTRUCT和NASTRAN中。

等效靜態(tài)載荷法是由Kang等[4]首先提出，并將之運用于線性動態(tài)系統(tǒng)的尺寸優(yōu)化和形狀優(yōu)化。隨后，Park等[5,6]證明了基于ESLM的優(yōu)化方法滿足Karush-Kuhn-Tucker條件，并且在數(shù)學(xué)意義上，ESLM和直接動態(tài)優(yōu)化方法是等價的。Jang等[7]將ESLM推廣到了線性動態(tài)拓?fù)鋬?yōu)化中。隨后，學(xué)者們在ESLM的計算效率[8-10]、周期性微結(jié)構(gòu)[11]、變幅載荷疲勞約束[12]、可靠性[13]和多目標(biāo)[14,15]動態(tài)優(yōu)化設(shè)計等領(lǐng)域進行拓展，并將其應(yīng)用于前端結(jié)構(gòu)[16]和車身[17]等汽車結(jié)構(gòu)的優(yōu)化設(shè)計中。文獻[18～20]在提高ESLM的等效精度方面進行了一些嘗試。

經(jīng)典的等效靜態(tài)載荷法遵循位移等效原則，即等效靜態(tài)載荷能夠產(chǎn)生和動態(tài)分析完全一致的位移響應(yīng)。當(dāng)結(jié)構(gòu)處于彈性范圍內(nèi)時，等效靜態(tài)載荷的計算是準(zhǔn)確的，這主要是因為結(jié)構(gòu)發(fā)生彈性變形后可以完全恢復(fù)，所以下一時刻的結(jié)構(gòu)變形和上一時刻完全無關(guān)?，F(xiàn)有的等效靜態(tài)載荷法能夠很好地處理線彈性材料的動態(tài)拓?fù)鋬?yōu)化問題，但該方法在計算等效靜態(tài)載荷時完全忽略了材料的非線性，在此方面，Lee等[21]試圖考慮材料的非線性因素，但是在計算等效靜態(tài)載荷時仍采用線剛度矩陣。在結(jié)構(gòu)發(fā)生塑性變形后，由于塑性變形不可恢復(fù)，發(fā)生塑性變形后會一直影響以后時刻的變形，因此結(jié)構(gòu)剛度具有時效性。而當(dāng)結(jié)構(gòu)應(yīng)力超過屈服極限時，若仍采用線剛度矩陣計算等效靜態(tài)載荷，由于位移等效原則勢必導(dǎo)致計算載荷偏大。

另外，根據(jù)有限元理論，剛度矩陣是按照節(jié)點編號順序和x,y，z三個方向的自由度排列的，載荷向量也是按照節(jié)點編號順序和x,y，z三個方向排列。由于等效靜態(tài)載荷法會將作用于某點的載荷轉(zhuǎn)化為x,y和z三個方向的分力，結(jié)構(gòu)在三個方向上的塑性變形不同會導(dǎo)致結(jié)構(gòu)剛度在三個方向的變化也不同，這就意味著某時刻等效后三個方向的分力形成的合力與真實合力在方向上存在誤差。所以，如果用線彈性矩陣進行塑性變形中的載荷等效，不僅載荷數(shù)值偏大，而且也會造成載荷方向和真實的載荷方向的偏差。

綜上所述，ESLM無法處理塑性變形材料非線性導(dǎo)致的剛度時效和方向異性。其難點在于如何合理地將塑性變形對結(jié)構(gòu)的影響反映到等效靜態(tài)載荷的計算公式中，使之即能保證等效精度，又不會過多地增加計算成本。針對上述問題，本文基于有限單元法，推導(dǎo)了不同時間點楊氏模量的數(shù)學(xué)表達式和修正的等效靜態(tài)載荷方程，研究了塑性變形對等效靜態(tài)載荷法的影響和能有效解決彈塑性結(jié)構(gòu)動態(tài)拓?fù)鋬?yōu)化問題的等效靜態(tài)載荷法，以期為材料非線性的結(jié)構(gòu)動態(tài)響應(yīng)優(yōu)化提供有效的解決方法。

2 等效靜態(tài)載荷法的塑性修正

2.1 模量比率因子

ESLM不適于非線性材料的主要原因是發(fā)生塑性變形以后的結(jié)構(gòu)剛度發(fā)生改變。所以，通過變形前后結(jié)構(gòu)楊氏模量的改變情況來對等效靜態(tài)載荷進行修正是合理的。由此，提出模量比率因子MRC(Modulus Ratio Coefficient)的概念，模量比率因子的含義是某時刻結(jié)構(gòu)楊氏模量與初始時刻結(jié)構(gòu)楊氏模量的相對變化。它反映了材料的非線性特性對系統(tǒng)響應(yīng)隨時間的影響。

在線彈性階段，結(jié)構(gòu)的應(yīng)力-應(yīng)變的比值是一個定值。受力超出屈服極限后，結(jié)構(gòu)的總應(yīng)變?yōu)閺椥詰?yīng)變和塑性應(yīng)變之和，應(yīng)力-總應(yīng)變的比值將不再是一個確定的數(shù)，而是一個和應(yīng)力應(yīng)變相關(guān)的值。某時刻三維結(jié)構(gòu)沿x，y和z方向的楊氏模量分別表示為

(1)

可以看出，計算楊氏模量能夠在數(shù)學(xué)形式上真實反映某時刻結(jié)構(gòu)的應(yīng)力應(yīng)變關(guān)系。對于發(fā)生塑性變形的結(jié)構(gòu)，下一時刻的結(jié)構(gòu)變形永遠(yuǎn)受上一時刻的影響。因此，將動力學(xué)分析中的整個時間歷程等分成若干份，計算每個時刻的計算楊氏模量，并計算相鄰兩時刻的計算楊氏模量的比值βj i(t)，可表示為

(i=x,y,z;t≥1;j=1,2,…,n)(2)

(i=x,y,z;t=0,2,4,…,T)(3)

式中λi(T)為T時刻i方向上的的模量比率因子，n為T時刻結(jié)構(gòu)中塑性變形不為0的節(jié)點數(shù)量，βj i(t)=1為第j個節(jié)點在t時刻與t-1時刻i方向上的計算楊氏模量比值。

2.2 基于模量比率因子的等效靜態(tài)載荷法

在x,y和z方向上都應(yīng)計算模量比率因子，得到等效靜態(tài)載荷的計算修正方程為

(i=x,y,z)(4)

基于等效靜態(tài)載荷法的結(jié)構(gòu)優(yōu)化流程如圖1所示。由靜態(tài)響應(yīng)優(yōu)化的內(nèi)層循環(huán)和設(shè)計變量更新的外層循環(huán)組成，靜態(tài)響應(yīng)優(yōu)化由靜態(tài)優(yōu)化模型保證，將優(yōu)化的結(jié)構(gòu)變量輸出給外循環(huán)，判斷是否滿足約束條件。外層循環(huán)進行動態(tài)響應(yīng)和等效載荷的計算。

圖1 基于等效靜態(tài)載荷法的結(jié)構(gòu)優(yōu)化循環(huán)過程

以柔度為目標(biāo)函數(shù)，MRC-ESLM的變密度拓?fù)鋬?yōu)化列式為

基于ESLM的變密度拓?fù)鋬?yōu)化，由于拓?fù)鋬?yōu)化的設(shè)計變量為每個單元的密度，其收斂準(zhǔn)則比較復(fù)雜，可表示為

(6)

3 算例分析

3.1 算例1

為了體現(xiàn)在結(jié)構(gòu)發(fā)生塑性變形時MRC-ESLM 的效果，采用link1桿單元建立如圖2所示的桁架結(jié)構(gòu)有限元模型，桁架的彈性模量為207 GPa，泊松比為0.3，密度為7850 kg/m3，屈服極限為235 MPa，剪切模量為79 GPa。在節(jié)點2和節(jié)點3上施加沿y軸負(fù)向的周期為1 s幅值為150 kN的正弦載荷，加載時長為0.5 s。進行非線性分析時，設(shè)置載荷步數(shù)為200，子步數(shù)為50。利用式(4)計算0.25 s修正后的等效靜態(tài)載荷。

圖2 桁架結(jié)構(gòu)有限元模型

利用計算得到的等效靜態(tài)載荷進行靜態(tài)分析，提取各節(jié)點的位移，并與瞬態(tài)分析中各節(jié)點的位移做對比，如圖3所示，桁架的所有節(jié)點都產(chǎn)生了塑性變形。可以明顯看到，在x方向上，節(jié)點2、節(jié)點3和節(jié)點6的位移很小，幾乎為0，而節(jié)點7和節(jié)點5的位移比較大。在y方向上，各個節(jié)點都有較大位移，出現(xiàn)這種情況主要是因為結(jié)構(gòu)約束及外載荷的對稱性。所以本文主要分析節(jié)點5和節(jié)點7的x方向位移，以及各個節(jié)點在y方向上的位移。

進一步由圖3可知，結(jié)構(gòu)出現(xiàn)塑性變形時，瞬態(tài)分析和靜態(tài)分析的節(jié)點位移有較大差別，得到的系統(tǒng)響應(yīng)差別較大，等效效果不好，尤其是y方向的位移響應(yīng)。而本算例中，施加在桁架上的載荷的方向是y方向，主要產(chǎn)生的位移是y方向位移，進一步說明了原ESLM法的局限性。而利用修正的等效靜態(tài)載荷法得到的位移響應(yīng)和瞬態(tài)分析幾乎相同。

圖3 動靜態(tài)分析結(jié)果對比

3.2 算例2

懸臂梁結(jié)構(gòu)有限元模型如圖4所示，懸臂梁的尺寸為100 mm×200 mm，網(wǎng)格尺寸為5 mm。完全約束板左邊的位移，在板的右邊中間的點上施加半正弦載荷，載荷的周期和加載時長均為 0.5 s。懸臂梁的彈性模量為207 GPa，泊松比為0.3，密度為7850 kg/m3，屈服極限為235 MPa，剪切模量為 79 GPa。對右邊中間節(jié)點施加載荷幅值為1300 kN的半正弦載荷。然后對該懸臂板進行動力學(xué)分析，提取節(jié)點位移和結(jié)構(gòu)剛度矩陣。為了比較修正前后的ESLM在整個時間歷程上的修正效果，分別提取0.2 s,0.25 s,0.3 s,0.35 s和 0.4 s 的節(jié)點位移。利用式(4)計算上述時間點修正后的等效靜態(tài)載荷，列入表1。

圖4 懸臂板有限元模型

利用表1的等效靜態(tài)載荷值進行有限元分析。根據(jù)有限元分析結(jié)果，當(dāng)載荷幅值為1300 kN、結(jié)構(gòu)發(fā)生塑性變形時，為了比較修正前后對所有節(jié)點位移的影響，分別計算ESLM和MRC-ESLM得到的位移響應(yīng)和動態(tài)響應(yīng)之間所有節(jié)點的誤差均值，計算結(jié)果如圖5所示。

表1 等效靜態(tài)載荷Tab.1 Equivalent static loads

圖5 修正前后所有節(jié)點位移誤差均值對比

可以看出，在x方向上，隨著時間的推移，雖然兩種方法得到的位移與動態(tài)位移之間的誤差都在逐漸增大，但是MRC-ESLM的誤差均值要小于ESLM的誤差均值，ESLM的最大均值誤差達到了 8.930%，而MRC-ESLM的最大均值誤差為 5.197%。

由于本算例大部分節(jié)點沒有發(fā)生塑性變形，因此選取節(jié)點2來說明模量比率因子的作用，如圖6所示?？梢婋m然節(jié)點2的塑性應(yīng)變值為0，但是由于受到結(jié)構(gòu)塑性變形的影響，ESLM 得到的位移響應(yīng)仍和動態(tài)分析得到的位移響應(yīng)有較大差別。相較于ESLM，不管節(jié)點塑性應(yīng)變值是否為0，MRC-ESLM得到的位移響應(yīng)都更接近動態(tài)分析得到的位移響應(yīng)，等效精度有了很大提高，證明所采用的模量比率因子有較好的修正效果。

圖6 節(jié)點2的位移

3.3 汽車膝緩沖器拓?fù)鋬?yōu)化

汽車膝緩沖器的主要作用是在汽車發(fā)生碰撞時盡可能地吸收能量，減少碰撞對乘員膝部的傷害。汽車膝緩沖器的幾何模型如圖7所示。設(shè)計域高120 mm，彈性模量E=70 GPa，泊松比為0.33，密度為2700 kg/m3，屈服強度為240 MPa，剪切模量為20.7 GPa。剛體的質(zhì)量為270 kg，位于膝緩沖器正上方10 mm處，初始速度為1.25 m/s。限制膝緩沖器的最大網(wǎng)格尺寸為2.5 mm，并保證總體沙漏不超過模型總體內(nèi)能的5%，對膝緩沖器進行非線性動力學(xué)分析可得，碰撞完成后膝緩沖器產(chǎn)生的等效塑性應(yīng)變?yōu)?.295%。

圖7 汽車膝緩沖器幾何模型

動態(tài)工況轉(zhuǎn)化為了21個靜態(tài)工況，其中第一個工況為即將接觸時刻的工況，此時膝緩沖器未發(fā)生變形。分別采用ESLM和MRC-ESLM對膝緩沖器進行拓?fù)鋬?yōu)化，優(yōu)化過程中以柔度為目標(biāo)函數(shù)，體積比約束為0.2，最大侵入量約束為2 mm。兩種方法的非線性動態(tài)分析次數(shù)分別為16次和15次，最后一次迭代時，第2-21個工況的模量比率因子如圖8所示。

圖8 各工況模量比率因子

可以看出，對于同一工況(即同一時刻點)，x，y和z方向的模量比率因子區(qū)別很大，膝緩沖器受力的主要方向,即z方向的模量比率因子一直遠(yuǎn)小于x和y方向，這意味著在受力主要方向上的修正幅度最大，而且修正后作用于節(jié)點的等效靜態(tài)載荷的方向也不同。而對于同一方向的不同工況來說，隨著時間的推移，整體呈現(xiàn)先下降后略微上升的趨勢，這表明隨著結(jié)構(gòu)塑性變形的增大，等效靜態(tài)載荷的修正幅度也在增大。同時，從圖8的最后一個工況在修正前后的等效靜態(tài)載荷可以看出，由于修正時各方向的模量比率因子不同，導(dǎo)致修正后三個方向的等效靜態(tài)載荷形成的合力與修正前在大小上和方向上均存在較大差異。由于在內(nèi)循環(huán)，即靜態(tài)拓?fù)鋬?yōu)化時選擇以柔度為目標(biāo)函數(shù)，所以該合力在方向上的修正會對靜態(tài)優(yōu)化結(jié)果產(chǎn)生較大影響，進而影響最終的優(yōu)化結(jié)果。優(yōu)化后模型的等效塑性應(yīng)變結(jié)果如圖9所示。

圖9 拓?fù)鋬?yōu)化結(jié)果

可以看出，ESLM法優(yōu)化結(jié)果的等效塑性應(yīng)變比原結(jié)果稍大，而采用MRC-ESLM法得到的優(yōu)化結(jié)果在減輕結(jié)構(gòu)質(zhì)量的基礎(chǔ)上還減小了塑性應(yīng)變。

4 結(jié) 論

(1) 針對塑性變形的動態(tài)響應(yīng)過程中結(jié)構(gòu)應(yīng)力應(yīng)變關(guān)系的時效性和方向異性，提出了一種基于模量比率因子的改進等效靜態(tài)載荷法。模量比率因子的含義是某時刻結(jié)構(gòu)楊氏模量與初始時刻結(jié)構(gòu)楊氏模量的相對變化比值，反映了材料的非線性特性對系統(tǒng)響應(yīng)的影響。

(2) 算例表明，通過計算動靜轉(zhuǎn)化時刻各個方向上的模量比率因子對等效靜態(tài)載荷計算公式進行修正，在產(chǎn)生塑性變形的結(jié)構(gòu)中，無論節(jié)點是否發(fā)生塑性變形，模量比率因子都可以達到很好的等效精度。

(3) 改進的等效靜態(tài)載荷法能夠穩(wěn)定有效地提高等效精度，降低彈塑性結(jié)構(gòu)動態(tài)拓?fù)鋬?yōu)化問題的計算成本，為動態(tài)拓?fù)鋬?yōu)化設(shè)計問題提供可靠的解決方案，且該方法也適用于動態(tài)尺寸優(yōu)化和形狀優(yōu)化，具有顯著的理論價值和工程意義。