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

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