基于整體-局部技術(shù)的結(jié)構(gòu)形狀優(yōu)化策略

吳志學(xué)

(揚(yáng)州大學(xué) 機(jī)械工程學(xué)院，揚(yáng)州 225127)

1 引 言

隨著現(xiàn)代工程設(shè)計(jì)方法和計(jì)算機(jī)輔助制造技術(shù)的快速發(fā)展，以提高結(jié)構(gòu)的承載能力、消除(或減少)材料的冗余和延長(zhǎng)使用壽命為目標(biāo)的結(jié)構(gòu)形狀優(yōu)化設(shè)計(jì)方法得到了迅速的發(fā)展和應(yīng)用[1,2]。各種優(yōu)化方法在減小應(yīng)力集中[3-6]、結(jié)構(gòu)輕量化設(shè)計(jì)[7-9]、提高疲勞斷裂強(qiáng)度[10-12]以及減小接觸應(yīng)力[13,14]等方面得到了廣泛應(yīng)用。但由于一般結(jié)構(gòu)力學(xué)響應(yīng)具有很強(qiáng)的非線性，再加上大量的迭代計(jì)算過(guò)程，使得結(jié)構(gòu)形狀優(yōu)化方法的應(yīng)用推廣一直受到很大限制。即使計(jì)算能力大幅度提高，為提高計(jì)算精度而精確離散計(jì)算模型也是不現(xiàn)實(shí)的，特別是針對(duì)大尺度結(jié)構(gòu)和復(fù)合材料結(jié)構(gòu)。

整體-局部技術(shù)(又稱(chēng)為子模型技術(shù))是從整體模型的局部區(qū)域中獲得更加精確解的有限單元技術(shù)。整體-局部技術(shù)基本思想是選取整體模型的局部區(qū)域作為局部模型,并重新劃分單元網(wǎng)格，然后用切割邊界上的整體模型分析結(jié)果作為局部模型的邊界條件并重新計(jì)算，從而得到更為精確的局部分析結(jié)果。該方法廣泛應(yīng)用于各種復(fù)雜結(jié)構(gòu)的分析計(jì)算，如應(yīng)用于復(fù)雜零件的疲勞損傷分析[15]、功能梯度材料的計(jì)算分析[16]、壓力容器表面裂紋的應(yīng)力強(qiáng)度因子計(jì)算[17]以及冷卻塔非線性結(jié)構(gòu)分析[18]等。顯然，將整體-局部技術(shù)與結(jié)構(gòu)形狀優(yōu)化設(shè)計(jì)結(jié)合起來(lái)，將大量的迭代計(jì)算過(guò)程局限于局部模型，這為有效解決求解效率與精度之間的矛盾提供了可能。但目前關(guān)于這方面的研究報(bào)道較少。理論上只要切割邊界上的邊界條件足夠準(zhǔn)確，局部模型的計(jì)算結(jié)果精度就有保證。然而，形狀優(yōu)化過(guò)程中，局部形狀的變化不可避免地會(huì)導(dǎo)致切割邊界的邊界條件的變化，進(jìn)而影響到局部模型的計(jì)算結(jié)果精度。因此，如何處理切割邊界的邊界條件變化是將整體-局部技術(shù)應(yīng)用于結(jié)構(gòu)形狀優(yōu)化的關(guān)鍵。文獻(xiàn)[11]在對(duì)某飛機(jī)油孔形狀優(yōu)化時(shí)采用了整體-局部技術(shù)，將由整體模型分析得到的載荷與位移直接移植到切割邊界作為固定的邊界條件處理，但沒(méi)有對(duì)優(yōu)化結(jié)果精度進(jìn)行討論分析。

本文將以無(wú)限大板中心孔形狀優(yōu)化為例，分析并討論有無(wú)切割邊界的邊界條件來(lái)修正兩種優(yōu)化方案對(duì)計(jì)算效率和精度的影響。在此基礎(chǔ)上，探索一種將整體-局部技術(shù)與優(yōu)化法相結(jié)合的結(jié)構(gòu)形狀優(yōu)化策略,以解決大尺度結(jié)構(gòu)局部形狀優(yōu)化設(shè)計(jì)中的求解效率與精度之間的矛盾。

2 理論基礎(chǔ)與問(wèn)題描述

圖1為在遠(yuǎn)場(chǎng)載荷作用下的帶缺口板，其中L0為初始缺口形狀，稱(chēng)之為設(shè)計(jì)邊界。優(yōu)化目標(biāo)是在給定的設(shè)計(jì)域Ω范圍內(nèi),通過(guò)優(yōu)化設(shè)計(jì)邊界形狀獲得最小應(yīng)力集中系數(shù)。形狀優(yōu)化是典型的力學(xué)反問(wèn)題。顯然，設(shè)計(jì)邊界形狀與沿切割邊界Lc的應(yīng)力或應(yīng)變分布之間有著確定的對(duì)應(yīng)關(guān)系。也就是說(shuō)，優(yōu)化過(guò)程中設(shè)計(jì)邊界形狀的變化會(huì)導(dǎo)致切割邊界Lc的應(yīng)力或應(yīng)變分布改變。因此，如果將整體-局部技術(shù)應(yīng)用到結(jié)構(gòu)形狀優(yōu)化中，必須解決切割邊界的邊界條件變化的影響問(wèn)題。

圖1 整體-局部模型的形狀優(yōu)化

解決這一問(wèn)題有兩種可能的方案。一種是不考慮切割邊界的邊界條件變化，直接應(yīng)用整體-局部技術(shù)進(jìn)行形狀優(yōu)化。如果設(shè)計(jì)邊界形狀的變化引起的切割邊界的邊界條件變化很小，進(jìn)而對(duì)優(yōu)化結(jié)果影響足夠小，則這種方案是可行的。本文以無(wú)限大板中心孔形狀優(yōu)化為例，分析由局部形狀的變化所導(dǎo)致切割邊界的邊界條件變化情況。圖2 為在遠(yuǎn)場(chǎng)拉伸載荷作用下的帶中心孔平板模型。假設(shè)局部模型的切割邊界為半徑等于Ri的圓(Ri>R,R為圓孔半徑)，在極坐標(biāo)系下，點(diǎn)(Ri,θ)沿圓周方向的正應(yīng)力分量(σθ)r =Ri為[19]

圖2 受遠(yuǎn)場(chǎng)拉伸載荷作用的帶中心孔平板

作為一階近似，由孔形狀改變而導(dǎo)致的沿切割邊界的應(yīng)力分布變化量可以表述為

(2)

(3)

由式(1～3)可以看出，類(lèi)似于孔對(duì)應(yīng)力分布的影響，孔形狀改變導(dǎo)致的沿切割邊界的應(yīng)力分布變化具有明顯的局部性特征。這意味著只要Ri足夠大，孔形狀改變導(dǎo)致的沿切割邊界的應(yīng)力分布變化就可以忽略。上述情況下不考慮邊界條件修正而直接應(yīng)用整體-局部技術(shù)進(jìn)行形狀優(yōu)化是可行的。顯然，該方案的代價(jià)是局部模型的切割邊界要足夠大，計(jì)算效率受其影響明顯。

另一種解決方案是在應(yīng)用整體-局部技術(shù)進(jìn)行形狀優(yōu)化時(shí)考慮邊界條件修正。整體模型的單元數(shù)量比局部模型多，特別是對(duì)于大尺度結(jié)構(gòu)，因此，計(jì)算效率主要取決于對(duì)整體模型分析的次數(shù)。由于邊界條件值是通過(guò)對(duì)整體模型分析得到的，所以如何減少局部模型的邊界條件修正次數(shù)是提高優(yōu)化效率的關(guān)鍵。

3 基于整體-局部技術(shù)的形狀優(yōu)化策略

3.1 結(jié)構(gòu)形狀優(yōu)化方法

根據(jù)對(duì)設(shè)計(jì)邊界形狀的描述，形狀優(yōu)化可以分為參數(shù)化和非參數(shù)化兩種方法。參數(shù)化方法用長(zhǎng)度、半徑或樣條曲線的控制點(diǎn)坐標(biāo)等參數(shù)來(lái)描述設(shè)計(jì)邊界形狀。該方法的優(yōu)點(diǎn)是設(shè)計(jì)變量相對(duì)較少，優(yōu)化算法也比較成熟；其缺點(diǎn)是解空間受限于所選擇的參數(shù)方程，一般情況下只能獲得局部最優(yōu)解。非參數(shù)化方法是選擇有限元模型設(shè)計(jì)邊界上的節(jié)點(diǎn)坐標(biāo)作為設(shè)計(jì)變量。與參數(shù)化方法相比，非參數(shù)化方法的解空間僅受限于單元離散化程度，有望獲得全局最優(yōu)解；其缺點(diǎn)也非常明顯，由于設(shè)計(jì)變量過(guò)多，如何提高優(yōu)化效率是關(guān)鍵。近年來(lái)提出的一種以模擬生物生長(zhǎng)為基礎(chǔ)的無(wú)梯度優(yōu)化算法為提高優(yōu)化效率提供了有效途徑[3-5,8-11]。該方法的收斂速度與設(shè)計(jì)變量數(shù)目無(wú)關(guān)，因此非常適于求解復(fù)雜的三維優(yōu)化問(wèn)題。綜合考慮解的精度與計(jì)算效率，本文采用以?xún)?yōu)化準(zhǔn)則為基礎(chǔ)的非參數(shù)化無(wú)梯度優(yōu)化法求解形狀優(yōu)化問(wèn)題。

對(duì)于基于應(yīng)力的形狀優(yōu)化問(wèn)題，根據(jù)優(yōu)化準(zhǔn)則，應(yīng)力最小化問(wèn)題轉(zhuǎn)變?yōu)榍蠼庋卦O(shè)計(jì)邊界的均勻分布的應(yīng)力問(wèn)題，表達(dá)式為

Min Max|σ-σref|

Subject toLi∈Ω

andXi(Pi)∈Li

(4)

式中σ與σref分別為沿設(shè)計(jì)邊界Li的應(yīng)力分量和參考應(yīng)力；Xi(Pi)為設(shè)計(jì)變量，本文選擇有限元模型設(shè)計(jì)邊界上的節(jié)點(diǎn)坐標(biāo)作為設(shè)計(jì)變量。在每次迭代中，節(jié)點(diǎn)i的位移量(沿設(shè)計(jì)邊界法向)按表達(dá)式(5)計(jì)算得到

(5)

式中σi為節(jié)點(diǎn)i上的應(yīng)力；C為步長(zhǎng)系數(shù)，其初始值可根據(jù)經(jīng)驗(yàn)或試算選取，在隨后的迭代計(jì)算過(guò)程中，通過(guò)檢測(cè)應(yīng)力差值來(lái)實(shí)時(shí)調(diào)整C以保證各點(diǎn)的di相等；σref為參考應(yīng)力。本文將設(shè)計(jì)邊界固定點(diǎn)上的應(yīng)力作為參考應(yīng)力[20]。

根據(jù)式(5)調(diào)整所有節(jié)點(diǎn)的位置，相應(yīng)地改變了設(shè)計(jì)邊界的形狀。對(duì)具有新邊界形狀的模型重新進(jìn)行有限元分析，并且重復(fù)上述過(guò)程，直至設(shè)計(jì)邊界上的應(yīng)力分布達(dá)到均勻狀態(tài)。采用如下性能指標(biāo)ER來(lái)控制優(yōu)化邊界上的應(yīng)力均勻程度,

(6)

式中σmax和σmin分別為設(shè)計(jì)邊界上的最大和最小應(yīng)力值。實(shí)際應(yīng)用中通過(guò)給定ER的許用值τ來(lái)控制優(yōu)化解的精度，并以此作為優(yōu)化結(jié)束的判據(jù)。如當(dāng)τ=0.01時(shí)，表示按算術(shù)平均應(yīng)力估算得到的均勻應(yīng)力與真實(shí)最優(yōu)值的相對(duì)誤差小于1%。

3.2 基于整體-局部技術(shù)的結(jié)構(gòu)形狀優(yōu)化方法

設(shè)Ld為最終優(yōu)化邊界如圖1所示，相應(yīng)的沿切割Lc的應(yīng)力分布函數(shù)為σ(u)*。顯然，當(dāng)次優(yōu)模型的設(shè)計(jì)邊界Li趨近于Ld時(shí)，優(yōu)化模型中沿Lc的應(yīng)力分布函數(shù)σ(u)i將收斂于σ(u)*。從整體-局部技術(shù)角度來(lái)看，常規(guī)形狀優(yōu)化相當(dāng)于每一次優(yōu)化迭代過(guò)程都對(duì)局部模型切割邊界的邊界條件修正一次。為了減少對(duì)整體模型分析的次數(shù)，本文提出了一個(gè)雙循環(huán)優(yōu)化程序： ① 局部?jī)?yōu)化循環(huán),采用局部模型進(jìn)行優(yōu)化獲得設(shè)計(jì)邊界Li的形狀； ② 整體優(yōu)化循環(huán),在整體模型中將設(shè)計(jì)邊界修正為L(zhǎng)i，并獲得次優(yōu)模型的優(yōu)化結(jié)果以及新的局部模型切割邊界條件。重復(fù)上述過(guò)程，直至設(shè)計(jì)邊界上的應(yīng)力分布達(dá)到均勻狀態(tài)?；谡w-局部技術(shù)的結(jié)構(gòu)形狀優(yōu)化流程如圖3所示。

圖3 基于整體-局部技術(shù)的結(jié)構(gòu)形狀優(yōu)化流程

4 實(shí)例分析與討論

在本節(jié)的計(jì)算中，所有的材料常數(shù)均為楊氏模量E=200 GPa，泊松比v=0.3。設(shè)定τ=0.01。

4.1 無(wú)切割邊界條件修正的優(yōu)化

以遠(yuǎn)場(chǎng)受拉伸載荷作用下無(wú)限大平板中心孔優(yōu)化問(wèn)題為例討論無(wú)切割邊界條件修正優(yōu)化方法的優(yōu)缺點(diǎn)。為模擬帶孔無(wú)限大平板，將孔直徑D(D=2R)設(shè)為有限板寬W的1/100，并取L=1.5W，如圖2所示。有限元模型采用8節(jié)點(diǎn)四邊形等參單元，基于線彈性和平面應(yīng)力假設(shè)。板左右邊界±L/2處施加均勻拉伸應(yīng)力σ0=1 MPa。在這種情況下，有限元計(jì)算得到圓孔附近凈截面應(yīng)力集中系數(shù)為2.97，與這一問(wèn)題的理論解3.0非常接近。

首先，采用整體模型對(duì)中心孔形狀進(jìn)行優(yōu)化，經(jīng)過(guò)27次迭代得到應(yīng)力集中系數(shù)為2.112，并且沿優(yōu)化曲線獲得了滿(mǎn)足給定精度要求的均勻應(yīng)力分布，如圖5(b)所示。為了研究孔形狀變化對(duì)沿不同切割邊界應(yīng)力分布的影響，選擇不同半徑Ri的圓弧作為切割邊界，如圖2所示。圖4(a)給出了初始圓弧孔模型(σ(Ri,θ))和優(yōu)化后孔模型(σ*(Ri,θ))在不同切割邊界上的應(yīng)力分布(1/4對(duì)稱(chēng)模型)?？梢钥闯觯诓煌那懈钸吔缟?，σ(Ri,θ)與σ*(Ri,θ)之間的差值是不同的。本文采用沿切割邊界的最大應(yīng)力差來(lái)近似表示由于孔形狀改變引起的應(yīng)力分布變化。圖4(b)給出了無(wú)量綱化的最大應(yīng)力差max{|σ(Ri,θ)-σ*(Ri,θ)|/σ0}隨無(wú)量綱化半徑Ri/R的變化情況。可以看出，孔形狀的變化對(duì)沿切割邊界的應(yīng)力分布變化的影響具有局部化特征。當(dāng)切割邊界的尺寸大于孔的三倍直徑時(shí)，由孔形狀變化引起的應(yīng)力最大誤差不超過(guò)10%。這種應(yīng)力變化的局部化特性使得直接應(yīng)用整體-局部?jī)?yōu)化設(shè)計(jì)方法成為可能。

圖4 初始圓弧孔模型和優(yōu)化后孔模型在不同切割邊界上的應(yīng)力分布

其次，在無(wú)切割邊界位移修正的條件下，采用基于整體-局部策略的優(yōu)化方法得到優(yōu)化解。表1給出了應(yīng)力集中系數(shù)的優(yōu)化結(jié)果。由表1可知，即使當(dāng)切割邊界非常接近設(shè)計(jì)邊界時(shí)，優(yōu)化也可以降低應(yīng)力集中。隨著切割邊界的增大，應(yīng)力集中逐漸逼近整體模型的最優(yōu)解。需要強(qiáng)調(diào)的是，由于邊界條件的不精確，單純從局部模型來(lái)計(jì)算應(yīng)力集中系數(shù)是不正確的。

表1 無(wú)切割邊界位移修正的優(yōu)化結(jié)果

圖5分別給出了孔邊優(yōu)化形狀和相應(yīng)的應(yīng)力分布(1/4對(duì)稱(chēng))曲線。由于結(jié)果之間過(guò)于接近導(dǎo)致難以區(qū)別，圖中沒(méi)有給出Ri/R=6-10的結(jié)果。可以看出，隨著Ri/R的增加，設(shè)計(jì)邊界的形狀向整體模型的最優(yōu)形狀收斂，應(yīng)力分布更加均勻。需要注意的是，當(dāng)Ri/R大于3時(shí)，應(yīng)力集中的誤差(見(jiàn)表1)以及孔形狀優(yōu)化結(jié)果誤差都很小。這意味著當(dāng)切割邊界的尺寸足夠大時(shí)，如果直接將整體-局部策略應(yīng)用到優(yōu)化方法中，優(yōu)化結(jié)果可以具有相當(dāng)高的精度。

圖5 不同切割邊界條件下的優(yōu)化結(jié)果(1/4對(duì)稱(chēng))

4.2 有切割邊界條件修正的優(yōu)化

分別選擇Ri/R=1.5和2的圓弧作為切割邊界，采用如圖3所示的優(yōu)化程序進(jìn)行求解。對(duì)于Ri/R=1.5和2的切割邊界，優(yōu)化后得到的應(yīng)力集中系數(shù)分別為2.119和2.117。應(yīng)力集中系數(shù)隨整體優(yōu)化循環(huán)次數(shù)的變化如圖6所示。可以看出，雙循環(huán)優(yōu)化方法具有收斂速度快的優(yōu)點(diǎn)。而且，當(dāng)設(shè)計(jì)邊界的形狀逐漸趨近于最優(yōu)形狀時(shí)，局部?jī)?yōu)化循環(huán)的迭代次數(shù)也迅速減少。值得一提的是，由于收斂速度快，僅在第一次修改邊界條件后就能得到較好的結(jié)果。這意味著，只需要增加有限計(jì)算量，應(yīng)用雙循環(huán)優(yōu)化程序就能獲得局部形狀優(yōu)化的滿(mǎn)意結(jié)果。

圖6 應(yīng)力集中系數(shù)隨整體優(yōu)化循環(huán)次數(shù)的變化(括號(hào)中數(shù)值分別為Ri/R =1.5和2時(shí)局部?jī)?yōu)化循環(huán)的迭代次數(shù))

4.3 開(kāi)式壓力機(jī)加強(qiáng)筋板曲線的優(yōu)化設(shè)計(jì)

圖7為100 t開(kāi)式壓力機(jī)結(jié)構(gòu)。壓力機(jī)加強(qiáng)筋板是關(guān)鍵區(qū)域，該位置局部應(yīng)力過(guò)大，使用中容易產(chǎn)生疲勞裂紋。設(shè)計(jì)目標(biāo)是找到加強(qiáng)筋板曲線的最優(yōu)形狀使其具有最小的應(yīng)力集中。

圖7 開(kāi)式壓力機(jī)結(jié)構(gòu)圖和局部模型

首先，采用整體模型對(duì)壓力機(jī)進(jìn)行分析。有限元模型采用20節(jié)點(diǎn)六面體單元。由于x-y平面對(duì)稱(chēng)，只分析結(jié)構(gòu)的1/2。有限元網(wǎng)格包括實(shí)體單元112560個(gè)，節(jié)點(diǎn)223959個(gè)，加載邊界條件如 圖8(a)所示。在2.93 GHz CPU和2GB RAM內(nèi)存的計(jì)算機(jī)上，整體模型的計(jì)算時(shí)間為778 s。除了對(duì)加強(qiáng)筋板區(qū)域的最大應(yīng)力進(jìn)行初步計(jì)算外，該分析結(jié)果主要提供了局部模型沿切割邊界的邊界條件。

選擇a=190,b=330,c=190 mm的六面體作為局部模型的切割邊界，如圖7所示。局部模型有9947個(gè)單元和46705個(gè)節(jié)點(diǎn)，如圖8(b)所示。在同一臺(tái)計(jì)算機(jī)上，局部模型的計(jì)算時(shí)間為28 s，遠(yuǎn)小于整體模型的計(jì)算時(shí)間。在這種情況下，計(jì)算得到加強(qiáng)筋板表面的最大von Mises應(yīng)力為205.8 MPa。

圖8 有限元網(wǎng)格和邊界條件

一般情況下，應(yīng)力沿加強(qiáng)筋板厚度方向分布是不均勻的。這里，選取加強(qiáng)筋板設(shè)計(jì)表面中間線段作為設(shè)計(jì)線，如圖8(b)所示。設(shè)計(jì)結(jié)束后，沿加強(qiáng)筋板厚度方向掃描設(shè)計(jì)線就可得到約束在設(shè)計(jì)域內(nèi)的設(shè)計(jì)邊界(曲面)。將設(shè)計(jì)線上節(jié)點(diǎn)坐標(biāo)定義為設(shè)計(jì)變量。采用如圖3所示的雙循環(huán)優(yōu)化程序進(jìn)行求解，整個(gè)優(yōu)化過(guò)程的計(jì)算時(shí)間為6312 s(約為1.75 h)。

圖9給出了最大von Mises應(yīng)力隨優(yōu)化循環(huán)次數(shù)的變化情況。結(jié)果表明，雙循環(huán)優(yōu)化策略實(shí)現(xiàn)了快速收斂，并且隨著設(shè)計(jì)邊界形狀接近最優(yōu)輪廓，局部?jī)?yōu)化循環(huán)的迭代次數(shù)迅速減少。計(jì)算總量包括5次整體模型計(jì)算和69次局部模型計(jì)算。因此，基于整體-局部策略的優(yōu)化方法可以在合理的計(jì)算代價(jià)下獲得滿(mǎn)意的結(jié)果。

圖9 最大von Mises應(yīng)力隨優(yōu)化循環(huán)次數(shù)的變化(括號(hào)中的數(shù)值為局部?jī)?yōu)化循環(huán)的迭代次數(shù))

加強(qiáng)筋板的最終優(yōu)化形狀及其沿設(shè)計(jì)線的von Mises應(yīng)力分布如圖10所示(圖10(b)橫坐標(biāo)為無(wú)量綱化處理的x方向坐標(biāo)值x/x0，x0為設(shè)計(jì)線x方向總長(zhǎng))?？梢钥闯觯瑑?yōu)化后的形狀與初始形狀相比，沿設(shè)計(jì)邊界的應(yīng)力分布極為均勻。該優(yōu)化形狀的最大von Mises應(yīng)力為143.9 MPa，比初始最大應(yīng)力降低了30%。

圖10 加強(qiáng)筋板的優(yōu)化結(jié)果

5 結(jié) 論

本文研究了采用整體-局部策略進(jìn)行結(jié)構(gòu)形狀優(yōu)化的可行性，目的在于解決大尺度結(jié)構(gòu)局部形狀優(yōu)化設(shè)計(jì)中的效率與精度之間的矛盾。對(duì)有無(wú)切割邊界位移修正條件下的優(yōu)化策略進(jìn)行了分析和討論。提出了一種結(jié)合整體-局部策略的非參數(shù)無(wú)梯度形狀優(yōu)化方法來(lái)優(yōu)化缺口形狀的雙循環(huán)優(yōu)化方案。設(shè)計(jì)實(shí)例表明，雙循環(huán)優(yōu)化程序在結(jié)構(gòu)局部形狀優(yōu)化中可以很好地兼顧求解效率與精度。另外，由于常見(jiàn)的工程結(jié)構(gòu)都不可避免地存在缺口，本文提出的形狀優(yōu)化方法可以很容易地推廣到類(lèi)似大型結(jié)構(gòu)的形狀設(shè)計(jì)中。