基于應(yīng)力和基頻的漸進(jìn)優(yōu)化法

胡瑞姣

（重慶大學(xué)機(jī)械工程學(xué)院，重慶 400044）

漸進(jìn)結(jié)構(gòu)優(yōu)化方法（evolutionary structural optimization，ESO）［1］是1993年由 Xie和 Steven提出的，其基本思想是通過逐漸將無效或低效的材料刪除，剩下的結(jié)構(gòu)將趨于優(yōu)化。在優(yōu)化迭代中，方法采用固定的有限元網(wǎng)格，對存在的材料單元設(shè)其材料數(shù)編號為1，而對不存在的材料單元設(shè)其材料數(shù)編號為0。方法設(shè)計變量少，每個單元只需一個標(biāo)明單元存在與否的狀態(tài)變量。當(dāng)計算結(jié)構(gòu)剛度矩陣等特性時，不計材料數(shù)編號為0的單元特性。因此隨著大量“死”單元的去掉，結(jié)構(gòu)特性參數(shù)分析和重分析的方程求解數(shù)大大減少。由于概念上的簡潔性和應(yīng)用上的有效性，ESO具有較大的理論價值和很好的工程應(yīng)用前景，目前已應(yīng)用于應(yīng)力、位移、剛度、頻率、穩(wěn)定性和臨界應(yīng)力約束問題的優(yōu)化。關(guān)于ESO方法的研究可參見文獻(xiàn)［2-8］，由此可見ESO的理論不斷得到完善。

ESO進(jìn)行拓?fù)鋬?yōu)化的關(guān)鍵在于用一個合適的準(zhǔn)則來評價每個單元對結(jié)構(gòu)特定行為的貢大小，以刪除一些貢獻(xiàn)較小的單元。而以前的準(zhǔn)則研究大都是基于靜力問題的，然而結(jié)構(gòu)的動態(tài)特性又是影響產(chǎn)品性能的關(guān)鍵因素之一。為了使ESO獲得更具有工程應(yīng)用價值的優(yōu)化結(jié)果，在優(yōu)化過程中需同時考慮應(yīng)力和頻率方面的約束。又由于低階模態(tài)特性基本決定了產(chǎn)品的動態(tài)性能，因此利用ESO進(jìn)行基于應(yīng)力和頻率約束的拓?fù)鋬?yōu)化，即從滿足應(yīng)力約束和頻率約束的結(jié)構(gòu)中系統(tǒng)地刪除低效材料，得到重量最小化的拓?fù)鋬?yōu)化。

1 ESO應(yīng)力優(yōu)化

ESO首先是針對應(yīng)力優(yōu)化而提出的，應(yīng)力優(yōu)化的準(zhǔn)則為:逐漸去掉結(jié)構(gòu)中的低應(yīng)力材料，是剩下的結(jié)構(gòu)更有效地承擔(dān)載荷，從而應(yīng)力分布更加均勻。所以將單元的VonMises應(yīng)力σvme與整個結(jié)構(gòu)的最大VonMises應(yīng)力進(jìn)行比較，便可確定每個單元的應(yīng)力水平。

在每一次有限元分析完成后，從模型中刪除滿足下列條件的所有單元式中RRi是當(dāng)前刪除率（RR），通過使用相同的RRi值，重復(fù)執(zhí)行有限元分析與單元刪除，直至達(dá)到穩(wěn)定狀態(tài)，也即是在當(dāng)前步已不能進(jìn)一步刪除材料。這時在材料刪除率中引入進(jìn)化率（ER）

通過增加刪除率，再一次進(jìn)行有限元分析與單元刪除，直至達(dá)到一個新的穩(wěn)定狀態(tài)。如此的一個進(jìn)化過程繼續(xù)至獲得期望的最佳結(jié)構(gòu)。

2 ESO頻率優(yōu)化

為了確定結(jié)構(gòu)修改的最好位置，需要進(jìn)行靈敏度分析來確定須刪除的單元。對有限元動平衡方程微分，可得到第i單元的特征值靈敏度。

在有限元分析中，結(jié)構(gòu)動特性的一般特征值問題為

其中:［K］為總剛度矩陣;［M］為質(zhì)量矩陣ωi為第i階固有頻率;{φi}為相應(yīng)于ωi的特征矢量。固有頻率ωi以及相對應(yīng)的特征矢量{φi}通過瑞利商相互關(guān)聯(lián)

其中模態(tài)剛度ki和模態(tài)質(zhì)量mi定義為

假定在當(dāng)前結(jié)構(gòu)中刪除單元j，對（2）式微分得到由于刪除該單元引起的第i階頻率變化量

假設(shè)單元增減前后，特征矢量{φi}不變，其中}是j單元的振型，因此

其中［ke］和［me］分別是刪除單元j的單元剛度矩陣和單元質(zhì)量矩陣，將式（8）和（9）帶入式（7）

在模態(tài)分析中，常將特征向量{φei}相對于質(zhì)量矩陣［M］正則化，使得模態(tài)質(zhì)量mi等于1，這種情況下靈敏度計算公式可簡化為

由此可發(fā)現(xiàn)，只要在有限元分析之后得到特征參數(shù)，敏度值就可以在單元層次上求得。

頻率靈敏度是在假設(shè)刪除單元前后振型不發(fā)生改變的基礎(chǔ)上推到得到的，具有一定的數(shù)值不穩(wěn)定性，以至于棋盤格問題嚴(yán)重，為了抑制棋盤格，采用靈敏度再分配技術(shù)。即將單元靈敏度分配到節(jié)點，節(jié)點I的靈敏度等于相關(guān)單元靈敏度的平均，即

再將節(jié)點靈敏度返回單元，也就是說單元的再分配靈敏度為本單元節(jié)點靈敏度的平均值。

頻率優(yōu)化有最大化或最小化指定的某階或某幾節(jié)頻率、最大化指定兩頻率間隔等問題。要增大某階固有頻率，優(yōu)化迭代時只需刪除靈敏度較大的單元即可，相反，欲減小某階固有頻率刪除靈敏度值為負(fù)的較小單元。

3 ESO應(yīng)力和頻率優(yōu)化

3.1 考慮應(yīng)力和頻率的結(jié)構(gòu)拓?fù)鋬?yōu)化模型

在結(jié)構(gòu)拓?fù)鋬?yōu)化中，基于應(yīng)力和基頻約束的優(yōu)化模型建立如下:

式中bi為設(shè)計變量，取0和1（0表示刪除單元，1表示保留單元）。W為結(jié)構(gòu)總重量，wi是第i個單元的重量，σmax是最大應(yīng)力，σ*、f*分別是應(yīng)力、頻率約束值。

3.2 優(yōu)化步驟

算法的具體實現(xiàn)步驟:

（1）定義結(jié)構(gòu)能占有的最大允許物理區(qū)域，并作為初始設(shè)計結(jié)構(gòu);用有限元精細(xì)離散網(wǎng)格，給定基于應(yīng)力的初始刪除率RR和進(jìn)化率ER，基于頻率的刪除率RE。

（2）指定載荷，約束，單元特性等。

（3）靜力分析，按照式（1）選取單元集合N，集合N中單元個數(shù)為n。

（4）模態(tài)分析，計算集合N中每個單元的再分配靈敏度αi，其中αi＞0的單元個數(shù)為m。

（5）基于頻率刪除單元的個數(shù)為f，f=d*RE（其中d為當(dāng)前迭代步活的單元），比較n，m以及f的大小，取3個當(dāng)中最小值作為本次迭代要刪除的單元個數(shù)g。

（6）從集合N中刪除g個頻率靈敏度大于0的單元。

（7）孤點以及單鉸支檢驗，并刪除。

（8）如果沒有滿足式（1）的單元，結(jié)構(gòu)拓?fù)溥_(dá)到一個穩(wěn)態(tài)，按式（2）增加基于應(yīng)力的刪除率。

（9）重復(fù)3－8，直到滿足某一準(zhǔn)則限。

如果要減小基頻，只需將（4）改為:模態(tài)分析，按照式（11）計算集合N中每個單元的αi，其中αi＜0的單元個數(shù)為m。與基于應(yīng)力的ESO相比，第4、5、6步是對基于應(yīng)力的ESO的一個修正;與基于頻率的ESO相比，第3、5步是對基于應(yīng)力的ESO的一個修正。

4 數(shù)值算例

4.1 對角固支的矩形板的拓?fù)鋬?yōu)化設(shè)計

矩形板處于平面應(yīng)力狀態(tài)，一組對角點固定，底邊中點豎向載荷F=100 N，厚度t=0.01 m，彈性模量E=70 GPa，泊松比 ν=0.3，ρ=2 700 kg/m3。初始設(shè)計區(qū)域為0.15×0.1 m的矩形區(qū)域，模型劃分為45×30個四邊形有限元網(wǎng)格單元。

設(shè)定初始基于應(yīng)力的單元刪除率RR=0.01、單元進(jìn)化率ER=0.01，基于頻率的刪除率RE=0.01（迭代過程中保持不變）。進(jìn)化結(jié)構(gòu)優(yōu)化進(jìn)程結(jié)果見圖2所示:

圖1 矩形板拓?fù)鋬?yōu)化計算模型

圖2 矩形板進(jìn)化過程

圖2是相應(yīng)6個刪除率RRi下的穩(wěn)態(tài)結(jié)構(gòu)，經(jīng)過15個穩(wěn)態(tài)迭代，當(dāng)RR=0.15時最大應(yīng)力由最初的28 185 Pa增加到29 485 Pa，頻率由2 441 Hz增加到3 239 Hz，增加了33%。最優(yōu)拓?fù)浣Y(jié)構(gòu)與基于頻率的ESO［6］拓?fù)浣Y(jié)構(gòu)有微小的差異，但鑒于算法的是以應(yīng)力和頻率為約束的ESO，一定程度上說明了算法的合理性。但方法只允許刪除單元，隨著剩余單元的減少得到結(jié)構(gòu)平滑性越來越不好，并且邊界粗超不均勻，有明顯的棋盤格。

4.2 兩端固支的矩形板拓?fù)鋬?yōu)化設(shè)計

矩形板結(jié)構(gòu)初始設(shè)計的尺寸為2H×H（H=50 mm），彈性模量E=100 GPa，泊松比v=0.3，板厚t=1 mm，ρ=7 000 kg/m3，豎向荷載F=1 N。有限元分析選用四邊形單元將幾何模型分成50×25個正方形固定網(wǎng)格（圖3）。

設(shè)定初始基于應(yīng)力的單元刪除率RR=0.01、單元進(jìn)化率ER=0.01，基于頻率的刪除率RE=0.01（迭代過程中保持不變）。進(jìn)化結(jié)構(gòu)優(yōu)化進(jìn)程結(jié)果如圖4所示:

圖3 矩形板拓?fù)鋬?yōu)化計算模型

圖4 矩形板拓?fù)溥M(jìn)程

圖4給出了6個不同迭代步的拓?fù)浣Y(jié)構(gòu)。在優(yōu)化前期與基于應(yīng)力的ESO拓?fù)浣Y(jié)構(gòu)沒有區(qū)別，隨著單元的減小，在中后期與基于應(yīng)力的ESO拓?fù)浣Y(jié)構(gòu)呈現(xiàn)一定的差異。算例中，最大應(yīng)力由最初的327 497 Pa增加到459 547 Pa，頻率由4 710降低到1 559，減少了67%。算法的最大應(yīng)力比基于應(yīng)力的ESO［8］要大，不如基于應(yīng)力的ESO最優(yōu)拓?fù)浣Y(jié)構(gòu)應(yīng)力均勻，更說明了算法是靜態(tài)性能優(yōu)化和動態(tài)性能優(yōu)化的一種平衡，作為一種啟發(fā)式算法，可根據(jù)工程需要，選擇合適的拓?fù)浣Y(jié)構(gòu)。

5 結(jié)論

綜合考慮了影響結(jié)構(gòu)性能的關(guān)鍵因素，提出一種以應(yīng)力和第一階固有頻率為約束的ESO，算例說明了拓?fù)鋬?yōu)化時算法一定程度上平衡了靜態(tài)性能和動態(tài)性能，使靜態(tài)性能和動態(tài)性能同時得到優(yōu)化。但由于算法中被刪除的單元在后續(xù)迭代中不能恢復(fù)的缺點，以及頻率靈敏度基于振型不變的基礎(chǔ)上推到計算的，算法還需繼續(xù)研究改進(jìn)。

［1］XIE Y M，STEVEN O P.A simple evolutionary procedure for structural optimization ［J］.Comput.Struct.1993，（49）:885-896

［2］TANSKANEN P.The evolutionary structural optimization method:Theoretical aspects［J］.Comput Methods Appl Mech Engrg，2002，24（6）:548-552

［3］XIE Y M，STEVEN G P.A simple approach to structural frequency optimization［J］.Comput Struct，1994，53（6）:148-152

［4］CHU D N，XIE Y M，HIRA A，et al.Onvarious aspects of evolutionary st-ructural optimization for problemswith stiffness const raint［J］.Finite Elements Anal，1997，24（2）:197-201

［5］林家浩.結(jié)構(gòu)動力優(yōu)化設(shè)計發(fā)展綜述［J］.力學(xué)進(jìn)展，1983，13（4）:423-431

［6］楊振興，榮見華，傅建林.三維結(jié)構(gòu)的頻率拓?fù)鋬?yōu)化設(shè)計［J］.振動與沖擊，2006，25（3）:44-47

［7］榮見華，鄭健龍，徐飛鴻.結(jié)構(gòu)動力修改及優(yōu)化設(shè)計［M］.北京:人民教育出版社，2001

［8］孫圣權(quán).基于ESO的液壓機(jī)下橫梁結(jié)構(gòu)拓?fù)鋬?yōu)化及雙向進(jìn)化結(jié)構(gòu)優(yōu)化方法研究［D］.重慶:重慶大學(xué)，2008