楊雨豪 鄭 偉 王英俊，2

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

2.華中科技大學(xué)數(shù)字制造裝備與技術(shù)國家重點(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]。常見的傳統(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)化無縫融合，為實(shí)現(xiàn)結(jié)構(gòu)的設(shè)計(jì)、分析、優(yōu)化的一體化奠定了理論基礎(chǔ)。

近年來，等幾何拓?fù)鋬?yōu)化得到了廣泛的研究與發(fā)展。WANG等[8]和GAO等[10]分別對(duì)等幾何拓?fù)鋬?yōu)化進(jìn)行了較為全面的綜述。GAO等[11]提出了一種具有密度分布增強(qiáng)函數(shù)的等幾何拓?fù)鋬?yōu)化方法來顯示結(jié)構(gòu)的拓?fù)錁?gòu)型，該方法可用于解決一些數(shù)值優(yōu)化問題(如多材料結(jié)構(gòu)和機(jī)械材料的合理設(shè)計(jì)等)。在通用代碼方面，GAO等[12]公開了一套高效、緊湊的等幾何拓?fù)鋬?yōu)化MATLAB代碼，為學(xué)習(xí)者提供了便利。

雖然等幾何拓?fù)鋬?yōu)化可避免網(wǎng)格離散，且具有高精度、高效率等優(yōu)點(diǎn)，但優(yōu)化中的迭代過程仍需要消耗大量時(shí)間。隨著工業(yè)的發(fā)展，實(shí)際工程應(yīng)用中的優(yōu)化問題變得日益復(fù)雜，導(dǎo)致優(yōu)化問題的計(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ā)，開展等幾何拓?fù)鋬?yōu)化高效算法的相關(guān)研究。

拓?fù)鋬?yōu)化的求解過程需要反復(fù)迭代，而每次迭代都需完成一次結(jié)構(gòu)性能分析，因此，提高拓?fù)鋬?yōu)化的效率可從減小計(jì)算量和加快優(yōu)化收斂速度兩方面開展研究。在拓?fù)鋬?yōu)化過程中，有限元分析是耗時(shí)占比較大的一個(gè)步驟，如果能夠縮減有限元方程的自由度，可有效減小計(jì)算量，且自由度規(guī)模的減小也有利于優(yōu)化迭代收斂。近年來，國內(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)化過程中，迭代收斂時(shí)所需的次數(shù)和時(shí)間也是影響整體優(yōu)化時(shí)間的重要因素。為了加速拓?fù)鋬?yōu)化問題的收斂，可在優(yōu)化迭代過程中縮減部分設(shè)計(jì)變量以加速迭代過程。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é)合來提高計(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)化的一體化。

通常使用遞歸方式來定義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)化問題中，優(yōu)化問題的數(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)化過程中的總體積；V0為設(shè)計(jì)域的總體積；VF為體積比約束(VF∈[0，1])。

最優(yōu)化準(zhǔn)則方法(OC)是求解式(5)優(yōu)化問題的一種高效方法。在求解過程中，設(shè)計(jì)變量x可以不斷迭代，并在迭代過程中逐步逼近最優(yōu)解，其迭代過程依賴于下式所示的算法：

(6)

(7)

式中，M為最大步長(zhǎng)即設(shè)計(jì)變量每次迭代的最大變化量(M>0)；η為數(shù)值衰減系數(shù)(取η=1/2)；λ為拉格朗日算子，用于在迭代過程中保證體積約束，通?？捎啥址ㄓ?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)系示意圖見圖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)化過程中位移向量的變化情況，如圖 3所示，發(fā)現(xiàn)大部分自由度對(duì)應(yīng)的位移在迭代后期幾乎保持不變。由此可知，一旦某自由度的位移變化趨于穩(wěn)定后，它在后續(xù)的迭代過程中幾乎保持不變。因此，將保持位移不變的自由度進(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為在過去的幾次迭代中位移的相對(duì)變化量；l為當(dāng)前迭代次數(shù)；U(l)為當(dāng)前迭代的位移向量；L=5，即式(14)可以衡量過去10次迭代位移變化的波動(dòng)情況。

然后，通過對(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ī)模(10×10)方程為例，闡明迭代過程中基于位移變化縮減自由度的過程：

(15)

式中，Ki，j、Ui、Fi(i，j=1，2，…，10)分別為總體剛度矩陣中第i行第j列的元素、自由度i的位移和自由度i的外載荷。

假設(shè)通過式(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ù)的變化曲線見圖3。

2.1.2基于空單元的自由度縮減

本方法的首要問題是判斷哪些自由度在等幾何分析中是可以被縮減的。以二維的等幾何拓?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減方程。雖然縮減方程無法求解U1和U31的值，但是，由式(8)可知，無論控制點(diǎn)的位移如何，空單元e1、e2、e5和e6的單元靈敏度數(shù)值都為0，即U1和U31的數(shù)值不會(huì)影響優(yōu)化求解結(jié)果，無需計(jì)算。同理，圖5中的U2、U32、U7、U37、U8和U38也無需計(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ù)的變化曲線見圖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ì)變化量來判斷是否處在迭代早期，計(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)化過程仍處于迭代早期。只有當(dāng)cobj<δ時(shí)，才可以進(jìn)行基于設(shè)計(jì)變量縮減的迭代加速。

2.2.2基于灰度抑制的收斂加速

一個(gè)單元規(guī)模為128×64的二維懸臂梁等幾何拓?fù)鋬?yōu)化結(jié)果的控制點(diǎn)密度分布如圖7所示?？梢?，即便使用了密度懲罰方法，最終的優(yōu)化結(jié)果中仍有許多中間密度存在，這將影響等幾何拓?fù)鋬?yōu)化問題的收斂速度。

圖7 控制點(diǎn)密度分布

針對(duì)此問題，參考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è)置得過大，則可能影響到優(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í)使用多種加速方法可以獲得顯著的加速效果，具體的算法流程圖見圖9。圖9中，加速流程從等幾何拓?fù)鋬?yōu)化進(jìn)入迭代開始，到迭代收斂結(jié)束。圖中黃色虛線框中是基于空單元的自由度縮減算法流程，藍(lán)色虛線框中是基于位移變化的自由度縮減算法流程，綠色框中是基于設(shè)計(jì)變量縮減和灰度抑制的收斂加速算法。流程圖中，L0是預(yù)設(shè)的迭代早期迭代次數(shù)，在經(jīng)過至少L0次迭代之后開始進(jìn)行基于位移變化的自由度縮減和收斂加速，一般設(shè)置L0>10。

圖9 高效算法流程圖

4 算例

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

二維懸臂梁是經(jīng)典拓?fù)鋬?yōu)化算例，該優(yōu)化問題的設(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)化問題的設(shè)計(jì)域、約束和外加載荷

本例中，懸臂梁的體積比約束VF=0.5，設(shè)計(jì)域劃分為三種網(wǎng)格規(guī)模：128×64、256×128、512×256。二維懸臂梁等幾何拓?fù)鋬?yōu)化問題對(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ì)取消基于單元位移的自由度縮減加速?？梢酝ㄟ^調(diào)整誤差閾值t來調(diào)整迭代后期的加速效率和時(shí)間曲線。

(a)網(wǎng)格規(guī)模128×64

由圖12可知，基于自由度縮減的加速方法可以有效加速等幾何分析求解速度，且隨著迭代的進(jìn)行，加速前后的求解時(shí)間差逐漸增大，說明基于自由度縮減的加速方法對(duì)大規(guī)模的等幾何拓?fù)鋬?yōu)化問題尤為有效。需要注意的是，基于自由度縮減的加速方法需要一些前處理，包括識(shí)別求解空單元對(duì)應(yīng)的自由度、需要縮減的位移等，但是相較于自由度縮減所帶來的加速效果，這些前處理時(shí)間的影響可以忽略。

4.1.2迭代收斂加速

本節(jié)通過對(duì)比使用迭代收斂加速算法前后，二維懸臂梁等幾何拓?fù)鋬?yōu)化問題的迭代次數(shù)、目標(biāo)函數(shù)值等來討論迭代收斂加速方法的效果。在所有算例中，優(yōu)化問題的收斂條件均設(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)化問題在收斂加速前后的迭代次數(shù)、柔度值對(duì)比

4.1.3高效等幾何拓?fù)鋬?yōu)化

綜合利用基于自由度縮減的加速方法和迭代收斂加速方法，即可在保證結(jié)果準(zhǔn)確、穩(wěn)定的基礎(chǔ)上獲得良好的加速效果。本節(jié)通過二維懸臂梁算例對(duì)比高效等幾何拓?fù)鋬?yōu)化方法與傳統(tǒng)等幾何拓?fù)鋬?yōu)化方法，來驗(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)型基本相同，說明提出的高效等幾何拓?fù)鋬?yōu)化具有可靠性。同時(shí)，得益于灰度抑制算法對(duì)中間單元的縮減，加速后的優(yōu)化結(jié)果邊界更加清晰、明顯。

(a)128×64加速前 (b)128×64加速后

不同規(guī)模下高效算法應(yīng)用前后的迭代次數(shù)、柔度和總優(yōu)化時(shí)間對(duì)比如表2所示。可見，在基于自由度縮減的加速算法和收斂加速算法的作用下，二維懸臂梁的等幾何拓?fù)鋬?yōu)化求解速度得到了較大提高?？傮w優(yōu)化的加速比分別為1.56、2.85和6.02，證明了該方法具有較好的加速效果。此外，前兩個(gè)規(guī)模下加速后的柔度值下降，第三個(gè)規(guī)模下的柔度值僅增大了0.76%，處于可接受范圍內(nèi)，說明加速方法得到的目標(biāo)函數(shù)依然準(zhǔn)確。

表2 二維懸臂梁算例高效等幾何拓?fù)鋬?yōu)化加速效率對(duì)比

綜上所述，高效等幾何拓?fù)鋬?yōu)化方法可以在保證優(yōu)化結(jié)果準(zhǔn)確的情況下大幅加速等幾何拓?fù)鋬?yōu)化問題的求解和迭代過程。

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)化方法在三維問題中的效果。

三維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é)果拓?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í)間，且隨著問題規(guī)模的擴(kuò)大，加速效果愈加顯著。

5 結(jié)論

本文提出了一種高效等幾何拓?fù)鋬?yōu)化方法，包括兩種自由度縮減方法和兩種迭代收斂加速方法。其中，基于自由度縮減的加速方法通過縮減等幾何分析方程的規(guī)模顯著縮短每次迭代中等幾何分析求解時(shí)間，迭代收斂加速方法通過加快設(shè)計(jì)變量的二值化，顯著減少優(yōu)化迭代次數(shù)。本文以二維懸臂梁和三維MBB梁兩個(gè)經(jīng)典優(yōu)化問題為例，分別對(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)化問題和大規(guī)模問題的加速效果尤其顯著，驗(yàn)證了本文提出方法的有效性。另外，高效拓?fù)鋬?yōu)化方法在處理柔度最小問題時(shí)采用了灰度抑制算法，減少了中間灰度單元的數(shù)量，可獲得更加清晰的拓?fù)錁?gòu)型。將來，該方法可以與其他拓?fù)鋬?yōu)化方法相結(jié)合，并推廣到熱、振動(dòng)、流體等其他領(lǐng)域結(jié)構(gòu)優(yōu)化問題的高效求解。