一種自由度縮減和收斂加速的高效等幾何拓?fù)鋬?yōu)化方法

楊雨豪 鄭 偉 王英俊，2

1.華南理工大學(xué)聚合物新型成型裝備國(guó)家工程研究中心，廣州，510641

2.華中科技大學(xué)數(shù)字制造裝備與技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢，430074

0 引言

拓?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ì)算效率。

1 理論基礎(chǔ)

1.1 NURBS基函數(shù)

等幾何分析將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ù)。

1.2 基于SIMP方法的拓?fù)鋬?yōu)化

固體各向同性材料懲罰(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è)單元都有單位體積。

1.3 基于控制點(diǎn)的等幾何拓?fù)鋬?yōu)化

與傳統(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及其所控制的單元示意圖

2 高效等幾何拓?fù)鋬?yōu)化方法

2.1 自由度縮減方法

拓?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 收斂加速方法

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í)引入基于灰度抑制的收斂加速方法。

3 多重加速流程的算法實(shí)現(xiàn)

本文闡述了兩種等幾何拓?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 高效算法流程圖

4 算例

4.1 二維懸臂梁等幾何拓?fù)鋬?yōu)化算例

二維懸臂梁是經(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ò)程。

4.2 三維MBB梁等幾何拓?fù)鋬?yōu)化算例

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ò)大，加速效果愈加顯著。

5 結(jié)論

本文提出了一種高效等幾何拓?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)題的高效求解。