考慮承載面最大位移約束的結(jié)構(gòu)拓?fù)鋬?yōu)化方法

龍凱，陳卓，谷春璐，王選

1. 華北電力大學(xué) 新能源電力系統(tǒng)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 102206 2. 合肥工業(yè)大學(xué) 土木與水利工程學(xué)院，合肥 230009

連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)指在給定的區(qū)域內(nèi)，尋求結(jié)構(gòu)的某種布局(如結(jié)構(gòu)有無(wú)孔洞、孔洞位置、數(shù)量和結(jié)構(gòu)連接方式等)，使其能夠在滿(mǎn)足指定約束條件下的設(shè)計(jì)目標(biāo)最優(yōu)[1-2]。拓?fù)鋬?yōu)化在早期研究中以全局性性能如結(jié)構(gòu)柔順度、振動(dòng)頻率為設(shè)計(jì)指標(biāo)，發(fā)展到考慮局部應(yīng)力、疲勞損傷等性能，從解決工程實(shí)際問(wèn)題的角度出發(fā)，以改善結(jié)構(gòu)綜合性能[3]。

位移或變形量是結(jié)構(gòu)產(chǎn)品剛度的重要衡量指標(biāo)。Rong等[4-5]在拓?fù)鋬?yōu)化設(shè)計(jì)中引入多點(diǎn)位移約束，實(shí)現(xiàn)了對(duì)結(jié)構(gòu)變形的控制。Huang和Xie[6]提出采用柔順度目標(biāo)和節(jié)點(diǎn)位移約束以實(shí)現(xiàn)對(duì)全局和局部剛度的控制。Zuo和Xie[7]提出一個(gè)全局化位移指標(biāo)，采用雙向進(jìn)化式拓?fù)鋬?yōu)化方法實(shí)現(xiàn)了結(jié)構(gòu)剛度設(shè)計(jì)。Qiao和Liu[8]提出了幾何平均位移指標(biāo)用于衡量拓?fù)鋬?yōu)化結(jié)構(gòu)的剛度。Li和Kim[9]從工程實(shí)際出發(fā)，提出了多相材料結(jié)構(gòu)位移約束下的輕量化設(shè)計(jì)列式和求解算法。Long等[10]提出了多點(diǎn)位移約束拓?fù)鋬?yōu)化問(wèn)題的二次規(guī)劃列式與求解算法。針對(duì)局部翹曲變形，朱繼宏等[11]提出了多點(diǎn)保形的概念，用來(lái)衡量多點(diǎn)間相對(duì)位移變化程度和變形協(xié)調(diào)。基于保形指標(biāo)的拓?fù)鋬?yōu)化設(shè)計(jì)方法在考慮變形方向性、諧響應(yīng)振動(dòng)和幾何非線(xiàn)性問(wèn)題中得到進(jìn)一步拓展[12-14]。

在工程問(wèn)題中，有一類(lèi)工程結(jié)構(gòu)承受分布力載荷作用，例如受風(fēng)載荷的建筑、承受氣動(dòng)力的飛行翼面、舵體等[3]。此類(lèi)結(jié)構(gòu)承載表面的節(jié)點(diǎn)變形量需要滿(mǎn)足剛度設(shè)計(jì)要求。由于承載面節(jié)點(diǎn)數(shù)目龐大，并且在優(yōu)化迭代過(guò)程中，最大位移發(fā)生的位置會(huì)發(fā)生變化，若采用多點(diǎn)位移約束的處理方式，不僅敏度求解中的伴隨工況數(shù)量龐大，有限元分析計(jì)算量大，而且會(huì)造成約束方程數(shù)目增多而難以?xún)?yōu)化求解。

針對(duì)上述問(wèn)題，提出考慮承載面最大位移約束的拓?fù)鋬?yōu)化模型。遵循獨(dú)立連續(xù)映射(Independent Continuous Mapping, ICM)建模方法[15-16]，引入倒變量中間變量，將原有拓?fù)鋬?yōu)化問(wèn)題轉(zhuǎn)換為一系列嚴(yán)格正定的二次規(guī)劃問(wèn)題，并采用高效、穩(wěn)健的序列二次規(guī)劃算法求解。將提出模型獲得的優(yōu)化結(jié)果與傳統(tǒng)的體積比約束下柔順度最小化列式結(jié)果、多點(diǎn)位移約束下體積最小化結(jié)果做對(duì)比，說(shuō)明提出列式與求解方法的可行性和優(yōu)越性。

1 體積比約束下的拓?fù)鋬?yōu)化列式

在現(xiàn)有的拓?fù)鋬?yōu)化方法中，以體積比為約束、柔順度最小化為目標(biāo)的拓?fù)鋬?yōu)化列式最為典型，數(shù)學(xué)表達(dá)式為[1]

find:ρ

minimize:c=FTu

(1)

式中:ρ為密度矢量;分量ρe為單元e的相對(duì)密度值，用于表征單元e的有無(wú)；c為系統(tǒng)柔順度值；K、u和F分別為總剛度陣、位移列陣和載荷列陣；V和V0分別為優(yōu)化后體積和初始結(jié)構(gòu)體積；f為設(shè)定體積比；密度下限ρmin用于避免結(jié)構(gòu)分析中的數(shù)值奇異性，這里取值ρmin=10-3；NE為設(shè)計(jì)域內(nèi)單元總數(shù)。

SIMP插值模型通過(guò)建立彈性模量(或剛度陣)與密度變量的假設(shè)關(guān)系，驅(qū)使在優(yōu)化過(guò)程中密度變量向0或1靠近，其表達(dá)式為[1]

(2)

柔順度c對(duì)密度變量ρe的一階導(dǎo)數(shù)為[1]

(3)

式中:ue為單元e的位移列陣。拓?fù)鋬?yōu)化列式(1)在得到結(jié)構(gòu)響應(yīng)量敏度的前提下，將敏度數(shù)值代入到連續(xù)型優(yōu)化算法，即可實(shí)現(xiàn)優(yōu)化求解。在拓?fù)鋬?yōu)化研究領(lǐng)域中，常見(jiàn)的優(yōu)化求解算法有準(zhǔn)則法(Optimality Criteria, OC)[17]、凸線(xiàn)性化(Convex Linearization, CONLIN)[18]、移動(dòng)漸進(jìn)線(xiàn)法(Method of Moving Asymptotes, MMA)[19-20]等。

2 考慮承載面最大位移約束的拓?fù)鋬?yōu)化列式

2.1 拓?fù)鋬?yōu)化列式

盡管列式(1)在拓?fù)鋬?yōu)化研究領(lǐng)域得到了廣泛的應(yīng)用，但存在著以下缺點(diǎn)[1]：① 柔順度值缺乏工程意義；② 體積比是人為給定的數(shù)值，但實(shí)際工程應(yīng)用中，拓?fù)鋬?yōu)化結(jié)果僅用于結(jié)構(gòu)傳力路徑的確定，工程可用的結(jié)構(gòu)需要進(jìn)一步通過(guò)尺寸優(yōu)化、形狀優(yōu)化等后續(xù)工作完成。

對(duì)于承受分布載荷的工程結(jié)構(gòu)，承載面的變形量需要進(jìn)行嚴(yán)格控制。假設(shè)承載面上每一個(gè)節(jié)點(diǎn)位移滿(mǎn)足約束限值，以結(jié)構(gòu)輕量化設(shè)計(jì)為目標(biāo)的拓?fù)鋬?yōu)化列式為

find:ρ

minimize:V/V0

(4)

為了解決上述困難，這里提出的解決方案是采用承載面上最大節(jié)點(diǎn)位移約束取代多點(diǎn)位移約束方程，即

find:ρ

minimize:V/V0

(5)

式(5)中的最大值函數(shù)具有不可導(dǎo)性，從而無(wú)法利用現(xiàn)有的連續(xù)型優(yōu)化算法求解。受最大應(yīng)力約束拓?fù)鋬?yōu)化問(wèn)題的啟發(fā)[21-22]，此問(wèn)題可采用包絡(luò)函數(shù)包括KS函數(shù)、p范數(shù)(p-norm)和p平均(p-mean)函數(shù)對(duì)最大值函數(shù)進(jìn)行連續(xù)光滑處理，這里選擇采用KS函數(shù)，最大值約束的代理模型為

(6)

式中:μ為KS包絡(luò)函數(shù)參數(shù)，并滿(mǎn)足:

(7)

根據(jù)數(shù)值經(jīng)驗(yàn)，當(dāng)μ取值較大時(shí)，容易引起優(yōu)化求解中的數(shù)值振蕩；當(dāng)μ取有限值時(shí)，包絡(luò)函數(shù)值與真實(shí)最大值之間具有一定的差異性。為了消除這種差異性，引入修正系數(shù)以調(diào)整兩者間的關(guān)系，即[23-24]

(8)

式中:cKS為調(diào)整系數(shù)，數(shù)學(xué)表達(dá)式為

(9)

2.2 敏度分析

對(duì)式(6)求偏導(dǎo)可得

(10)

由鏈?zhǔn)角髮?dǎo)法則得

(11)

自由度j對(duì)應(yīng)的位移表達(dá)為位移向量的函數(shù)，即

(12)

將式(12)代入到式(11)可得

(13)

對(duì)平衡方程Ku=F兩邊求導(dǎo)可得

(14)

假設(shè)載荷F與密度變量無(wú)關(guān)，則

(15)

將式(15)代入到式(13)中可得

(16)

建立伴隨方程

(17)

式中：λ為伴隨向量，則敏度表達(dá)式為

(18)

由式(18)可得，當(dāng)采用提出的最大位移指標(biāo)時(shí)，伴隨方程對(duì)應(yīng)一次有限元分析，相比較多點(diǎn)位移約束問(wèn)題，可以大大減少伴隨方程求解計(jì)算量。

2.3 獨(dú)立連續(xù)映射法顯式化過(guò)程

與常見(jiàn)拓?fù)鋬?yōu)化方法不同，ICM以各類(lèi)結(jié)構(gòu)響應(yīng)量為約束，重量最小化為目標(biāo)。通常情況下，柔順度、節(jié)點(diǎn)位移和應(yīng)力等各類(lèi)結(jié)構(gòu)響應(yīng)量難以求其二階敏度，在獲取一階敏度信息的情況下，可獲取其一階近似表達(dá)。在引入倒變量情形下，目標(biāo)體積函數(shù)可以解析獲取對(duì)設(shè)計(jì)變量的二階敏度值。綜合目標(biāo)函數(shù)和約束函數(shù)的近似表達(dá)，可構(gòu)建出二次規(guī)劃問(wèn)題的模型，配合穩(wěn)健高效的序列二次規(guī)劃(Sequential Quadratic Programming, SQP)算法求解。本節(jié)將遵循ICM建模方法，實(shí)現(xiàn)目標(biāo)函數(shù)和約束函數(shù)的顯式化表達(dá)。

定義設(shè)計(jì)變量為密度變量ρe的倒數(shù)函數(shù)，即

(19)

由此推導(dǎo)體積函數(shù)對(duì)設(shè)計(jì)變量的一和二階敏度表達(dá)式為

(20a)

(20b)

式中:ve為單元e的體積。

由式(20b)組成的海塞矩陣(Hessian)的表達(dá)式為

H=

(21)

與優(yōu)化列式(4)對(duì)比可知，若直接采用ρe作為設(shè)計(jì)變量，體積函數(shù)對(duì)密度變量ρe的一階敏度為常數(shù)，二階敏度為0。當(dāng)采用倒變量形式的設(shè)計(jì)變量xe時(shí)，由式(21)易得，海塞矩陣非對(duì)角線(xiàn)元素為0且對(duì)角線(xiàn)元素恒定大于0，即海塞矩陣嚴(yán)格正定且具有可分離性質(zhì)。

利用鏈?zhǔn)角髮?dǎo)法則可得

(22)

式中:λe為單元伴隨向量。

基于上述敏度信息，位移約束函數(shù)采用一階泰勒展開(kāi)得到顯式表達(dá)式：

(23)

式中:上標(biāo)l代表第l輪優(yōu)化迭代。

體積目標(biāo)函數(shù)采用二階泰勒展開(kāi)得到顯式表達(dá)式，并忽略其中的常數(shù)項(xiàng)可得

find:x

(24)

式中:矩陣B與體積函數(shù)的一階敏度相關(guān)。優(yōu)化列式(24)是標(biāo)準(zhǔn)的二次規(guī)劃模型，為了克服拓?fù)鋬?yōu)化中設(shè)計(jì)變量龐大的困難，通常轉(zhuǎn)化為對(duì)偶模型并采用SQP求解[15-16, 25]，優(yōu)化迭代至與體積比相對(duì)變化率滿(mǎn)足收斂條件，即

|V(l)-V(l-1)|/V(l)≤ε

(25)

收斂率ε通常為1‰。

借鑒數(shù)字圖像處理中的過(guò)濾技術(shù)消除拓?fù)鋬?yōu)化中常見(jiàn)的棋盤(pán)格現(xiàn)象，具體過(guò)程可參考文獻(xiàn)[15-16，24]。值得注意的是，由于過(guò)濾平均的影響，將導(dǎo)致結(jié)構(gòu)邊界存在著一定數(shù)量的灰度單元。

3 數(shù)值算例與討論

本節(jié)采用數(shù)值算例來(lái)驗(yàn)證提出模型的可行性和有效性。在各算例中，材料的彈性模量和泊松比值分別為1 MPa和0.3，平面結(jié)構(gòu)厚度為1 mm，單元尺寸為1 mm。初始結(jié)構(gòu)的單元密度值取1，即設(shè)計(jì)區(qū)域全部為實(shí)體材料。

算例1如圖1所示的長(zhǎng)懸臂梁結(jié)構(gòu)，結(jié)構(gòu)尺寸為240 mm×84 mm，頂部4層單元為非設(shè)計(jì)區(qū)域，即單元密度值恒定為1。左端全部固定，頂部受到總載荷值為100 kN的力作用，載荷均勻分布在各個(gè)節(jié)點(diǎn)上。假設(shè)頂層y方向最大位移許用值為8 mm。KS函數(shù)參數(shù)μ=20。采用提出的拓?fù)鋬?yōu)化列式和求解算法，所得到的拓?fù)鋬?yōu)化結(jié)果如圖2所示。

圖1 長(zhǎng)懸臂梁結(jié)構(gòu)示意圖Fig.1 Illustration of long cantilever structure

由圖2可知，拓?fù)鋬?yōu)化結(jié)果頂部的變形盡可能保持水平且上下起伏，這表明可能多處節(jié)點(diǎn)的位移達(dá)到約束限值。最終結(jié)構(gòu)體積比為0.309，最大節(jié)點(diǎn)位移值為7.999 mm。

為了驗(yàn)證提出方法的有效性，采用拓?fù)鋬?yōu)化列式(5)和MMA算法求解，拓?fù)鋬?yōu)化結(jié)果如圖3所示。

由圖2和圖3對(duì)比可知，2種不同方法下的拓?fù)鋬?yōu)化結(jié)果及變形特點(diǎn)類(lèi)似。MMA算法求解下的最優(yōu)結(jié)構(gòu)體積比為0.310，最大節(jié)點(diǎn)位移值為8.001 mm。上述結(jié)果也表明，本文提出的優(yōu)化列式具有可行性。

圖2 所提方法的拓?fù)鋬?yōu)化結(jié)果Fig.2 Optimized topology obtained from the proposed method

圖3 采用拓?fù)鋬?yōu)化列式(5)和 MMA算法的拓?fù)浣Y(jié)果Fig.3 Optimized topology obtained from Eq.(5) and MMA algorithm

為了對(duì)比與傳統(tǒng)拓?fù)鋬?yōu)化結(jié)果的異同，以體積比0.309為約束條件，采用剛度最大化列式(1)，得到的拓?fù)鋬?yōu)化結(jié)果如圖4所示，優(yōu)化結(jié)構(gòu)頂部受載區(qū)域最大位移為10.167 mm，比設(shè)定的許用位移值8 mm大27.1%。頂面結(jié)構(gòu)變形保持為直線(xiàn)，朝著自由端傾斜，這與圖2和圖3中最優(yōu)結(jié)構(gòu)的變形特點(diǎn)明顯不同。上述結(jié)果表明，傳統(tǒng)的剛度最大化列式(1)不具有控制局部最大變形的能力。

圖4 拓?fù)鋬?yōu)化列式(1)下的拓?fù)鋬?yōu)化結(jié)果Fig.4 Optimized topology using Eq.(1)

算例2如圖5所示的平面結(jié)構(gòu)，結(jié)構(gòu)尺寸為180 mm×120 mm，底部4層單元為非設(shè)計(jì)區(qū)域，即單元密度值恒定為1,下端兩角點(diǎn)全約束。下端面受到合力F=100 kN作用，載荷均勻分配到下端面每個(gè)節(jié)點(diǎn)上，假設(shè)底層y方向最大位移許用值為8 mm，最優(yōu)拓?fù)渑c結(jié)構(gòu)變形如圖6所示。

圖5 平面結(jié)構(gòu)示意圖Fig.5 Illustration of plane structure

圖6 所提列式下的拓?fù)鋬?yōu)化結(jié)果與變形Fig.6 Optimized topology obtained from the proposed equation

為了對(duì)比與多點(diǎn)位移約束下拓?fù)鋬?yōu)化結(jié)果的異同，設(shè)置底部距離左端點(diǎn)1/6，2/6，3/6，4/6和5/6長(zhǎng)度位置點(diǎn)位移約束條件，形成多點(diǎn)位移約束拓?fù)鋬?yōu)化列式，同樣采用ICM方法求解[15-16]，最優(yōu)拓?fù)渑c結(jié)構(gòu)變形如圖7所示。位移約束設(shè)置位置對(duì)應(yīng)的節(jié)點(diǎn)位移均達(dá)到許用值8 mm，但底部最大節(jié)點(diǎn)位移為9.433 mm，這說(shuō)明即使采用多點(diǎn)位移約束，也無(wú)法嚴(yán)格控制承載面內(nèi)每一點(diǎn)的變形。值得注意的是，多點(diǎn)位移約束的敏度分析包含5個(gè)伴隨工況，即對(duì)應(yīng)5次結(jié)構(gòu)有限元計(jì)算，而本文提出的優(yōu)化列式將多點(diǎn)位移約束凝聚為單約束函數(shù)，僅包含1個(gè)伴隨工況。上述結(jié)果證明了提出方法的有效性和優(yōu)越性。

圖7 多點(diǎn)位移約束下的拓?fù)鋬?yōu)化結(jié)果Fig.7 Optimized topology obtained from multiple nodal displacement constraints

算例3三維結(jié)構(gòu)如圖8所示，結(jié)構(gòu)尺寸為90 mm×40 mm×28 mm，底部四角點(diǎn)全約束，頂部受到總載荷F=26 390 kN的力作用，載荷均勻分配到上端面每個(gè)節(jié)點(diǎn)上。頂部2層單元為非設(shè)計(jì)區(qū)域，即單元密度值恒定為1。設(shè)頂部區(qū)域y向最大許用位移為18 mm。采用提出的方法得到的拓?fù)鋬?yōu)化結(jié)果如圖9(a)所示，最優(yōu)結(jié)構(gòu)體積比為0.125，承載面最大節(jié)點(diǎn)位移為17.996 mm；設(shè)定0.125為體積比約束,采用柔順度最小優(yōu)化列式(1),拓?fù)鋬?yōu)化結(jié)果如圖9(b)所示，承載面最大節(jié)點(diǎn)位移為19.056 mm。

圖8 三維結(jié)構(gòu)示意圖Fig.8 Illustration of 3D structure

由圖9可知，不同的拓?fù)鋬?yōu)化列式下對(duì)應(yīng)的拓?fù)鋬?yōu)化結(jié)果有所不同。類(lèi)似于二維結(jié)構(gòu)，本文提出方法對(duì)對(duì)該三維結(jié)構(gòu)具有控制承載面最大位移的效果。即在相同的材料用量下，采用柔順度最小化目標(biāo)，優(yōu)化結(jié)構(gòu)的節(jié)點(diǎn)最大位移值超過(guò)許用位移值約5.89%。算例結(jié)果證明了提出列式和優(yōu)化求解算法在三維結(jié)構(gòu)優(yōu)化問(wèn)題中的適用性。

圖9 不同列式下的拓?fù)鋬?yōu)化結(jié)果Fig.9 Optimized topology obtained from different equations

4 結(jié) 論

針對(duì)工程中受分布載荷的連續(xù)體結(jié)構(gòu)，提出承載面最大位移約束下的拓?fù)鋬?yōu)化設(shè)計(jì)模型。該模型通過(guò)KS包絡(luò)函數(shù)實(shí)現(xiàn)了多點(diǎn)位移最大值的凝聚，推導(dǎo)了相應(yīng)的敏度表達(dá)式。基于ICM方法，形成了具有正定可分離海塞矩陣的二次規(guī)劃列式，采用序列二次規(guī)劃法優(yōu)化求解，得到以下結(jié)論：

1) 與傳統(tǒng)的體積比約束下柔順度最小化列式對(duì)比，提出方法具有控制局部區(qū)域最大位移的效果。優(yōu)化結(jié)果證明了提出方法的有效性。

2) 與多點(diǎn)位移約束對(duì)比，提出方法更精確、有效控制承載面的最大變形量。且將多點(diǎn)位移約束凝聚為單一約束，不僅減少了伴隨工況計(jì)算工作量，且有利于優(yōu)化求解。

3) 提出的優(yōu)化列式和相應(yīng)求解算法有望在幾何非線(xiàn)性位移約束問(wèn)題中進(jìn)一步拓展。