基于階躍函數(shù)的連續(xù)體拓?fù)鋬?yōu)化問(wèn)題

基于階躍函數(shù)的連續(xù)體拓?fù)鋬?yōu)化問(wèn)題

1基于邊界類條件的拓?fù)鋬?yōu)化作為結(jié)構(gòu)優(yōu)化的高水平領(lǐng)域，在國(guó)內(nèi)外科學(xué)家的共同努力下，做出了迅速的發(fā)展。目前在連續(xù)體拓?fù)鋬?yōu)化中占據(jù)主導(dǎo)地位的是均勻化法和SIMP(SolidIsotropicMaterialwithPenalization)法。因?yàn)樗鼈兌际侵饕芯繂卧挠泻蜔o(wú),所以把它們統(tǒng)稱為單元類方法。這類方法從本質(zhì)上講是一個(gè)0-1整數(shù)規(guī)劃問(wèn)題。但是由于這一問(wèn)題是不能保證存在最優(yōu)解的病態(tài)問(wèn)題,所以不能直接求解。為了求解,通常采用松弛設(shè)計(jì)變量和增加約束兩類辦法。均勻化方法是把拓?fù)鋬?yōu)化的基結(jié)構(gòu)分成許多帶有孔洞的胞元,以每一個(gè)胞元孔洞的幾何參數(shù)為設(shè)計(jì)變量,從而把拓?fù)鋬?yōu)化問(wèn)題轉(zhuǎn)化為尺寸優(yōu)化問(wèn)題來(lái)求解。SIMP法則是把每一個(gè)單元內(nèi)的材料屬性作為常數(shù),單元密度作為設(shè)計(jì)變量,通過(guò)對(duì)單元密度懲罰的方式來(lái)得到趨近0-1的拓?fù)浣Y(jié)構(gòu)。此外,Beckers等人通過(guò)引進(jìn)周長(zhǎng)約束直接求解整數(shù)0-1規(guī)劃,隋允康等人把拓?fù)渥兞繌某叽鐑?yōu)化變量中提取出來(lái)的獨(dú)立映射反演(ICM)方法也都屬于單元類方法。這一類方法在拓?fù)鋬?yōu)化的研究過(guò)程中起到了重大的作用,現(xiàn)在已經(jīng)被廣泛用于求解多種約束,多種物理場(chǎng)耦合等很多拓?fù)鋬?yōu)化新興領(lǐng)域之中。但是由于它們都是以結(jié)構(gòu)的單元為研究對(duì)象,所以不容易處理拓?fù)錁?gòu)型演化過(guò)程中的邊界問(wèn)題。而以邊界移動(dòng)來(lái)確定結(jié)構(gòu)構(gòu)型的邊界類算法在邊界的處理上則具有非常明顯的優(yōu)勢(shì)。但是標(biāo)準(zhǔn)的邊界類算法只是通過(guò)邊界形狀改變的靈敏度來(lái)確定構(gòu)型邊界的移動(dòng)速度,不能起到改變結(jié)構(gòu)拓?fù)涞淖饔?所以大多數(shù)情況下它只能作為一種形狀優(yōu)化算法來(lái)應(yīng)用。為了能夠更好地應(yīng)用邊界類算法的優(yōu)點(diǎn)來(lái)處理拓?fù)鋬?yōu)化中的邊界問(wèn)題,J.A.Sethian和AndreasWiegmann把用于描述材料邊界運(yùn)動(dòng)的水平集概念引入到結(jié)構(gòu)優(yōu)化之中。水平集算法的主要過(guò)程是由優(yōu)化的目標(biāo)和約束來(lái)確定邊界的移動(dòng)速度,然后由Hamilton-Jacobi方程求解水平集函數(shù),從而確定結(jié)構(gòu)的拓?fù)錁?gòu)型。由于水平集函數(shù)可以使邊界在移動(dòng)的同時(shí)產(chǎn)生拓?fù)涓淖?所以水平集方法在拓?fù)鋬?yōu)化中也取得了很好的數(shù)值結(jié)果。但是因?yàn)樗郊椒ú荒苁惯h(yuǎn)離邊界的結(jié)構(gòu)內(nèi)部的拓?fù)浒l(fā)生改變,所以當(dāng)初始構(gòu)型取的不好時(shí)不易收斂到最優(yōu)解。盡管水平集方法在拓?fù)鋬?yōu)化的應(yīng)用中還存在一些問(wèn)題,但是它可以使拓?fù)鋬?yōu)化構(gòu)型的邊界處理變得簡(jiǎn)單。最近,Belytschko和郭旭分別提出了基于水平集函數(shù)的求解方法。他們把有限元離散后的水平集函數(shù)節(jié)點(diǎn)值作為拓?fù)湓O(shè)計(jì)變量,通過(guò)階躍函數(shù)把水平集函數(shù)轉(zhuǎn)化為0-1形式建立優(yōu)化模型,然后通過(guò)有限元離散和優(yōu)化算法進(jìn)行求解,從而確定結(jié)構(gòu)的拓?fù)錁?gòu)型。其實(shí),這種算法從本質(zhì)上講是求解整數(shù)0-1規(guī)劃的一種連續(xù)化方法。數(shù)值算例證明了這種方法的有效性。但是在這兩篇文獻(xiàn)中也可以發(fā)現(xiàn)一對(duì)矛盾。在這兩篇文獻(xiàn)中都引進(jìn)了三角函數(shù)來(lái)光滑階躍函數(shù)。如果其中的光滑區(qū)域取的比較小,問(wèn)題的求解過(guò)程就接近于水平集方法。但是通過(guò)對(duì)拓?fù)鋵?dǎo)數(shù)的概念進(jìn)行分析后,可以看出這種接近于水平集的求解方法很難改變離邊界較遠(yuǎn)的結(jié)構(gòu)內(nèi)部的拓?fù)洹Ｈ绻饣瑓^(qū)域取的比較大,拓?fù)鋵?dǎo)數(shù)就可以在離邊界比較遠(yuǎn)的結(jié)構(gòu)內(nèi)部發(fā)揮作用,結(jié)構(gòu)的拓?fù)涓淖兙妥兊酶尤菀?優(yōu)化初始階段數(shù)值收斂速度也就越快,但是這樣做卻會(huì)影響收斂到最終拓?fù)錁?gòu)型的速度。為了解決這個(gè)矛盾,在文的基礎(chǔ)上,采用sigmoid函數(shù)來(lái)光滑階躍函數(shù),并通過(guò)對(duì)sigmoid函數(shù)中陡度參數(shù)的調(diào)整來(lái)更加方便地控制拓?fù)鋬?yōu)化的進(jìn)程。這樣既可以在拓?fù)鋬?yōu)化的初始階段采用較小的陡度參數(shù),使結(jié)構(gòu)內(nèi)部的拓?fù)涓淖兤饋?lái)更加容易,加快收斂速度;又可以在拓?fù)鋬?yōu)化的后期采用比較大陡度參數(shù)來(lái)確保以較快的速度收斂到最優(yōu)拓?fù)錁?gòu)型。2等效積分?jǐn)M合拓?fù)鋬?yōu)化模型在連續(xù)體拓?fù)鋬?yōu)化模型中,需要引進(jìn)一個(gè)0-1變量的特征函數(shù)χ(x)來(lái)確定基結(jié)構(gòu)內(nèi)的點(diǎn)是否屬于最優(yōu)拓?fù)錁?gòu)型χ(x)={0x∈D\Ω1x∈Ω(1)式中:D代表基結(jié)構(gòu);Ω代表拓?fù)錁?gòu)型所占區(qū)域。對(duì)于二維線彈性體,不考慮體積力,以整體柔順性為目標(biāo),材料體積為約束的連續(xù)體拓?fù)鋬?yōu)化模型P原為Μinimizeu,?L(u,χ)=∫?ΩitudS=∫DEijklεij(u)εkl(u)χ(x)dΩsubjectto:a(u,v,χ)=L(v,χ),u|?Dn=u0,?v∈UV(χ)=∫Dχ(x)dΩ≤Vmax(2)其中a(u,v,χ)=∫DEijklεij(u)εkl(v)χ(x)dΩ,L(v,χ)=∫?ΩitvdS,L(u,χ)=∫DEijklεij(u)εkl(u)χ(x)dΩ為由最小位能原理推出的柔順性在整個(gè)基結(jié)構(gòu)上的表達(dá),a(u,v,χ)=L(v,χ)為平衡方程的等效積分弱形式,V(χ)≤Vmax為結(jié)構(gòu)的體積約束,?Ωt和?Ωu分別為結(jié)構(gòu)的力和位移邊界。在本文中,采用水平集函數(shù)來(lái)建立這個(gè)拓?fù)鋬?yōu)化模型,水平集函數(shù)的形式為?(x)＞0?x∈Ω\?Ω?(x)=0?x∈?Ω?(x)＜0?x∈D\Ω(3)即?(x)=0代表拓?fù)錁?gòu)型的邊界,?(x)>0代表拓?fù)錁?gòu)型的內(nèi)部,?(x)<0代表基結(jié)構(gòu)內(nèi)不屬于拓?fù)錁?gòu)型的部分。為了把水平集函數(shù)轉(zhuǎn)化為拓?fù)鋬?yōu)化模型中的0-1形式,需要引進(jìn)階躍函數(shù)Η(?(x))={0,?＜01,?≥0(4)此外,在后面求靈敏度時(shí)還需要用到階躍函數(shù)的導(dǎo)數(shù)(δ函數(shù))Η′(?(x))=δ(?(x)){=0,?≠0≠0,?=0(5)綜上,用水平集函數(shù)表述的連續(xù)體拓?fù)鋬?yōu)化模型PTDF為MinimizeL(u,?)=∫?ΩtudS=∫DEijklεij(u)εkl(u)H(?)dΩsubjectto:a(u,v,?)=L(v,?),u|?Ω=u0,?v∈UV(?)=∫DH(?)dΩ≤Vmax(6)其中a(u,v,?)=∫DEijklεij(u)εkl(v)H(?)dΩ,L(v,?)=∫?ΩtvdS。3基結(jié)構(gòu)有限單元的剛度變化本文采用雙線性四節(jié)點(diǎn)單元建立有限元模型。位移場(chǎng)函數(shù)u的插值形式為uh(x)=4∑l=1Νl(x)ul(7)這里:Nl(x)為雙線性四節(jié)點(diǎn)單元的插值形函數(shù);ul為節(jié)點(diǎn)I上的位移值。對(duì)于水平集函數(shù),采用相同的插值形函數(shù)進(jìn)行離散?h(x)=4∑l=1Νl(x)?l(8)式中?l為水平集函數(shù)在節(jié)點(diǎn)I上的值。由PTDF中平衡方程的等效積分弱形式∫DEijklεij(u)εkl(v)H(?)dΩ-∫?ΩitvdS=0和式(7)、式(8)兩個(gè)離散形式,可以推出有限元方程ΚuΝ=fΝ(9)式中:uN為節(jié)點(diǎn)位移;fN為節(jié)點(diǎn)力;Κ=ΝE∑e-1ke?ΝE為基結(jié)構(gòu)中劃分的單元總數(shù),ke為單元?jiǎng)偠染仃?它的形式為ke=∫ΩeBΤDBΗ(?)dΩ(10)這里B為應(yīng)變-位移矩陣,D為材料彈性常數(shù)張量的矩陣形式,Ωe為單元所占區(qū)域。但是用式(10)直接做數(shù)值積分求得的單元?jiǎng)偠染仃囋谟?jì)算中不能使結(jié)構(gòu)收斂到0-1形式。在這個(gè)問(wèn)題的處理上,文獻(xiàn)和文獻(xiàn)都采用了類似SIMP法中罰的方式。文獻(xiàn)中對(duì)彈性能和體積約束中的階躍函數(shù)選取了不同的光滑函數(shù),其中的乘方關(guān)系體現(xiàn)的就是罰的意思。文獻(xiàn)則是把每個(gè)單元水平集函數(shù)的節(jié)點(diǎn)平均值看作單元的相對(duì)密度,然后采用和SIMP法相同的懲罰方式。為了簡(jiǎn)單起見(jiàn),這里采用了和文獻(xiàn)相同的方式,取單元的彈性矩陣為De=[4∑l=1(Η(?el)/4)]pD(11)上式中p為罰因子,這里取為3。這樣就可以得到變化后的單元?jiǎng)偠染仃噆e=[4∑l=1(Η(?el)/4)]p∫ΩeBΤDBdΩ=[4∑l=1(Η(?el)/4)]pke0(12)式中ke0=∫ΩBTDBdΩ。同理,目標(biāo)函數(shù)的柔順性可以表示為L(zhǎng)(uΝ,?Ν)=ΝE∑e=1[4∑l=1(Η(?el)/4)]pueΤΝke0ueΝ(13)式中:?l表示基結(jié)構(gòu)有限單元節(jié)點(diǎn)上的水平集函數(shù)值;ueΝ表示單元四個(gè)節(jié)點(diǎn)的位移列陣。綜上所述,有限元離散后的拓?fù)鋬?yōu)化模型PhΤDF為MinimizeL(uN,?N)=ΝE∑e=1[4∑l=1(Η(?el)/4)]pueΤΝke0ueΝ,subjectto:ΚuΝ=fΝ,uΝ|?Ωe=u0?V(?Ν)=ΝE∑e=1∫ΩeΗ(4∑l=1Νel(x)?el)dΩ≤Vmax(14)4sigmod函數(shù)的圖像PhΤDF優(yōu)化模型中,柔順性目標(biāo)函數(shù)和體積約束的靈敏度分別為?L??m=-14ΝRE∑e=1p[4∑l=1(Η(?el)/4)]p-1δ(?m)ueΤΝke0ueΝ(15)?V??m=ΝRE∑e=1∫Ωeδ(4∑l=1Νel(x)?el)Νem(x)dΩ(16)這里NRE表示擁有節(jié)點(diǎn)m的單元總數(shù)。由于H(?)和δ(?)均為非光滑函數(shù),需要光滑后才能進(jìn)行數(shù)值結(jié)算。而且PTDF本身并沒(méi)有得到松弛,仍然是一個(gè)病態(tài)問(wèn)題,所以需要對(duì)H(?)進(jìn)行松弛。不同于文獻(xiàn)和,引進(jìn)無(wú)窮次連續(xù)可微的sigmoid函數(shù)作為松弛光滑函數(shù)Ηs(?)=11+e-α?(17)式中α為陡度參數(shù),它的值取的越大sigmoid函數(shù)就越逼近階躍函數(shù)。sigmoid函數(shù)的圖像參見(jiàn)圖1。它的圖像參見(jiàn)圖2。由圖像可知α取的越大sigmoid函數(shù)的導(dǎo)數(shù)就越逼近δ函數(shù)。sigmoid函數(shù)的導(dǎo)數(shù)為Η′s(?)=δs(?)=αeα?+e-α?+2(18)從式(5)δ函數(shù)的定義,可以看出式(15)的?L??m只有在?m=0時(shí)才不為零,也就是說(shuō)目標(biāo)函數(shù)只對(duì)拓?fù)錁?gòu)型邊界的改變敏感,這樣就非常不利于拓?fù)鋬?yōu)化的進(jìn)行。而由圖(2)又可以看出,當(dāng)α比較小的時(shí)候,H′s在較大的范圍內(nèi)不為零,而隨著α的逐漸增大,這個(gè)范圍變得越來(lái)越小,直到最后接近δ函數(shù)。因此,當(dāng)把sigmoid函數(shù)的導(dǎo)數(shù)代入式(15)中代替δ(?)就會(huì)使?L??m不為零的范圍隨著α的變化而變化。由文獻(xiàn)中關(guān)于拓?fù)鋵?dǎo)數(shù)的描述,可以看出拓?fù)鋵?dǎo)數(shù)通常與應(yīng)變能密度成正比。而式(15)又恰好是應(yīng)變能密度正比的表現(xiàn)形式。根據(jù)sigmoid函數(shù)和δ函數(shù)的性質(zhì),當(dāng)α比較小的時(shí)候,式(15)可以在相對(duì)較大的結(jié)構(gòu)范圍內(nèi)發(fā)揮拓?fù)鋵?dǎo)數(shù)的作用,而當(dāng)α比較大的時(shí)候,這個(gè)范圍就會(huì)相應(yīng)的變小。因此可以在優(yōu)化的初始階段把α取的比較小,讓較大的構(gòu)型范圍內(nèi)都可以發(fā)生拓?fù)涓淖?達(dá)到加快收斂速度的目的。隨著優(yōu)化的進(jìn)行,逐漸增大α,使能夠發(fā)生拓?fù)涓淖兊膮^(qū)域越來(lái)越小,直到收斂到最優(yōu)拓?fù)錁?gòu)型。值得注意的是,如果α始終為一個(gè)較小的值,構(gòu)型便一直可以在較大的范圍內(nèi)發(fā)生拓?fù)涓淖?不但減慢了收斂速度,而且不易收斂到最優(yōu)構(gòu)型。由于sigmoid函數(shù)具有無(wú)窮次連續(xù)可微性,因此可以相對(duì)容易地選取一種求解算法。本文采用mma數(shù)學(xué)規(guī)劃法來(lái)求解PΤDFb。5視頻綱化處理為了和其他文獻(xiàn)中的結(jié)果進(jìn)行比較,本文采用了兩個(gè)經(jīng)典的拓?fù)鋬?yōu)化算例。由于只需要得到量化結(jié)果,所有參數(shù)都進(jìn)行了無(wú)量綱化處理。彈性模量E=1,泊松比ν=0.3。此外,例1和例2均采用設(shè)計(jì)域內(nèi)水平集函數(shù)?=1所對(duì)應(yīng)的sigmoid函數(shù)值作為均勻分布初始構(gòu)型。初始型拓?fù)浣Y(jié)構(gòu)如圖3所示,懸臂梁長(zhǎng)度為8,寬度為5,厚度為1,在自由端中點(diǎn)作用一個(gè)集中力P=1。初始構(gòu)型內(nèi)劃分了80×50的有限元網(wǎng)格。體積約束Vmax=16。陡度參數(shù)α從1逐漸增加到3。最優(yōu)拓?fù)錁?gòu)型如圖4所示。表1中列出了幾種算法的比較結(jié)果。此外,當(dāng)陡度參數(shù)α始終等于1時(shí),第35步和45步的目標(biāo)函數(shù)值分別為69.92和69.89。陡度參數(shù)對(duì)分型結(jié)構(gòu)的影響如圖五所示,算例二的初始構(gòu)型和有限元離散均與算例一相同,集中力載荷P=1作用在自由端邊界下端點(diǎn)處。體積約束Vmax=12。陡度參數(shù)α從1逐漸增加到3。最優(yōu)拓?fù)錁?gòu)型如圖6所示。表2中列出了幾種算法的比較結(jié)果。需要指出的是,當(dāng)有限元網(wǎng)格取為80×50時(shí),用SIMP法計(jì)算會(huì)產(chǎn)生細(xì)小的分支結(jié)構(gòu),而且拓?fù)錁?gòu)型不是很穩(wěn)定,所以表中僅列出了第100步時(shí)所得到的結(jié)果。此外,當(dāng)陡度參數(shù)α始終等于1時(shí),第42步和52步的目標(biāo)函數(shù)值分別為126.34和126.07。從上面的兩個(gè)算例中,可以看出,陡度參數(shù)α從1增加到3可以明顯加快構(gòu)型的收斂,而α始終等于1時(shí)不但收斂變慢,而且最終很難收斂到最優(yōu)解。這和前面由拓?fù)鋵?dǎo)數(shù)概念分析所得到的