基于rpim的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)

基于rpim的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)

1密度和接觸問題的數(shù)值計(jì)算結(jié)構(gòu)的連續(xù)優(yōu)化設(shè)計(jì)是結(jié)構(gòu)優(yōu)化領(lǐng)域之后的一個(gè)充滿挑戰(zhàn)的研究方向。連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化的優(yōu)點(diǎn)是能在不知道結(jié)構(gòu)拓?fù)湫螤畹那疤嵯?根據(jù)已知邊界條件和荷載條件確定比較合理的結(jié)構(gòu)形式,從而提出最佳形狀設(shè)計(jì)方案。連續(xù)體結(jié)構(gòu)的拓?fù)鋬?yōu)化本質(zhì)上是一種0-1離散變量的組合優(yōu)化問題。目前主要有三大類拓?fù)鋬?yōu)化方法:均勻化方法、變密度方法和漸進(jìn)結(jié)構(gòu)優(yōu)化方法。其中,變密度法是一種常用的拓?fù)鋬?yōu)化方法,其基本思想是不引入微結(jié)構(gòu),而是引入一種假想的相對密度在0～1之間可變的材料,人為假定相對密度和材料彈性模量之間的某種對應(yīng)關(guān)系,程序?qū)崿F(xiàn)簡單,計(jì)算效率高。變密度法中所指的密度是反映材料密度和材料特性之間對應(yīng)關(guān)系的一種偽密度。變密度法中常見的插值模型是,固體各向同性懲罰微結(jié)構(gòu)模型SIMP(SolidIsotropicMicrostructureswithPenalization)和材料特性的合理近似模型RAMP(RationalApproximationofMaterialProperties)。目前,連續(xù)體結(jié)構(gòu)的拓?fù)鋬?yōu)化問題基本上都是基于有限元數(shù)值計(jì)算方法。有限元法在處理網(wǎng)格畸變及網(wǎng)格移動(dòng)等問題時(shí),計(jì)算中需要不斷重新劃分和重構(gòu)網(wǎng)格以解決與原始網(wǎng)格線不一致的不連續(xù)和大變形問題。無網(wǎng)格方法(MeshlessMethod)是近年來迅速發(fā)展起來的一種新型的數(shù)值方法。這種方法采用基于點(diǎn)的近似,可以徹底或部分的消除網(wǎng)格,不需要網(wǎng)格的初始劃分和重構(gòu),不僅可以保證計(jì)算精度,而且可以減小計(jì)算的難度。目前已提出多種無網(wǎng)格方法主要有:光滑質(zhì)點(diǎn)流體動(dòng)力學(xué)方法(SPH)、無單元伽遼金法(EFG)、再生核粒子法(RKPM)以及無網(wǎng)格局部伽遼金法(MLPG)等。與有限元不同,無網(wǎng)格法中使用的近似函數(shù)大都是不具有插值特性,因此在基于Galerkin法的無網(wǎng)格中對于邊界條件的處理比較棘手。由Liu等提出的點(diǎn)插值方法(PIM)則較好的解決了這個(gè)問題。點(diǎn)插值方法的插值函數(shù)具有Delta函數(shù)性質(zhì),可以很方便的施加本質(zhì)邊界條件,不足之處是在計(jì)算插值函數(shù)時(shí)矩陣易于奇異。實(shí)際上,帶有多項(xiàng)式的徑向點(diǎn)插值法就可以有效地解決點(diǎn)插值法中出現(xiàn)的奇異性問題。本文利用無網(wǎng)格RPIM方法對二維線彈性連續(xù)體結(jié)構(gòu)的拓?fù)鋬?yōu)化問題進(jìn)行研究。在優(yōu)化過程中,選擇節(jié)點(diǎn)的相對密度作為設(shè)計(jì)變量,有效地抑制了點(diǎn)態(tài)棋盤格現(xiàn)象。最后通過算例分析證明了應(yīng)用無網(wǎng)格徑向點(diǎn)插值法進(jìn)行結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)的正確性和有效性。2彈性力學(xué)平面問題中非格徑點(diǎn)的插值方法2.1徑向基函數(shù)的引入在徑向點(diǎn)插值法中,計(jì)算域是用一系列點(diǎn)來離散的,每個(gè)點(diǎn)都有一定的影響域,某給定點(diǎn)處的位移是通過對該點(diǎn)的影響域中其他點(diǎn)處的位移進(jìn)行插值而得到。設(shè)二維域Ω中的任一函數(shù)u(x),可用徑向基和多項(xiàng)式基的線性組合表示為u(x)=n∑i=1Ri(x)ai+m∑j=1Ρj(x)bj=RΤ(x)a+ΡΤ(x)b(1)式中Ri(x)為徑向基函數(shù)(RBF),ai為Ri(x)的系數(shù),Pj(x)為坐標(biāo)xT=[x,y]中的單項(xiàng)式,bj為Pj(x)的系數(shù),n為RBF的項(xiàng)數(shù),m為多項(xiàng)式基函數(shù)的項(xiàng)數(shù)。為保證取得較好的穩(wěn)定性,通常取m<<n。在二維問題中,一般采用線性基:PT(x)=[1,x,y)]。徑向基的引入主要是為了避免單純采用多項(xiàng)式基時(shí)可能導(dǎo)致的剛度矩陣奇異性,目前常用的徑向基函數(shù)有以下四種形式:1)Multi-quadrics(MQ)Ri(x,y)=(r2i+(αcdc)2)q,αc≥02)Gaussian(EXP)Ri(x,y)=exp[-αc(ridc)2]3)ThinPlateSpline(TPS)Ri(x,y)=rηi4)LogarithmicRi(x,y)=rηilogri式中ri=√(x-xi)2+(y-yi)2?αc,dc,q和η為形狀參數(shù)。式(1)中的系數(shù)ai和bj可由u(x)通過影響域中的n個(gè)離散節(jié)點(diǎn)來確定,如第k個(gè)插值點(diǎn)為uk=u(xk,yk)=n∑i=1Ri(xk,yk)ai+m∑j=1Ρj(xk,yk)bjk=1,2,?,n(2)寫成矩陣形式為Us=R0a+Pmb(3)為了保證近似函數(shù)的唯一性,對系數(shù)要有附加條件,如n∑i=1Ρj(xi,yi)ai=ΡΤma=0j=1,2,?,m(4)聯(lián)立式(3)和式(4),得ˉUs={Us0}[R0ΡmΡΤm0][ab]=G{ab}(5)可將式(1)重新寫成:u(x)=RΤ(x)a+ΡΤ(x)b=[RΤ(x)ΡΤ(x)]{ab}(6)利用式(5),可得u(x)=[RΤ(x)ΡΤ(x)]G-1ˉUs=ˉΦΤ(x)ˉUs(7)式中ˉΦΤ(x)為RPIM的形函數(shù),可表示為ˉΦΤ(x)=[RΤ(x)ΡΤ(x)]G-1=[?1(x)?2(x)??n(x)?n+1(x)??n+m(x)](8)最終對應(yīng)于節(jié)點(diǎn)位移向量的RPIM形函數(shù)Φ(x)可表示為ΦΤ(x)=[?1(x)?2(x)??n(x)](9)式(7)可重寫成:u(x)=ΦΤ(x)Us=n∑i=1?iui(10)2.2單位外法線向量考慮彈性力學(xué)平面問題:{LΤσ+b=0在∈Ωσn=ˉt在∈Γtu=ˉu在∈Γu(11)式中L為應(yīng)變微分算子,ˉt為自然邊界上給定的表面力,ˉu為位移邊界上給定的位移,n為自然邊界上某點(diǎn)處的單位外法線向量。式(11)的變分形式可表示為∫Ω(Lδu)Τ(DLu)dΩ-∫ΩδuΤbdΩ-∫ΓtδuΤˉtdΓ=0(12)將式(10)代入上式,可得KU=F(13)式中節(jié)點(diǎn)的剛度矩陣KIJ和節(jié)點(diǎn)的荷載向量FI可分別表示為ΚΙJ=∫ΩBΤΙDBJdΩFΙ=∫ΩΦΤΙbdΩ+∫ΓtΦΤΙˉtdΓ其中BΙ=[?Ι,x00?Ι,y?Ι,y?Ι,x]D=E1-μ2[1μ0μ1000(1-μ)/2]3結(jié)構(gòu)拓?fù)鋬?yōu)化3.1密度場插值設(shè)計(jì)SIMP模型是工程中應(yīng)用最多的密度函數(shù)插值模型,通過引入懲罰因子使中間密度值向0、1兩端聚集,使連續(xù)變量的拓?fù)鋬?yōu)化模型能很好地逼近0-1離散變量的優(yōu)化模型,其數(shù)學(xué)模型為Eijkl(x)=ρp(x)E0ijkl(14)式中E0ijkl和Eijkl分別表示初始彈性模量和優(yōu)化后的彈性模量,p為懲罰因子。設(shè)計(jì)域內(nèi)任意一點(diǎn)的相對密度,可以用RPIM形函數(shù)對其影響域內(nèi)節(jié)點(diǎn)的相對密度插值得到,即ρg=np∑i=1Φiρi(15)式中ρi為第i個(gè)節(jié)點(diǎn)的相對密度,并確定為設(shè)計(jì)變量;Φi為RPIM形函數(shù),np為影響域內(nèi)的節(jié)點(diǎn)數(shù)目。這種插值方案保證了密度場的函數(shù)具有C0連續(xù)性,從而克服了點(diǎn)態(tài)棋盤格現(xiàn)象。本文以節(jié)點(diǎn)的相對密度作為優(yōu)化的設(shè)計(jì)變量,以柔度的最小化作為優(yōu)化的目標(biāo),以結(jié)構(gòu)整體的體積約束作為優(yōu)化的約束條件,基于SIMP模型建立線彈性結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)在靜力狀態(tài)下的數(shù)學(xué)模型,其形式如下:find.ρ(x),x∈Ωmin.c=FTUs.t.KU=FV=∫ΩρgdΩ=fV0,0<ρmin≤ρi≤1(16)式中K為整體剛度矩陣,U為位移列陣,F為荷載列陣,V為在設(shè)計(jì)變量狀態(tài)下的設(shè)計(jì)區(qū)域體積,V0為優(yōu)化前的設(shè)計(jì)區(qū)域體積,f為體積系數(shù)。為了避免計(jì)算中的奇異性,取密度下限值為ρmin=0.001。3.2oc法變量迭代格式求解拓?fù)鋬?yōu)化的算法有優(yōu)化準(zhǔn)則法(OC)、序列線性規(guī)劃法(SLP)和移動(dòng)漸近線法(MMA)等。利用優(yōu)化準(zhǔn)則法進(jìn)行優(yōu)化求解時(shí),最大的特點(diǎn)是對設(shè)計(jì)變量修改較大,收斂速度快,迭代次數(shù)少,且與結(jié)構(gòu)的大小和復(fù)雜程度無關(guān)。本文采用文獻(xiàn)中OC法的設(shè)計(jì)變量迭代格式,即ρnewi={max(ρmin,ρi-m)ifρiBηi≤max(ρmin,ρi-m)ρiBηiifmax(ρmin,ρi-m)<ρiBηi<min(1,ρi+m)min(1,ρi+m)ifmin(1,ρi+m)≤ρiBηi(17)式中Bi=-?c?ρi/λ?V?ρi其中λ為拉格朗日乘子,可以通過二分法來求解,m為移動(dòng)極限常數(shù),η為阻尼系數(shù)。移動(dòng)極限常數(shù)和阻尼系數(shù)的引入是為了保證優(yōu)化迭代的穩(wěn)定性。3.3各變量間的匹配本文參考有限元法的靈敏度分析方法,根據(jù)無網(wǎng)格法的離散原理和積分方法,利用伴隨靈敏度分析方法,求解目標(biāo)函數(shù)的靈敏度。通過添加零函數(shù)將目標(biāo)函數(shù)改寫為c=FΤU-?UΤ(ΚU-F)(18)式中?U為任意實(shí)向量。對式(18)關(guān)于節(jié)點(diǎn)密度求導(dǎo),得到?c?ρi=FΤ?U?ρi-?UΤ(?Κ?ρiU+Κ?U?ρi)=(FΤ-?UΤΚ)?U?ρi-?UΤ?Κ?ρiU(19)當(dāng)?U滿足伴隨方程FΤ-?UΤΚ=0時(shí),?U=U,式(19)轉(zhuǎn)化為?c?ρi=-UΤ?Κ?ρiU(20)這樣,目標(biāo)函數(shù)的靈敏度就轉(zhuǎn)化為求解剛度矩陣關(guān)于設(shè)計(jì)變量的靈敏度。由K=∫ΩρgpBTDBdΩ得?Κ?ρi=∫Ωpρgp-1ΦiBΤDBdΩ(21)體積約束對設(shè)計(jì)變量的靈敏度為?V?ρi=∫ΩΦidΩ(22)4計(jì)算在本小節(jié)中,將通過幾個(gè)經(jīng)典的拓?fù)鋬?yōu)化算例來說明本文所提出的優(yōu)化算法的可行性和有效性。4.1不同材料的pla和丙綸拓?fù)鋬?yōu)化設(shè)計(jì)考慮如圖1(a)所示的懸臂梁,長10m,寬10m,厚1m,左邊固定,右邊中點(diǎn)處作用一集中荷載F=1kN。將設(shè)計(jì)區(qū)域離散成441個(gè)節(jié)點(diǎn),根據(jù)節(jié)點(diǎn)的布置,將設(shè)計(jì)區(qū)域劃分為400個(gè)積分單元,每個(gè)積分單元設(shè)置2×2個(gè)高斯積分點(diǎn)。材料彈性模量E=3×108Pa,泊松比μ=0.3,取體積比為50%。利用本文的方法對圖1(a)所示模型進(jìn)行拓?fù)鋬?yōu)化設(shè)計(jì),迭代到47步時(shí)結(jié)果收斂,拓?fù)鋬?yōu)化結(jié)果如圖1(b)所示,拓?fù)鋬?yōu)化目標(biāo)函數(shù)值從0.4510減小到0.0478。圖1(c)為用有限元法進(jìn)行拓?fù)鋬?yōu)化的結(jié)果,迭代22次,拓?fù)鋬?yōu)化目標(biāo)函數(shù)值從0.4496減小到0.0489。圖1(d)為采用了靈敏度過濾技術(shù)的有限元拓?fù)鋬?yōu)化結(jié)果,迭代35次,拓?fù)鋬?yōu)化目標(biāo)函數(shù)值從0.4496減小到0.0535。從圖1可以看出,利用無網(wǎng)格徑向點(diǎn)插值法可以有效地對二維連續(xù)體結(jié)構(gòu)進(jìn)行拓?fù)鋬?yōu)化設(shè)計(jì),優(yōu)化結(jié)果更加趨于最優(yōu),但計(jì)算時(shí)間略有增加,同時(shí)選取節(jié)點(diǎn)的相對密度作為拓?fù)鋬?yōu)化設(shè)計(jì)變量可以有效地抑制有限元法中的棋盤格現(xiàn)象。4.2低影響單元拓?fù)鋬?yōu)化設(shè)計(jì)考慮如圖2(a)所示的懸臂梁,長10m,寬10m,厚1m,左邊固定,右邊下端點(diǎn)處作用一集中荷載F=1kN。仍將設(shè)計(jì)區(qū)域離散成441個(gè)節(jié)點(diǎn),將設(shè)計(jì)區(qū)域劃分為400個(gè)積分單元,每個(gè)積分單元設(shè)置2×2個(gè)高斯積分點(diǎn)。材料彈性模量E=3×108Pa,泊松比μ=0.3,取體積比為40%。利用本文的方法對圖2(a)所示模型進(jìn)行拓?fù)鋬?yōu)化設(shè)計(jì),迭代到89步時(shí)結(jié)果收斂,拓?fù)鋬?yōu)化結(jié)果如圖2(b)所示,拓?fù)鋬?yōu)化目標(biāo)函數(shù)值從1.9150減小到0.1036。圖2(c)為利用RPIM方法選取高斯點(diǎn)相對密度作為設(shè)計(jì)變量進(jìn)行拓?fù)鋬?yōu)化的結(jié)果,迭代59次,拓?fù)鋬?yōu)化目標(biāo)函數(shù)值從1.9150減小到0.0907。從圖2可以看出,用無網(wǎng)格徑向點(diǎn)插值法基于節(jié)點(diǎn)相對密度進(jìn)行結(jié)構(gòu)拓?fù)鋬?yōu)化可以有效地抑制基于高斯點(diǎn)密度進(jìn)行拓?fù)鋬?yōu)化中的點(diǎn)態(tài)棋盤格現(xiàn)象。5基于simp模型的結(jié)構(gòu)