楊雨豪 鄭 偉 王英俊,2
1.華南理工大學(xué)聚合物新型成型裝備國(guó)家工程研究中心,廣州,510641
2.華中科技大學(xué)數(shù)字制造裝備與技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室,武漢,430074
拓?fù)鋬?yōu)化是一種在給定載荷與約束的條件下尋找最優(yōu)材料分布使結(jié)構(gòu)目標(biāo)性能最佳的結(jié)構(gòu)優(yōu)化方法[1]。拓?fù)鋬?yōu)化設(shè)計(jì)可在設(shè)計(jì)域內(nèi)自由改變拓?fù)浣Y(jié)構(gòu),具有高設(shè)計(jì)自由度,目前已成為結(jié)構(gòu)設(shè)計(jì)、增材制造領(lǐng)域的研究熱點(diǎn),并廣泛應(yīng)用于船舶機(jī)械、航空航天、軍工及建筑土木等領(lǐng)域[2-3]。常見(jiàn)的傳統(tǒng)拓?fù)鋬?yōu)化方法有變密度法[4]、水平集方法(level set method,LSM)[5]、漸進(jìn)優(yōu)化法(evolutionary structural optimization,ESO)[6]、移動(dòng)變形組件方法(moving morphable component,MMC)[7]和等幾何拓?fù)鋬?yōu)化[8]等。其中,等幾何拓?fù)鋬?yōu)化是一種基于等幾何分析(isogeometric analysis,IGA)[9]的新型結(jié)構(gòu)拓?fù)鋬?yōu)化方法。與其他拓?fù)鋬?yōu)化方法相比,等幾何拓?fù)鋬?yōu)化將CAD、CAE與結(jié)構(gòu)拓?fù)鋬?yōu)化無(wú)縫融合,為實(shí)現(xiàn)結(jié)構(gòu)的設(shè)計(jì)、分析、優(yōu)化的一體化奠定了理論基礎(chǔ)。
近年來(lái),等幾何拓?fù)鋬?yōu)化得到了廣泛的研究與發(fā)展。WANG等[8]和GAO等[10]分別對(duì)等幾何拓?fù)鋬?yōu)化進(jìn)行了較為全面的綜述。GAO等[11]提出了一種具有密度分布增強(qiáng)函數(shù)的等幾何拓?fù)鋬?yōu)化方法來(lái)顯示結(jié)構(gòu)的拓?fù)錁?gòu)型,該方法可用于解決一些數(shù)值優(yōu)化問(wèn)題(如多材料結(jié)構(gòu)和機(jī)械材料的合理設(shè)計(jì)等)。在通用代碼方面,GAO等[12]公開(kāi)了一套高效、緊湊的等幾何拓?fù)鋬?yōu)化MATLAB代碼,為學(xué)習(xí)者提供了便利。
雖然等幾何拓?fù)鋬?yōu)化可避免網(wǎng)格離散,且具有高精度、高效率等優(yōu)點(diǎn),但優(yōu)化中的迭代過(guò)程仍需要消耗大量時(shí)間。隨著工業(yè)的發(fā)展,實(shí)際工程應(yīng)用中的優(yōu)化問(wèn)題變得日益復(fù)雜,導(dǎo)致優(yōu)化問(wèn)題的計(jì)算規(guī)模也相應(yīng)增大,因此,提高計(jì)算效率成為相關(guān)研究人員重點(diǎn)關(guān)注的內(nèi)容。目前主要從以下兩個(gè)方面提高拓?fù)鋬?yōu)化的計(jì)算效率:一是從硬件角度,引入CPU/GPU并行計(jì)算[13-14]以提高運(yùn)算速度;二是從高效算法角度。從算法層面提高優(yōu)化計(jì)算效率可從根本上減小拓?fù)鋬?yōu)化的時(shí)間復(fù)雜度[15],因此,本文主要從算法角度出發(fā),開(kāi)展等幾何拓?fù)鋬?yōu)化高效算法的相關(guān)研究。
拓?fù)鋬?yōu)化的求解過(guò)程需要反復(fù)迭代,而每次迭代都需完成一次結(jié)構(gòu)性能分析,因此,提高拓?fù)鋬?yōu)化的效率可從減小計(jì)算量和加快優(yōu)化收斂速度兩方面開(kāi)展研究。在拓?fù)鋬?yōu)化過(guò)程中,有限元分析是耗時(shí)占比較大的一個(gè)步驟,如果能夠縮減有限元方程的自由度,可有效減小計(jì)算量,且自由度規(guī)模的減小也有利于優(yōu)化迭代收斂。近年來(lái),國(guó)內(nèi)外學(xué)者提出了一系列用于縮減拓?fù)鋬?yōu)化中有限元方程自由度的方法。WANG等[16]實(shí)現(xiàn)了基于多層網(wǎng)格的等幾何拓?fù)鋬?yōu)化有限元求解,有效減小了有限元方程的規(guī)模;杜義賢等[17]將序列插值模型和多重網(wǎng)格方法應(yīng)用于拓?fù)鋬?yōu)化,縮減了有限元求解的自由度;GOGU[18]提出了一種動(dòng)態(tài)構(gòu)造簡(jiǎn)化基的方法,有效縮短有限元求解步驟的耗時(shí)。
在拓?fù)鋬?yōu)化過(guò)程中,迭代收斂時(shí)所需的次數(shù)和時(shí)間也是影響整體優(yōu)化時(shí)間的重要因素。為了加速拓?fù)鋬?yōu)化問(wèn)題的收斂,可在優(yōu)化迭代過(guò)程中縮減部分設(shè)計(jì)變量以加速迭代過(guò)程。KIM等[19]提出一種設(shè)計(jì)變量可縮減方法用于減少拓?fù)鋬?yōu)化迭代次數(shù),有效縮短迭代收斂所需時(shí)間。LIAO等[20]提出了一種局部更新策略,只更新實(shí)體區(qū)域及其鄰近邊界區(qū)域,在減小更新設(shè)計(jì)變量所需計(jì)算量的同時(shí)減少了優(yōu)化迭代次數(shù)。除了減少設(shè)計(jì)變量以加快收斂速度外,還有一些方法可以加速迭代收斂,如廉睿超等[21]提出一種灰度單元雙重懲罰方法,加速了中間密度單元的二值化;許小奎等[22]提出了兩種形式的對(duì)比度增強(qiáng)算子,驅(qū)動(dòng)中間密度向兩級(jí)進(jìn)行分化。此外,拓?fù)鋬?yōu)化也可以和機(jī)器學(xué)習(xí)等熱門方法相結(jié)合來(lái)提高計(jì)算效率[23-24]。
本文從性能分析方程的自由度縮減和優(yōu)化收斂加速兩個(gè)方面提高等幾何拓?fù)鋬?yōu)化的計(jì)算效率。
等幾何分析將NURBS引入CAE描述分析物理場(chǎng),實(shí)現(xiàn)了CAD建模與CAE分析的統(tǒng)一表達(dá)。在拓?fù)鋬?yōu)化領(lǐng)域,用NURBS基函數(shù)插值優(yōu)化變量場(chǎng),可進(jìn)一步將CAD建模、CAE分析與結(jié)構(gòu)拓?fù)鋬?yōu)化納入統(tǒng)一框架,有望實(shí)現(xiàn)結(jié)構(gòu)設(shè)計(jì)、分析與優(yōu)化的一體化。
通常使用遞歸方式來(lái)定義B樣條的基函數(shù)。Ξ=(ξ1,ξ2,…,ξn+p+1)是一個(gè)參數(shù)空間中的非遞減實(shí)數(shù)數(shù)列,被稱作節(jié)點(diǎn)向量,其中,n是控制點(diǎn)和基函數(shù)的數(shù)量,p是樣條的次數(shù)。由Cox-de Boor遞推公式[25]定義p次B樣條基函數(shù):
(1)
其中,Bi,p(ξ)為第i個(gè)p次B樣條基函數(shù)。向每個(gè)B樣條基函數(shù)引入一個(gè)權(quán)重ωi,便可得到單變量的NURBS基函數(shù):
(2)
根據(jù)張量積性質(zhì),可以構(gòu)造二維或多維的NURBS基函數(shù)Ni,p(ξ)如下:
(3)
p=(p1,p2,…,pdp)
式中,ξ為參數(shù)坐標(biāo);i為張量積結(jié)構(gòu)中的索引位置;p為所有參數(shù)方向的多項(xiàng)式次數(shù);Nim,pm為沿參數(shù)方向m的B樣條基函數(shù),可由式(2)計(jì)算得到;dp為參數(shù)空間的維度;Ξ(m)為該參數(shù)空間內(nèi)的節(jié)點(diǎn)向量,表示參數(shù)坐標(biāo);nm為對(duì)應(yīng)的函數(shù)個(gè)數(shù);pm為沿參數(shù)方向m的多項(xiàng)式次數(shù)。
固體各向同性材料懲罰(SIMP)法是密度法拓?fù)鋬?yōu)化中最為常用的材料插值方法。在SIMP方法中,每個(gè)單元被賦予一個(gè)介于0、1之間的單元密度,單元密度xe和彈性模量Ee之間的關(guān)系如下:
(4)
式中,Emin為彈性模量最小值(Emin>0),給空單元賦予該最小值可防止剛度矩陣發(fā)生奇異;E0為單元密度為1的實(shí)心單元的彈性模量;β為懲罰系數(shù),用于實(shí)現(xiàn)中間密度的二值化,通常取β= 3。
在經(jīng)典最小柔度目標(biāo)優(yōu)化問(wèn)題中,優(yōu)化問(wèn)題的數(shù)學(xué)模型如下:
(5)
式中,C(x)為目標(biāo)函數(shù),此處表示結(jié)構(gòu)柔度(是實(shí)際柔度值的兩倍);x為設(shè)計(jì)變量;n為設(shè)計(jì)域中設(shè)計(jì)變量的個(gè)數(shù);U為總體位移向量;F為總體受力(載荷)向量;K為總體剛度矩陣;ue為單元e的位移向量;ke為單元e的單元?jiǎng)偠染仃?;V(x)為優(yōu)化過(guò)程中的總體積;V0為設(shè)計(jì)域的總體積;VF為體積比約束(VF∈[0,1])。
最優(yōu)化準(zhǔn)則方法(OC)是求解式(5)優(yōu)化問(wèn)題的一種高效方法。在求解過(guò)程中,設(shè)計(jì)變量x可以不斷迭代,并在迭代過(guò)程中逐步逼近最優(yōu)解,其迭代過(guò)程依賴于下式所示的算法:
(6)
(7)
式中,M為最大步長(zhǎng)即設(shè)計(jì)變量每次迭代的最大變化量(M>0);η為數(shù)值衰減系數(shù)(取η=1/2);λ為拉格朗日算子,用于在迭代過(guò)程中保證體積約束,通??捎啥址ㄓ?jì)算獲得[26]。
目標(biāo)函數(shù)C和體積V對(duì)xe的導(dǎo)數(shù)分別為
(8)
(9)
其中,式(9)建立的前提是每個(gè)單元都有單位體積。
與傳統(tǒng)的SIMP方法不同,等幾何拓?fù)鋬?yōu)化的數(shù)值計(jì)算是在控制點(diǎn)上進(jìn)行的,因此,基于等參的思想,在等幾何拓?fù)鋬?yōu)化中參數(shù)坐標(biāo)為ξ的設(shè)計(jì)變量x可由附近的控制點(diǎn)求出:
(10)
式中,Ni為影響該參數(shù)坐標(biāo)的第i個(gè)控制點(diǎn)的NURBS基函數(shù);xi為控制點(diǎn)密度。
在基于密度法的等幾何拓?fù)鋬?yōu)化中,單元的剛度矩陣ke可由下式計(jì)算:
(11)
在等幾何拓?fù)鋬?yōu)化中,控制點(diǎn)密度x(ic)可由下式計(jì)算:
(12)
式中,xi為設(shè)計(jì)域中編號(hào)為i的控制點(diǎn)的密度值;ic對(duì)應(yīng)編號(hào)為i的單元的中心;ci為影響單元i的所有控制點(diǎn);cij為影響單元i的第j個(gè)控制點(diǎn);Ncij(ic)為控制點(diǎn)cij的NURBS基函數(shù)。
單元i及對(duì)應(yīng)的控制點(diǎn)關(guān)系如圖1所示,在次數(shù)為2的二維控制點(diǎn)網(wǎng)格中,單元i受到周圍9個(gè)控制點(diǎn)(圖中紅色控制點(diǎn))的影響。
圖1 單元i及其控制點(diǎn)的對(duì)應(yīng)關(guān)系示意圖
因此,傳統(tǒng)SIMP法中求解靈敏度的式(8)可以改寫為
(13)
式中,ei為被控制點(diǎn)i影響的所有單元;eij為ei中的第j個(gè)單元。
一個(gè)關(guān)于控制點(diǎn)i及所控制的單元的位置關(guān)系示意圖見(jiàn)圖2,圖中以二維控制點(diǎn)網(wǎng)絡(luò)為例,紅色控制點(diǎn)i的影響范圍是它周圍的9個(gè)單元。
圖2 控制點(diǎn)i及其所控制的單元示意圖
拓?fù)鋬?yōu)化中的關(guān)鍵一步是求解式(5)中的方程KU=F,并得到位移向量U。本節(jié)介紹兩種自由度縮減方法以加速等幾何拓?fù)鋬?yōu)化中的控制點(diǎn)位移求解,分別是基于位移變化的自由度縮減和基于空單元的自由度縮減。
2.1.1基于位移變化的自由度縮減
以單元規(guī)模為64×32的懸臂梁等幾何拓?fù)鋬?yōu)化為例,觀察優(yōu)化過(guò)程中位移向量的變化情況,如圖 3所示,發(fā)現(xiàn)大部分自由度對(duì)應(yīng)的位移在迭代后期幾乎保持不變。由此可知,一旦某自由度的位移變化趨于穩(wěn)定后,它在后續(xù)的迭代過(guò)程中幾乎保持不變。因此,將保持位移不變的自由度進(jìn)行縮減,能減小性能分析方程的求解規(guī)模,從而加快優(yōu)化進(jìn)程。
圖3 隨機(jī)選取的自由度位移隨迭代次數(shù)變化曲線
基于位移變化的自由度縮減可分為兩個(gè)步驟:①判斷各個(gè)自由度的位移是否已經(jīng)保持穩(wěn)定;②對(duì)位移已穩(wěn)定的自由度進(jìn)行縮減。
首先,本文提出了一種多次迭代內(nèi)自由度位移相對(duì)變化的量化表示:
(14)
式中,D為在過(guò)去的幾次迭代中位移的相對(duì)變化量;l為當(dāng)前迭代次數(shù);U(l)為當(dāng)前迭代的位移向量;L=5,即式(14)可以衡量過(guò)去10次迭代位移變化的波動(dòng)情況。
然后,通過(guò)對(duì)比D與設(shè)置的閾值ε的大小,即可判斷某個(gè)自由度對(duì)應(yīng)的位移是否已經(jīng)保持穩(wěn)定,即當(dāng)D(i)<ε時(shí),可認(rèn)為在i處的位移已經(jīng)保持穩(wěn)定,可進(jìn)行自由度縮減,反之不能。通常情況下,隨著迭代的進(jìn)行,可縮減的自由度會(huì)逐漸增加。
一旦判斷某自由度處的位移在數(shù)次迭代處于穩(wěn)定狀態(tài),即可在方程求解過(guò)程中移除該自由度。下面以小規(guī)模(10×10)方程為例,闡明迭代過(guò)程中基于位移變化縮減自由度的過(guò)程:
(15)
式中,Ki,j、Ui、Fi(i,j=1,2,…,10)分別為總體剛度矩陣中第i行第j列的元素、自由度i的位移和自由度i的外載荷。
假設(shè)通過(guò)式(14)判斷U9在最近的數(shù)次迭代中基本保持不變,對(duì)上面的分析方程進(jìn)行自由度縮減,可得到縮減規(guī)模(9×9)的方程:
(16)
基于位移變化的自由度縮減流程中,部分自由度位移被直接代入下一次迭代計(jì)算中,可能會(huì)引入一定誤差。隨著迭代的進(jìn)行,累計(jì)誤差可能會(huì)大到足以影響等幾何拓?fù)鋬?yōu)化結(jié)果,因此本文引入一種誤差的量化表示方法:
(17)
式中,er為等幾何分析求解誤差;K、F為自由度縮減之前的全局剛度矩陣和力向量;U*為規(guī)??s減后等式的求解結(jié)果位移向量。
當(dāng)er大于設(shè)定的誤差閾值時(shí),下次迭代將不再進(jìn)行基于位移變化的自由度縮減,以保證等幾何分析求解的誤差在設(shè)定閾值內(nèi)。隨機(jī)選取的自由度位移隨迭代次數(shù)的變化曲線見(jiàn)圖3。
2.1.2基于空單元的自由度縮減
本方法的首要問(wèn)題是判斷哪些自由度在等幾何分析中是可以被縮減的。以二維的等幾何拓?fù)鋬?yōu)化為例,圖4所示為一個(gè)由12個(gè)矩形單元、30個(gè)控制點(diǎn)和60個(gè)自由度描述(前30個(gè)為各控制點(diǎn)橫向自由度,后30個(gè)為各控制點(diǎn)縱向自由度)的等幾何拓?fù)鋬?yōu)化設(shè)計(jì)域。圖4中,每個(gè)單元的編號(hào)ei寫在單元內(nèi)部上方,每個(gè)控制點(diǎn)的編號(hào)寫在控制點(diǎn)的上方,括號(hào)中是該控制點(diǎn)對(duì)應(yīng)的兩個(gè)自由度編號(hào)。該設(shè)計(jì)域?qū)?yīng)的等幾何分析方程為
圖4 由12個(gè)矩形單元、30個(gè)控制點(diǎn)和60個(gè)自由度描述的二維等幾何設(shè)計(jì)域
(18)
基于SIMP法的等幾何拓?fù)鋬?yōu)化中,控制點(diǎn)的密度會(huì)隨著迭代的進(jìn)行向0或1靠近。假設(shè)在某次迭代之后,圖4中影響e1、e2、e5和e6單元的所有控制點(diǎn)的控制點(diǎn)密度均更新為0,如圖5所示。如果忽略賦給空單元的非常小的彈性模量Emin,那么單元e1、e2、e5和e6的單元?jiǎng)偠染仃嚳梢暈榱憔仃?,控制點(diǎn)1影響的所有單元均為空單元。根據(jù)總體剛度矩陣的組裝方法可推出,整體剛度矩陣中控制點(diǎn)1所對(duì)應(yīng)的第1、31行和1、31列的元素均為0,同理,控制點(diǎn)2、7、8對(duì)應(yīng)的2、32、7、37、8、38自由度所在行與列的元素也為0。僅考慮1、31自由度對(duì)應(yīng)的等幾何分析方程情況,可以將式(18)寫為
圖5 可縮減控制點(diǎn)示意圖
(19)
鑒于總體剛度矩陣的第1、31行和1、31列均為零元素,可將這些行、列移除,同時(shí)移除U和F中的對(duì)應(yīng)元素,從而得到結(jié)果與原方程相同的規(guī)??s減方程。雖然縮減方程無(wú)法求解U1和U31的值,但是,由式(8)可知,無(wú)論控制點(diǎn)的位移如何,空單元e1、e2、e5和e6的單元靈敏度數(shù)值都為0,即U1和U31的數(shù)值不會(huì)影響優(yōu)化求解結(jié)果,無(wú)需計(jì)算。同理,圖5中的U2、U32、U7、U37、U8和U38也無(wú)需計(jì)算。因此,當(dāng)一個(gè)控制點(diǎn)影響的所有單元均為空單元時(shí),該控制點(diǎn)對(duì)應(yīng)的自由度可以縮減。隨著迭代的進(jìn)行,將出現(xiàn)更多的空單元,即有更多的自由度可以縮減。
2.2.1基于設(shè)計(jì)變量縮減的收斂加速
受拓?fù)鋬?yōu)化加速算法研究[16,19-20]的啟發(fā),本文引入了一種基于設(shè)計(jì)變量縮減的收斂加速方法,具體的收斂策略是:在等幾何拓?fù)鋬?yōu)化迭代中,當(dāng)某設(shè)計(jì)變量連續(xù)多次迭代保持基本不變時(shí),可以停止更新該設(shè)計(jì)變量,以加速收斂。
以單元規(guī)模為64×32的二維懸臂梁等幾何拓?fù)鋬?yōu)化為例,隨機(jī)選取的幾個(gè)控制點(diǎn)處的控制點(diǎn)密度隨迭代次數(shù)的變化曲線見(jiàn)圖6。從圖6中可以看出,大多數(shù)控制點(diǎn)密度在接近0或1后即保持基本不變,因此在之后的迭代中可以不更新這些已保持穩(wěn)定的設(shè)計(jì)變量。
圖6 二維等幾何拓?fù)鋬?yōu)化在64×32規(guī)模下的幾個(gè)控制點(diǎn)密度變化曲線
判斷控制點(diǎn)密度是否需要更新的標(biāo)準(zhǔn)如下:
max(x(k),x(k-1),…,x(k-Q+2),x(k-Q+1))<0.01
或
min(x(k),x(k-1),…,x(k-Q+2),x(k-Q+1))>0.99
(20)
式中,x(k)為第k次迭代的設(shè)計(jì)變量值;Q為一個(gè)整數(shù),本文將Q設(shè)置為5。
式(20)的判斷條件表示,當(dāng)某控制點(diǎn)的密度連續(xù)Q次迭代保持大于0.99或小于0.01時(shí),即可在之后的迭代中不更新該變量。另外,對(duì)于不再更新的變量,將控制點(diǎn)密度大于0.99的賦值為1,而將控制點(diǎn)密度小于0.01的賦值為0,以加速收斂。
由圖6可發(fā)現(xiàn),索引為755和777的設(shè)計(jì)變量可能在迭代早期尚未穩(wěn)定,為保證基于設(shè)計(jì)變量縮減的收斂加速方法的準(zhǔn)確性,應(yīng)在迭代的中后期使用該收斂加速方法。可以利用目標(biāo)函數(shù)值的相對(duì)變化量來(lái)判斷是否處在迭代早期,計(jì)算公式如下:
(21)
式中,cobj為目標(biāo)函數(shù)的相對(duì)變化量;c(k)為第k次迭代的目標(biāo)函數(shù)值;N為一個(gè)整數(shù),本文將N設(shè)置為5。
當(dāng)目標(biāo)函數(shù)的相對(duì)變化量cobj大于設(shè)定的閾值δ時(shí),可以認(rèn)為優(yōu)化過(guò)程仍處于迭代早期。只有當(dāng)cobj<δ時(shí),才可以進(jìn)行基于設(shè)計(jì)變量縮減的迭代加速。
2.2.2基于灰度抑制的收斂加速
一個(gè)單元規(guī)模為128×64的二維懸臂梁等幾何拓?fù)鋬?yōu)化結(jié)果的控制點(diǎn)密度分布如圖7所示。可見(jiàn),即便使用了密度懲罰方法,最終的優(yōu)化結(jié)果中仍有許多中間密度存在,這將影響等幾何拓?fù)鋬?yōu)化問(wèn)題的收斂速度。
圖7 控制點(diǎn)密度分布
針對(duì)此問(wèn)題,參考GROENWOLD等[27]的研究,本文引入了基于灰度抑制的設(shè)計(jì)變量更新方法,在最優(yōu)化準(zhǔn)則方法的基礎(chǔ)上繼續(xù)進(jìn)行灰度抑制。將sigmoid函數(shù)向x軸正向平移0.5單位長(zhǎng)度可得
(22)
圖8 不同參數(shù)a對(duì)應(yīng)的平移后的sigmoid函數(shù)圖像
本文中,迭代中的參數(shù)a由下式求得:
a=t/cobj
(23)
其中,t是一個(gè)常數(shù)。在等幾何拓?fù)鋬?yōu)化中,隨著迭代的進(jìn)行,cobj將逐漸減小,a逐漸增大,即加快設(shè)計(jì)變量二值化。如果常數(shù)t設(shè)置得太小,則加速不明顯;如果t設(shè)置得過(guò)大,則可能影響到優(yōu)化結(jié)果的準(zhǔn)確性。經(jīng)測(cè)試,本文將t設(shè)置為6。由于拓?fù)鋬?yōu)化迭代早期的拓?fù)錁?gòu)型還不穩(wěn)定,應(yīng)根據(jù)式(21)判斷何時(shí)引入基于灰度抑制的收斂加速方法。
本文闡述了兩種等幾何拓?fù)鋬?yōu)化加速方法,分別是自由度縮減方法和迭代收斂加速方法,在等幾何拓?fù)鋬?yōu)化中同時(shí)使用多種加速方法可以獲得顯著的加速效果,具體的算法流程圖見(jiàn)圖9。圖9中,加速流程從等幾何拓?fù)鋬?yōu)化進(jìn)入迭代開(kāi)始,到迭代收斂結(jié)束。圖中黃色虛線框中是基于空單元的自由度縮減算法流程,藍(lán)色虛線框中是基于位移變化的自由度縮減算法流程,綠色框中是基于設(shè)計(jì)變量縮減和灰度抑制的收斂加速算法。流程圖中,L0是預(yù)設(shè)的迭代早期迭代次數(shù),在經(jīng)過(guò)至少L0次迭代之后開(kāi)始進(jìn)行基于位移變化的自由度縮減和收斂加速,一般設(shè)置L0>10。
圖9 高效算法流程圖
二維懸臂梁是經(jīng)典拓?fù)鋬?yōu)化算例,該優(yōu)化問(wèn)題的設(shè)計(jì)域、邊界條件和外加載荷如圖10所示。本算例的優(yōu)化目標(biāo)是在保證體積約束的條件下,求得柔度最小的拓?fù)錁?gòu)型。本文在三個(gè)不同控制點(diǎn)網(wǎng)格規(guī)模下,分別對(duì)比加速前、自由度縮減后、收斂加速后和整體高效方法的優(yōu)化時(shí)間與優(yōu)化結(jié)果,進(jìn)而研究各個(gè)算法的加速效果。
圖10 二維懸臂梁優(yōu)化問(wèn)題的設(shè)計(jì)域、約束和外加載荷
本例中,懸臂梁的體積比約束VF=0.5,設(shè)計(jì)域劃分為三種網(wǎng)格規(guī)模:128×64、256×128、512×256。二維懸臂梁等幾何拓?fù)鋬?yōu)化問(wèn)題對(duì)應(yīng)的控制點(diǎn)網(wǎng)格如圖11所示,以10×5的單元網(wǎng)格(12×7的控制點(diǎn)網(wǎng)格)為例,控制點(diǎn)用綠色表示。
圖11 二維懸臂梁等幾何拓?fù)鋬?yōu)化設(shè)計(jì)域的控制點(diǎn)網(wǎng)格
4.1.1基于自由度縮減的加速
本例中,自由度縮減算法的參數(shù)設(shè)置如下:誤差閾值t= 0.04,位移變化閾值ε=0.0004。對(duì)比三種網(wǎng)格規(guī)模下等幾何拓?fù)鋬?yōu)化的等幾何分析求解時(shí)間,前100次迭代的等幾何分析求解時(shí)間-迭代次數(shù)曲線如圖12所示。在迭代早期兩條曲線基本重合,當(dāng)引入基于自由度縮減的加速方法后,求解時(shí)間明顯縮短,而且隨著網(wǎng)格規(guī)模的增大,加速效果也愈加顯著。圖中藍(lán)色曲線在迭代后期出現(xiàn)一些波動(dòng),這是由于觸發(fā)了誤差閾值限制,當(dāng)誤差值大于設(shè)定的閾值時(shí),下一次迭代會(huì)取消基于單元位移的自由度縮減加速??梢酝ㄟ^(guò)調(diào)整誤差閾值t來(lái)調(diào)整迭代后期的加速效率和時(shí)間曲線。
(a)網(wǎng)格規(guī)模128×64
由圖12可知,基于自由度縮減的加速方法可以有效加速等幾何分析求解速度,且隨著迭代的進(jìn)行,加速前后的求解時(shí)間差逐漸增大,說(shuō)明基于自由度縮減的加速方法對(duì)大規(guī)模的等幾何拓?fù)鋬?yōu)化問(wèn)題尤為有效。需要注意的是,基于自由度縮減的加速方法需要一些前處理,包括識(shí)別求解空單元對(duì)應(yīng)的自由度、需要縮減的位移等,但是相較于自由度縮減所帶來(lái)的加速效果,這些前處理時(shí)間的影響可以忽略。
4.1.2迭代收斂加速
本節(jié)通過(guò)對(duì)比使用迭代收斂加速算法前后,二維懸臂梁等幾何拓?fù)鋬?yōu)化問(wèn)題的迭代次數(shù)、目標(biāo)函數(shù)值等來(lái)討論迭代收斂加速方法的效果。在所有算例中,優(yōu)化問(wèn)題的收斂條件均設(shè)置為:連續(xù)兩次迭代的控制點(diǎn)密度的最大差值小于0.01。
迭代收斂加速的參數(shù)設(shè)置如下:目標(biāo)函數(shù)改變的閾值δ=0.01,灰度抑制算法的參數(shù)a=6/cobj。在三種網(wǎng)格規(guī)模下,算例在收斂加速前后的迭代次數(shù)、目標(biāo)函數(shù)值的對(duì)比如表1所示。由表1可知,迭代收斂加速算法可以顯著加速等幾何拓?fù)鋬?yōu)化進(jìn)程,減少大量的迭代次數(shù),減小冗余計(jì)算量。得益于灰度抑制算法,加速后結(jié)果的柔度通常優(yōu)于加速前的優(yōu)化結(jié)果。
表1 二維懸臂梁等幾何拓?fù)鋬?yōu)化問(wèn)題在收斂加速前后的迭代次數(shù)、柔度值對(duì)比
4.1.3高效等幾何拓?fù)鋬?yōu)化
綜合利用基于自由度縮減的加速方法和迭代收斂加速方法,即可在保證結(jié)果準(zhǔn)確、穩(wěn)定的基礎(chǔ)上獲得良好的加速效果。本節(jié)通過(guò)二維懸臂梁算例對(duì)比高效等幾何拓?fù)鋬?yōu)化方法與傳統(tǒng)等幾何拓?fù)鋬?yōu)化方法,來(lái)驗(yàn)證高效等幾何拓?fù)鋬?yōu)化方法的高效性。
高效等幾何拓?fù)鋬?yōu)化算法的參數(shù)設(shè)置與前述相同,以128×64、256×128、512×256三種網(wǎng)格規(guī)模劃分設(shè)計(jì)域。圖13所示為不同規(guī)模下使用高效等幾何拓?fù)鋬?yōu)化算法前后的優(yōu)化結(jié)果模型,可以看出,使用加速算法前后,拓?fù)錁?gòu)型基本相同,說(shuō)明提出的高效等幾何拓?fù)鋬?yōu)化具有可靠性。同時(shí),得益于灰度抑制算法對(duì)中間單元的縮減,加速后的優(yōu)化結(jié)果邊界更加清晰、明顯。
(a)128×64加速前 (b)128×64加速后
不同規(guī)模下高效算法應(yīng)用前后的迭代次數(shù)、柔度和總優(yōu)化時(shí)間對(duì)比如表2所示??梢?jiàn),在基于自由度縮減的加速算法和收斂加速算法的作用下,二維懸臂梁的等幾何拓?fù)鋬?yōu)化求解速度得到了較大提高??傮w優(yōu)化的加速比分別為1.56、2.85和6.02,證明了該方法具有較好的加速效果。此外,前兩個(gè)規(guī)模下加速后的柔度值下降,第三個(gè)規(guī)模下的柔度值僅增大了0.76%,處于可接受范圍內(nèi),說(shuō)明加速方法得到的目標(biāo)函數(shù)依然準(zhǔn)確。
表2 二維懸臂梁算例高效等幾何拓?fù)鋬?yōu)化加速效率對(duì)比
綜上所述,高效等幾何拓?fù)鋬?yōu)化方法可以在保證優(yōu)化結(jié)果準(zhǔn)確的情況下大幅加速等幾何拓?fù)鋬?yōu)化問(wèn)題的求解和迭代過(guò)程。
MBB梁也是拓?fù)鋬?yōu)化的經(jīng)典基準(zhǔn)算例,取它一半的設(shè)計(jì)域、邊界條件和外加載荷如圖14a所示。圖中,左側(cè)面的X、Z方向受位移約束,右下方的活動(dòng)鉸支座表示其Y、Z方向受位移約束。以體積比約束VF= 0.3的三維MBB梁等幾何拓?fù)鋬?yōu)化為例,驗(yàn)證自由度縮減和收斂加速的高效等幾何拓?fù)鋬?yōu)化方法在三維問(wèn)題中的效果。
三維MBB梁等幾何拓?fù)鋬?yōu)化的目標(biāo)是在保證體積約束的情況下獲得柔度最小的拓?fù)浣Y(jié)構(gòu),設(shè)計(jì)域?qū)?yīng)的控制點(diǎn)網(wǎng)格如圖14b所示,為了展示方便,圖中的網(wǎng)格規(guī)模被縮小為16×9×4,綠色為控制點(diǎn),紅色為單元網(wǎng)格線。
(a)設(shè)計(jì)域與邊界條件
本節(jié)在三種網(wǎng)格規(guī)模上(30×10×4,60×20×6和90×30×8)對(duì)比驗(yàn)證加速前后迭代次數(shù)、優(yōu)化結(jié)果和優(yōu)化時(shí)間。本例中,高效等幾何拓?fù)鋬?yōu)化方法的參數(shù)設(shè)置為:自由度縮減的誤差閾值t= 0.04、位移變化閾值ε=0.005、目標(biāo)函數(shù)改變的閾值δ=0.3、灰度抑制算法的參數(shù)a=6/cobj。如表3所示,加速后的優(yōu)化迭代次數(shù)明顯減少,目標(biāo)函數(shù)柔度值也有所減小,整體迭代優(yōu)化時(shí)間明顯縮短,不同規(guī)模下的加速比分別為2.03、2.69和5.64。理論上,高效等幾何拓?fù)鋬?yōu)化算法在更大的網(wǎng)格規(guī)模下會(huì)獲得更好的加速效果。
表3 三維MBB梁算例高效等幾何拓?fù)鋬?yōu)化加速效率對(duì)比
不同規(guī)模下加速前后的三維MBB梁等幾何拓?fù)鋬?yōu)化結(jié)果構(gòu)型如圖15所示??梢?jiàn),加速前后的結(jié)果拓?fù)錁?gòu)型基本相同,因此高效等幾何拓?fù)鋬?yōu)化算法在加速計(jì)算的同時(shí)不影響最終的優(yōu)化構(gòu)型。由于引入了灰度抑制算法,使得圖15右側(cè)的加速后結(jié)果構(gòu)型的中間灰度單元更少、邊界更加清晰。
(a)30×10×4加速前 (b)30×10×4加速后
綜上所述,本文所述的高效等幾何拓?fù)鋬?yōu)化方法在三維算例中依然有著較好表現(xiàn),可大幅減少優(yōu)化迭代次數(shù)、縮短優(yōu)化時(shí)間,且隨著問(wèn)題規(guī)模的擴(kuò)大,加速效果愈加顯著。
本文提出了一種高效等幾何拓?fù)鋬?yōu)化方法,包括兩種自由度縮減方法和兩種迭代收斂加速方法。其中,基于自由度縮減的加速方法通過(guò)縮減等幾何分析方程的規(guī)模顯著縮短每次迭代中等幾何分析求解時(shí)間,迭代收斂加速方法通過(guò)加快設(shè)計(jì)變量的二值化,顯著減少優(yōu)化迭代次數(shù)。本文以二維懸臂梁和三維MBB梁兩個(gè)經(jīng)典優(yōu)化問(wèn)題為例,分別對(duì)比了迭代次數(shù)、目標(biāo)函數(shù)值和迭代時(shí)間,驗(yàn)證了高效等幾何拓?fù)鋬?yōu)化的加速效果及其對(duì)優(yōu)化結(jié)果的影響。結(jié)果表明,本文提出的高效等幾何拓?fù)鋬?yōu)化方法可以顯著加速等幾何拓?fù)鋬?yōu)化,計(jì)算時(shí)間成本大幅縮減,對(duì)難以收斂的優(yōu)化問(wèn)題和大規(guī)模問(wèn)題的加速效果尤其顯著,驗(yàn)證了本文提出方法的有效性。另外,高效拓?fù)鋬?yōu)化方法在處理柔度最小問(wèn)題時(shí)采用了灰度抑制算法,減少了中間灰度單元的數(shù)量,可獲得更加清晰的拓?fù)錁?gòu)型。將來(lái),該方法可以與其他拓?fù)鋬?yōu)化方法相結(jié)合,并推廣到熱、振動(dòng)、流體等其他領(lǐng)域結(jié)構(gòu)優(yōu)化問(wèn)題的高效求解。