魏鵬 蔣子潤(rùn)

(1.華南理工大學(xué) 土木與交通學(xué)院∥亞熱帶建筑科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，廣東 廣州 510640；2.廣州大學(xué) 廣東省地震工程與應(yīng)用技術(shù)重點(diǎn)實(shí)驗(yàn)室，廣東 廣州 510006)

至今為止，變密度法[13-14]仍是應(yīng)用最為廣泛的拓?fù)鋬?yōu)化方法，常被用來求解柔順度最小化問題[15]、柔順機(jī)構(gòu)優(yōu)化問題[16]，可靠性優(yōu)化問題[17]等。變密度法采用0-1之間連續(xù)分布的單元密度表示結(jié)構(gòu)中的材料分布，0和1分別表示孔洞和實(shí)體材料，優(yōu)化過程中中間密度逐漸趨近于0和1，使材料分布最終收斂于0-1解并得到最優(yōu)結(jié)果。與水平集方法相比，由于中間密度的存在，變密度法的優(yōu)化過程具有更高的連續(xù)性，對(duì)初始設(shè)計(jì)的依賴性更小，也更容易得到最優(yōu)解。改進(jìn)的參數(shù)化水平集方法[12]克服了傳統(tǒng)水平集方法不能自動(dòng)開孔的缺陷，但在優(yōu)化過程中材料密度根據(jù)水平集曲面與零水平集的切割關(guān)系嚴(yán)格服從0-1分布，容易造成結(jié)構(gòu)邊界附近的拓?fù)渫蛔?，降低演化過程的連續(xù)性，導(dǎo)致優(yōu)化結(jié)果難以收斂到最優(yōu)解。當(dāng)改變初始設(shè)計(jì)時(shí)(孔洞數(shù)量和分布)，常常得到不一樣的優(yōu)化結(jié)果，即優(yōu)化結(jié)果對(duì)初始設(shè)計(jì)有依賴性。

為提高結(jié)構(gòu)拓?fù)溲莼倪B續(xù)性，Wei等[18]提出了水平集帶方法，即在零水平集附近引入一個(gè)水平集帶區(qū)域，區(qū)域內(nèi)的材料密度通過水平集函數(shù)進(jìn)行插值，區(qū)域外的材料密度則采用0-1分布。在優(yōu)化過程中通過減小水平集帶寬度逐漸去除區(qū)域內(nèi)的中間材料密度，使目標(biāo)函數(shù)與約束函數(shù)在優(yōu)化結(jié)果收斂到0-1解前保持連續(xù)性，從而提高優(yōu)化過程的穩(wěn)定性，使優(yōu)化結(jié)果更容易趨近于最優(yōu)解。本文工作在文獻(xiàn)[18]基礎(chǔ)上，針對(duì)水平集方法中常見的初始設(shè)計(jì)依賴性問題進(jìn)行了研究，結(jié)果表明水平集帶的引入可以顯著減輕水平集方法的初始設(shè)計(jì)依賴性。隨后探討了帶寬減小速度對(duì)最優(yōu)拓?fù)涞挠绊?，最后通過算例驗(yàn)證了該方法在采用不規(guī)則設(shè)計(jì)域和非結(jié)構(gòu)化網(wǎng)格優(yōu)化問題中的有效性。

1 引入水平集帶的參數(shù)化水平集方法

本節(jié)將首先簡(jiǎn)單介紹參數(shù)化水平集方法，隨后詳細(xì)說明在參數(shù)化水平集方法中引入水平集帶的拓?fù)鋬?yōu)化方法。

1.1 參數(shù)化水平集方法

在水平集方法中，通過引入高一維的水平集函數(shù)的零等值面描述實(shí)體材料區(qū)域邊界，即結(jié)構(gòu)邊界。采用水平集方法對(duì)實(shí)體材料區(qū)域內(nèi)部、外部和邊界的定義可以表示為：

(1)

式中，x∈D?{(x,y)|x,y∈}為設(shè)計(jì)域D中的任意一點(diǎn)，t是時(shí)間，φ(x,t)為水平集函數(shù)，Ω和?Ω分別表示實(shí)體材料的內(nèi)部區(qū)域和邊界。

在鄉(xiāng)村美食方面，里昂的聯(lián)合檢查組與旅游辦公室合作認(rèn)證的里昂傳統(tǒng)餐廳的代名詞“Bouchon Lyonnais”，可以證明餐廳不僅具有質(zhì)量保證，還能烹制最典型和最正宗的里昂菜，對(duì)推動(dòng)鄉(xiāng)村飲食的規(guī)范化經(jīng)營(yíng)起到了重要作用。

在傳統(tǒng)水平集方法中，需要通過求解以下Hamilton-Jacobi偏微分方程更新水平集函數(shù)：

(2)

(3)

式中，gk(x)為第k個(gè)徑向基函數(shù)，αk(t)為相應(yīng)的展開系數(shù)。將式(3)代入式(2)中，則更新水平集函數(shù)的Hamilton-Jacobi偏微分方程可表示為更新展開系數(shù)的常微分方程：

α(ti+1)=α(ti)+G-1ΔtVnδ(φ)

(4)

1.2 在參數(shù)化水平集方法中引入水平集帶

由于水平集方法側(cè)重于邊界演化，有限元分析模型缺少單元中間密度，基于水平集方法的拓?fù)鋬?yōu)化過程容易因?yàn)橥負(fù)渥兓鼓繕?biāo)函數(shù)與約束函數(shù)發(fā)生突變，導(dǎo)致優(yōu)化過程不穩(wěn)定，引起嚴(yán)重的局部最優(yōu)問題。為解決此問題，Wei等[18]提出了水平集帶方法。在傳統(tǒng)水平集方法的零水平集附近引入一個(gè)包含一定寬度的帶狀區(qū)域，在這個(gè)區(qū)域內(nèi)，材料參數(shù)由水平集函數(shù)值的插值函數(shù)確定。如圖1所示，圖1(a)描述了水平集函數(shù)曲面和水平集帶，圖1(b)顯示了水平集帶上下邊界與水平集函數(shù)曲面的交線及材料密度分布。上下邊界之間為水平集帶區(qū)域，其中材料密度的分布與密度法類似，為0-1之間連續(xù)變化的數(shù)值，由插值函數(shù)與水平集函數(shù)值決定。水平集帶的帶寬是上邊界和下邊界之間的寬度，表達(dá)式如下：

圖1 水平集帶方法示意圖Fig.1 Schematic of the level set band method

φb=|φu-φl|

(5)

式中，φu、φl分別代表水平集帶的上邊界和下邊界。在引入水平集帶的水平集方法中，單元密度通過密度插值函數(shù)得到，第i個(gè)單元的密度值為：

ρi=H(φi)

(6)

文中采用光滑的Heaviside函數(shù)，表達(dá)式如下：

(7)

式中,h=φu=-φl，α和ε是趨近于0的正常數(shù)，本文算例均取0.001。

圖2是密度插值方案的示意圖。圖2(a)表示當(dāng)水平集帶的帶寬較大時(shí)，整個(gè)水平集函數(shù)曲面均處于水平集帶的上下邊界之內(nèi)，則單元密度均處于ε和1之間；圖2(b)表示當(dāng)水平集帶的帶寬減小時(shí)，處于上邊界以上的水平集曲面對(duì)應(yīng)的單元密度等于1，處于下邊界以下的水平集曲面對(duì)應(yīng)的單元密度等于ε，上下邊界之間的水平集曲面對(duì)應(yīng)的單元密度處于ε和1之間。隨著優(yōu)化過程的進(jìn)行，可以通過減小帶寬，使灰度區(qū)域逐漸減??；圖2(c)表示當(dāng)水平集帶的帶寬下降為零時(shí)，即上下邊界重合，整個(gè)結(jié)構(gòu)的單元密度接近0-1分布，優(yōu)化結(jié)果呈現(xiàn)出清晰、光滑的設(shè)計(jì)邊界。

圖2 基于水平集帶方法的密度插值方法Fig.2 Density interpolation scheme based on the level set band method

2 拓?fù)鋬?yōu)化模型

本文中僅考慮柔順度最小化問題，目標(biāo)函數(shù)和約束函數(shù)分別表示如下：

(8)

(9)

(10)

(11)

式中，b表示體力。驅(qū)動(dòng)邊界運(yùn)動(dòng)的法向速度場(chǎng)Vn由邊界擴(kuò)展到整個(gè)設(shè)計(jì)域，這由基于形狀導(dǎo)數(shù)的靈敏度分析推導(dǎo)得到[4-5]，其表達(dá)式如下：

Vn=-ε(u):C:ε(u)-ε(u):C:ε(v)+

b·v+(·v)·n+κ(·v)-

(12)

式中，κ=·n為邊界曲率，為控制體積約束的拉格朗日乘子，可以由拉格朗日乘子法或二分法計(jì)算得到。在本文的柔順度最小化問題中，不考慮體力和表面分布力，所以法向速度場(chǎng)可以簡(jiǎn)化為下式：

(13)

通過自然速度擴(kuò)展，可將上式計(jì)算得到的邊界法向速度場(chǎng)擴(kuò)展到整個(gè)設(shè)計(jì)域，代入式(4)中求出更新后的α值，再采用式(3)計(jì)算出水平集函數(shù)值。迭代過程中通過引入并逐漸改變水平集帶帶寬φb實(shí)現(xiàn)水平集帶方法，得到不同區(qū)域材料的本構(gòu)關(guān)系：

C=ρ(φ)C0

(14)

式中，C0是材料的原始本構(gòu)關(guān)系。在每次更新帶寬之后，需要通過式(7)重新計(jì)算結(jié)構(gòu)單元密度分布，直至帶寬減小為0，最終得到具有清晰、光滑邊界的優(yōu)化結(jié)果。優(yōu)化過程的流程如圖3所示。

圖3 優(yōu)化過程的流程圖Fig.3 Flow chart of the optimization process

3 數(shù)值算例

本節(jié)采用數(shù)值算例驗(yàn)證引入水平集帶的參數(shù)化水平集方法可以減輕初始設(shè)計(jì)對(duì)優(yōu)化結(jié)果的影響，并將該方法拓展到不規(guī)則設(shè)計(jì)域優(yōu)化問題。幾個(gè)重要的參數(shù)如下：實(shí)體區(qū)域材料的楊氏模量E=1，孔洞的楊氏模量E=10-9，泊松比υ=0.3,當(dāng)最終的目標(biāo)函數(shù)值與前9個(gè)迭代步的目標(biāo)函數(shù)值均相差小于0.1%時(shí)達(dá)到收斂，優(yōu)化過程終止。算例中采用的徑向基函數(shù)的類型為MultiQuadric(MQ)樣條，其位置分布與有限元網(wǎng)格的節(jié)點(diǎn)重合。

3.1 算例1

如圖4所示的簡(jiǎn)支梁設(shè)計(jì)域的長(zhǎng)和寬分別是80和20，有限元網(wǎng)格劃分為80×20，單元類型為四節(jié)點(diǎn)雙線性正方形單元，體積分?jǐn)?shù)設(shè)置為50%。徑向基函數(shù)的節(jié)點(diǎn)數(shù)量為81×21，與有限元網(wǎng)格的節(jié)點(diǎn)重合。邊界條件為兩端下方簡(jiǎn)支，并在上端中點(diǎn)作用一個(gè)F=1N的集中力。本算例通過對(duì)5種不同的初始設(shè)計(jì)得到的優(yōu)化結(jié)果進(jìn)行比較，研究引入水平集帶方法對(duì)初始設(shè)計(jì)依賴性的改進(jìn)效果。5種不同的初始設(shè)計(jì)如表1所示，其中第5種初始設(shè)計(jì)表示隨機(jī)分布的水平集函數(shù)初始值。水平集帶的初始帶寬為6，每一個(gè)迭代步減小0.1，最終帶寬是0。

表1 不同初始設(shè)計(jì)對(duì)應(yīng)的優(yōu)化結(jié)果Table 1 Optimization results corresponding to different initial designs

圖4 簡(jiǎn)支梁的設(shè)計(jì)域和邊界條件Fig.4 Design domain and boundary condition of simply supported beam

表2給出了第1種初始設(shè)計(jì)對(duì)應(yīng)的優(yōu)化迭代過程。第2、3、4列分別對(duì)應(yīng)零水平集、材料密度分布和水平集曲面。水平集函數(shù)值的初始范圍是[-3,3]，初始水平集曲面全部處于水平集帶的上下邊界之內(nèi)，因此初始的密度分布圖是充滿灰度的。優(yōu)化過程中，水平集函數(shù)值及水平集帶的帶寬發(fā)生改變，部分水平集曲面逐漸處于水平集帶的上下邊界之外，這部分水平集曲面對(duì)應(yīng)的單元密度值等于1或ε，密度分布逐漸趨近于黑白分明。當(dāng)?shù)?0步時(shí)，水平集帶的帶寬減小為0，當(dāng)?shù)?7步時(shí)，達(dá)到收斂條件，優(yōu)化過程結(jié)束。

表2 第1種初始設(shè)計(jì)對(duì)應(yīng)的優(yōu)化過程Table 2 Optimization process corresponding to the first initial design

表1展示了不同的初始設(shè)計(jì)對(duì)應(yīng)的優(yōu)化結(jié)果，從左到右分別為初始設(shè)計(jì)、引入水平集帶方法的優(yōu)化結(jié)果密度分布、引入水平集帶方法的優(yōu)化結(jié)果零水平集、未引入水平集帶方法的優(yōu)化結(jié)果零水平集。由第2、3列可知，水平集帶寬度為0時(shí)，零水平集和密度分布一致。從第3列可知，即使初始設(shè)計(jì)不同，引入水平集帶后也可以得到形狀和結(jié)構(gòu)形式相似的優(yōu)化結(jié)果，而從第4列可知，未引入水平集帶時(shí)，得到的優(yōu)化結(jié)果相差較大，因此可以看出引入水平集帶方法可以有效地減輕初始設(shè)計(jì)對(duì)優(yōu)化結(jié)果的影響。表3給出了引入和未引入水平集帶方法時(shí)優(yōu)化結(jié)果的目標(biāo)函數(shù)值。由表3可知，引入水平集帶方法時(shí)，不同初始設(shè)計(jì)對(duì)應(yīng)的目標(biāo)函數(shù)值更加接近，而且大部分情況下小于未引入水平集帶方法時(shí)對(duì)應(yīng)的目標(biāo)函數(shù)值，說明引入水平集帶方法減小了參數(shù)化水平集方法對(duì)初始設(shè)計(jì)的依賴性并容易得到更優(yōu)的優(yōu)化結(jié)果。

表3 引入和未引入水平集帶方法的優(yōu)化目標(biāo)函數(shù)值Table 3 Objective function values with and without the level set band

3.2 算例2

在3.1節(jié)中，算例1的優(yōu)化結(jié)果說明引入水平集帶的參數(shù)化水平集方法可以有效地減輕初始設(shè)計(jì)的影響，本算例將研究不同的水平集帶寬度減小速度對(duì)優(yōu)化結(jié)果的影響。仍然采用算例1中的簡(jiǎn)支梁優(yōu)化問題，邊界條件、荷載、單元類型和數(shù)量與算例1相同，采用表2中的前四種初始設(shè)計(jì)，在每一種初始設(shè)計(jì)中設(shè)置不同的水平集帶寬度減小速度，初始帶寬均為6，考慮3種迭代步長(zhǎng)每步分別減小0.1、0.2、0.4，最終帶寬均為0。

表4展示了不同初始設(shè)計(jì)在不同的水平集帶寬度減小速度下得到的最優(yōu)拓?fù)洹？梢钥闯?，除了虛線框內(nèi)的兩個(gè)最優(yōu)拓?fù)洳灰恢?，其他結(jié)果均具有相似的形狀和拓?fù)湫问剑f明最優(yōu)結(jié)果對(duì)水平集帶寬度減小的速度并不敏感。結(jié)合算例1中的結(jié)論可知，水平集帶的存在有效地減輕了初始設(shè)計(jì)對(duì)優(yōu)化結(jié)果的影響。從表5可知，當(dāng)水平集帶寬度減小速度不一致時(shí)，不同初始設(shè)計(jì)對(duì)應(yīng)的目標(biāo)函數(shù)值變化幅度很小，大部分結(jié)果小于未引入水平集帶時(shí)的目標(biāo)函數(shù)值，進(jìn)一步論證了優(yōu)化結(jié)果的合理性。

表4 不同的水平集帶寬度減小速度對(duì)應(yīng)的最優(yōu)拓?fù)銽able 4 Optimal topologies corresponding to the different speed reduction of the level set band

表5 不同水平集帶帶寬減小速度對(duì)應(yīng)的目標(biāo)函數(shù)值Table 5 Objective function values corresponding to the different speed reduction of the level set band

3.3 彎曲懸臂梁

參數(shù)化水平集方法的優(yōu)勢(shì)之一是可以方便地求解不規(guī)則設(shè)計(jì)域的優(yōu)化問題，因此，本節(jié)將把水平集帶方法應(yīng)用到采用不規(guī)則設(shè)計(jì)域的優(yōu)化問題中。圖5給出了彎曲懸臂梁的設(shè)計(jì)域和邊界條件，彎曲懸臂梁的左端固定，在右下角作用一個(gè)F=1N 的集中力。由于設(shè)計(jì)域是不規(guī)則的，所以采用線性三角形單元，網(wǎng)格數(shù)量為5 970，節(jié)點(diǎn)數(shù)量為3 113，體積分?jǐn)?shù)為55%。徑向基函數(shù)的節(jié)點(diǎn)數(shù)量為3 113，與有限元網(wǎng)格的節(jié)點(diǎn)重合。水平集函數(shù)初始值為1，水平集帶的初始帶寬為6，每迭代步減小0.2，最終帶寬是0。

圖5 彎曲懸臂梁的設(shè)計(jì)域和邊界條件Fig.5 Design domain and boundary condition of the serpentine beam

圖6是彎曲懸臂梁算例的非結(jié)構(gòu)化有限元網(wǎng)格和拓?fù)鋬?yōu)化結(jié)果。邊界光滑的優(yōu)化結(jié)果類似于由若干粗細(xì)不一的桿件組成的桁架結(jié)構(gòu)。圖7是彎曲懸臂梁算例的目標(biāo)函數(shù)和體積分?jǐn)?shù)的迭代曲線，展示了良好的收斂性能。該結(jié)果與文獻(xiàn)[20]和文獻(xiàn)[21]中的結(jié)果相似，表明本文的優(yōu)化結(jié)果是合理的，驗(yàn)證了水平集帶方法可以直接拓展到非結(jié)構(gòu)化網(wǎng)格中用于求解實(shí)際工程問題。

圖6 彎曲懸臂梁的最優(yōu)拓?fù)銯ig.6 Optimal topology of the serpentine beam

圖7 彎曲懸臂梁優(yōu)化問題的迭代曲線Fig.7 Iterative curves of the optimization problem for the serpentine beam

4 結(jié)語

文中分析了水平集方法中拓?fù)渥兓牟贿B續(xù)性問題，并提出通過引入水平集帶方法可以有效地減輕初始設(shè)計(jì)對(duì)優(yōu)化結(jié)果的影響，并能一定程度上得到更優(yōu)的目標(biāo)函數(shù)值。數(shù)值算例驗(yàn)證了水平集帶方法可以有效地減輕水平集方法對(duì)初始設(shè)計(jì)的依賴性,并顯示優(yōu)化結(jié)果對(duì)水平集帶寬度的減小速度不敏感。最后采用彎曲懸臂梁算例說明該方法可以直接拓展到不規(guī)則設(shè)計(jì)域優(yōu)化問題中用于求解復(fù)雜的實(shí)際工程問題，并展示出良好的收斂性能。

文中的數(shù)值算例沒有考慮表面力的作用，由于在水平集帶的帶寬減小到0之前，結(jié)構(gòu)沒有清晰的邊界，因此表面力的存在會(huì)對(duì)帶寬的大小有一定的限制。此外，本文方法還可用于處理帶有中間密度的材料如孔隙材料的優(yōu)化問題，這些也是后續(xù)將要開展的研究工作。