胡 清 元， 李 小 漁， 梁 緣

(1.江南大學(xué) 理學(xué)院，江蘇 無(wú)錫 214122；2.大連理工大學(xué) 運(yùn)載工程與力學(xué)學(xué)部，遼寧 大連 116024 )

0 引 言

結(jié)構(gòu)拓?fù)鋬?yōu)化旨在根據(jù)給定的約束條件，尋求最優(yōu)的材料分布，發(fā)揮材料最高性價(jià)比[1]．基于傳統(tǒng)有限元分析的拓?fù)鋬?yōu)化方法發(fā)展成熟，涌現(xiàn)出了眾多的優(yōu)秀方法，例如基于單元懲罰密度的SIMP方法[2-3]、基于演化結(jié)構(gòu)的ESO方法[4]、基于雙向漸進(jìn)結(jié)構(gòu)優(yōu)化的BESO方法[5]、基于水平集的level-set優(yōu)化方法[6-7]，以及基于顯示邊界描述的MMC/MMV方法[8-9]等．這些方法仍在不斷發(fā)展，也被用于指導(dǎo)實(shí)際工程問(wèn)題[10]，取得了不同程度的成效．

然而，就目前廣泛被采用的基于有限元SIMP方法的拓?fù)鋬?yōu)化來(lái)說(shuō)，仍存在兩個(gè)亟待解決的問(wèn)題．首先，受限于傳統(tǒng)有限元單元的以直代曲特性，拓?fù)鋬?yōu)化所得結(jié)果邊界往往不夠光滑，需要網(wǎng)格足夠細(xì)密或經(jīng)可視化后處理才能得到光滑邊界；其次，SIMP方法通常會(huì)產(chǎn)生灰度單元，導(dǎo)致該處材料布局結(jié)果模糊．這些問(wèn)題導(dǎo)致優(yōu)化結(jié)果難以指導(dǎo)后續(xù)設(shè)計(jì)加工．

等幾何分析于2005年被Hughes等提出[11]，該方法直接采用CAD模型中天然內(nèi)含的曲邊單元進(jìn)行有限元分析，將設(shè)計(jì)和分析融為一體．得益于等幾何分析的優(yōu)勢(shì)，Gao等[12]采用密度分布函數(shù)定義結(jié)構(gòu)拓?fù)?，以得到光滑性和連續(xù)性可控的優(yōu)化結(jié)果；董振宇[13]基于Catmull-Clark體細(xì)分方法，開(kāi)展了多分辨率等幾何拓?fù)鋬?yōu)化的一系列研究；楊佳明等[14]采用混合B樣條進(jìn)行材料密度分布和等幾何分析，研究了實(shí)體模型的等幾何優(yōu)化．

本文基于等幾何分析，采用離散變量拓?fù)鋬?yōu)化方法，針對(duì)平面問(wèn)題展開(kāi)理論算法研究．基于序列近似整數(shù)規(guī)劃和正則松弛算法的離散變量拓?fù)鋬?yōu)化算法由Liang等[15]提出，該方法將單元密度嚴(yán)格限制為0或1，因此可以得到完全黑白的材料布局．將離散變量拓?fù)鋬?yōu)化方法引入等幾何分析框架，一方面可以利用等幾何分析曲邊單元的優(yōu)勢(shì)，另一方面可以優(yōu)化出清晰的結(jié)構(gòu)邊界，使得優(yōu)化結(jié)果對(duì)后續(xù)設(shè)計(jì)加工更具指導(dǎo)性．

1 等幾何分析簡(jiǎn)介

等幾何分析采用非均勻有理B樣條(NURBS)作為基函數(shù)同時(shí)插值幾何和位移場(chǎng)，其本質(zhì)可理解為采用NURBS作為形函數(shù)的等參有限元．因此，本章重點(diǎn)介紹NURBS及由此帶來(lái)的等幾何分析特性．

從一維單變量B樣條函數(shù)出發(fā)，B樣條函數(shù)可由預(yù)先給定的節(jié)點(diǎn)矢量Ξ生成：

Ξ=(ξ1ξ2…ξn+p+1)

(1)

節(jié)點(diǎn)下標(biāo)i=1，2，…，n+p+1，最后一個(gè)下標(biāo)編號(hào)給出了所生成的B樣條函數(shù)的個(gè)數(shù)n和階次p．接下來(lái)，利用Cox-de Boor遞推公式

(2)

生成B樣條函數(shù)Ni，p，如圖1(a)所示．

由于B樣條函數(shù)無(wú)法描述圓錐曲線，故對(duì)基函數(shù)進(jìn)行有理化處理，給定權(quán)重wi，即可得到NURBS函數(shù)：

(3)

當(dāng)權(quán)重都為1時(shí)，NURBS函數(shù)退化為B樣條函數(shù)．搭配相應(yīng)控制點(diǎn)pi，可生成曲線

(4)

通常，CAD和等幾何分析所采用的節(jié)點(diǎn)矢量以0開(kāi)始、以1結(jié)束，并要求所生成的NURBS函數(shù)是開(kāi)的，即對(duì)于p次NURBS函數(shù)節(jié)點(diǎn)矢量Ξ，其開(kāi)始的0和結(jié)束的1被重復(fù)p+1次．圖1(b)給出了生成圓弧曲線的示例．

生成二維曲面需要擴(kuò)充維度，為此，給定另一組節(jié)點(diǎn)矢量Ψ和權(quán)重，得到NURBS函數(shù)Rj，q(η)，其中η為該維度的方向參數(shù)，j為函數(shù)個(gè)數(shù)，q表示函數(shù)階次．搭配相應(yīng)控制點(diǎn)pi，j，得到張量積曲面

(5)

與幾何模型的插值類似，等幾何分析采用相同的NURBS基函數(shù)近似位移場(chǎng)：

(6)

(a)NURBS基函數(shù)

其中ui，j=(ui，jvi，j)，為對(duì)應(yīng)控制點(diǎn)pi，j上的自由度，ui,j和vi,j分別表示該點(diǎn)沿x軸和y軸的位移．

等幾何分析后續(xù)步驟與等參有限元相同．需要注意的是，與傳統(tǒng)的等參有限元相比，等幾何分析直接采用原始的幾何模型進(jìn)行分析，無(wú)須劃分網(wǎng)格，也不產(chǎn)生離散誤差．但實(shí)際上等幾何分析仍然依賴網(wǎng)格，該網(wǎng)格在幾何建模的過(guò)程中由節(jié)點(diǎn)矢量Ξ×Ψ自然地生成．

2 離散變量拓?fù)鋬?yōu)化方法

考慮在給定材料用量的約束條件下，求結(jié)構(gòu)最小柔順性的問(wèn)題．將結(jié)構(gòu)離散為M個(gè)單元，令ρ為單元密度，問(wèn)題的優(yōu)化列式為

(7)

為求得目標(biāo)函數(shù)c(ρ)的最小值，需要求出c關(guān)于ρi的變化率，即目標(biāo)函數(shù)靈敏度．借鑒SIMP方法，假設(shè)單元彈性模量與單元密度之間的關(guān)系為

E=E0ρP

(8)

其中P=3是懲罰參數(shù)．目標(biāo)函數(shù)的靈敏度為

(9)

根據(jù)Liang等[15]的方法，將式(7)轉(zhuǎn)化為求解如下的整數(shù)線性規(guī)劃問(wèn)題：

(10)

由于整數(shù)線性規(guī)劃的約束條件以及目標(biāo)函數(shù)的組合復(fù)雜性，很難獲得全局最優(yōu)解．通過(guò)轉(zhuǎn)化整數(shù)松弛變量約束以及引入拉格朗日乘子λ和σ，得到如下的對(duì)偶關(guān)系：

(11)

同時(shí)有

(12)

以及

(13)

這里β=2×104，是攝動(dòng)參數(shù)．給定λ初始值，通過(guò)求解式(13)得到對(duì)偶變量σi，結(jié)合對(duì)偶關(guān)系式(11)得到原設(shè)計(jì)變量ρi，如果更新后的設(shè)計(jì)變量ρi滿足收斂準(zhǔn)則，即認(rèn)為達(dá)到了近似解，跳出迭代；否則通過(guò)式(12)更新λ，繼續(xù)迭代．

3 算法實(shí)現(xiàn)細(xì)節(jié)

本研究所述算法的代碼實(shí)現(xiàn)基于開(kāi)源等幾何分析框架IGAFEM[16]和離散變量拓?fù)鋬?yōu)化開(kāi)源代碼DVTOPCRA[17]．在將兩者有機(jī)融合的基礎(chǔ)上，需要注意以下細(xì)節(jié)：

首先，需要注意的是，Liang等提出的離散變量?jī)?yōu)化算法[15]是基于密度的．考慮到如果將等幾何分析的控制點(diǎn)密度作為設(shè)計(jì)變量，那么即使單個(gè)控制點(diǎn)的密度是非黑即白的，通過(guò)NURBS函數(shù)插值后所得曲線意義也將變得不明確，所涉及單元可能仍是灰度的．因此，本研究將單元密度作為設(shè)計(jì)變量．在參數(shù)空間Ξ×Ψ得到0或1的密度后，通過(guò)式(5)映射到物理空間，即得到包含曲邊單元的材料布局，如圖2所示．

圖2 參數(shù)空間單元映射到物理空間曲邊單元Fig.2 Mapping from elements in parameter space to curved elements in physical space

此外，為避免出現(xiàn)棋盤格現(xiàn)象，需對(duì)靈敏度進(jìn)行過(guò)濾．傳統(tǒng)有限元中的靈敏度過(guò)濾方法高度依賴結(jié)構(gòu)化網(wǎng)格，需要對(duì)鄰近單元進(jìn)行搜索．由于等幾何分析中存在曲邊單元，本文采用基于Helmh?ltz 方程的隱式過(guò)濾方法[18]，該方法無(wú)須搜索鄰近單元，易于處理復(fù)雜網(wǎng)格．過(guò)濾公式為

(14)

(15)

為標(biāo)量問(wèn)題的等幾何剛度陣，Re為基函數(shù)插值矩陣，?Re為基函數(shù)梯度矩陣，

(16)

(17)

最終得到過(guò)濾后的單元靈敏度列陣

(18)

由于優(yōu)化算法求解的是近似的子問(wèn)題，需要嚴(yán)格控制迭代過(guò)程，保證近似子問(wèn)題的解緩慢收斂到真實(shí)解．這里采用體分比縮減因子χ，用于逐漸減少材料用量約束，在當(dāng)前體分比約束下計(jì)算收斂后，再按χ縮減體分比約束繼續(xù)計(jì)算．根據(jù)數(shù)值測(cè)試，參數(shù)χ取值范圍為[0.95，1.00)．本文令體分比縮減因子χ=0.98，即假定每一輪迭代結(jié)構(gòu)中只有2%的材料發(fā)生改變，以限制設(shè)計(jì)變量的變化范圍，起到運(yùn)動(dòng)極限的作用，從而保證基于靈敏度的整數(shù)規(guī)劃子問(wèn)題(10)的近似精度．在未來(lái)工作中可以使用Liang等提出的基于信賴域的方法精確施加運(yùn)動(dòng)極限[19]．

4 數(shù)值算例

4.1 懸臂梁

圖3 懸臂梁Fig.3 Cantilever beam

本文方法所得優(yōu)化結(jié)果見(jiàn)圖4，迭代過(guò)程中的結(jié)構(gòu)柔順性變化和體分比v變化見(jiàn)圖5．目標(biāo)函數(shù)值為3 210 N·m．圖6展示了優(yōu)化后結(jié)構(gòu)的應(yīng)力分布，可見(jiàn)應(yīng)力分布均勻，整體處于較低水平．

圖4 懸臂梁優(yōu)化結(jié)果Fig.4 Optimization result of the cantilever beam

圖5 懸臂梁迭代數(shù)據(jù)Fig.5 Iteration data of the cantilever beam

圖6 懸臂梁應(yīng)力分布Fig.6 Stress distribution of the cantilever beam

為了驗(yàn)證本文算法的有效性，在ABAQUS中采用了細(xì)密的有限元網(wǎng)格，圖7中展示了ABAQUS的優(yōu)化結(jié)果，目標(biāo)函數(shù)值為3 394 N·m．由于連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化依賴于網(wǎng)格劃分，而ABAQUS采用的有限元網(wǎng)格更加細(xì)密，因此本文方法所得結(jié)果與ABAQUS所得結(jié)果在構(gòu)形上有所不同．離散變量下拓?fù)鋬?yōu)化問(wèn)題的數(shù)學(xué)本質(zhì)是偏微分方程約束的非線性整數(shù)規(guī)劃，因此不同的離散方法也會(huì)得到不同的局部最優(yōu)解．圖4與圖7兩種構(gòu)形的目標(biāo)函數(shù)值相近，則一定程度上可以說(shuō)明兩種優(yōu)化結(jié)果構(gòu)形的受力特性相似．

圖7 懸臂梁ABAQUS優(yōu)化結(jié)果Fig.7 Optimization result of the cantilever beam by ABAQUS

采用基于傳統(tǒng)有限元的SIMP方法計(jì)算，網(wǎng)格劃分與本文方法(圖4)相同，優(yōu)化結(jié)果見(jiàn)圖8．可見(jiàn)SIMP方法優(yōu)化結(jié)果邊界產(chǎn)生大量灰度單元，邊界模糊，而本文方法(圖4)邊界清晰，體現(xiàn)了本文方法的優(yōu)勢(shì)．

圖8 懸臂梁SIMP法優(yōu)化結(jié)果Fig.8 Optimization result of the cantilever beam by SIMP method

4.2 四分之一圓環(huán)

圖9 四分之一圓環(huán)Fig.9 One-quarter annulus

采用本文方法得到的優(yōu)化結(jié)果見(jiàn)圖10，目標(biāo)函數(shù)和材料用量迭代曲線見(jiàn)圖11，目標(biāo)函數(shù)值為28.6 N·m．應(yīng)力分布圖12表明，優(yōu)化后的構(gòu)形應(yīng)力分布較為平均．圖13為ABAQUS的優(yōu)化結(jié)果，目標(biāo)函數(shù)值為30.5 N·m．本文方法所得優(yōu)化結(jié)果與ABAQUS結(jié)果相似，兩種構(gòu)形目標(biāo)函數(shù)相近，驗(yàn)證了方法的有效性．

圖10 四分之一圓環(huán)優(yōu)化結(jié)果Fig.10 Optimization result of the one-quarter annulus

圖11 四分之一圓環(huán)迭代數(shù)據(jù)Fig.11 Iteration data of the one-quarter annulus

圖12 四分之一圓環(huán)應(yīng)力分布Fig.12 Stress distribution of the one-quarter annulus

圖13 四分之一圓環(huán)ABAQUS優(yōu)化結(jié)果Fig.13 Optimization result of the one-quarter annulus by ABAQUS

此外，ABAQUS分析所用自由度數(shù)為1.9×105，本文所用自由度數(shù)為6.7×104，相比ABAQUS節(jié)省了計(jì)算量．但本文所得到的結(jié)果中存在大量曲邊單元，例如，將圖10中畫框部分放大為圖14(a)，可見(jiàn)在本文方法的優(yōu)化結(jié)果中，對(duì)局部某些包含曲線邊界的區(qū)域，采用較少單元即可使結(jié)構(gòu)邊界清晰、光滑；而ABAQUS單元對(duì)應(yīng)部分(圖14(b))雖然單元更密，但幾乎所有曲線邊界均呈鋸齒狀，通常需要采用黑線描邊的方式進(jìn)行可視化后處理．

(a)本文

5 結(jié) 語(yǔ)

針對(duì)基于傳統(tǒng)有限元SIMP拓?fù)鋬?yōu)化方法存在的問(wèn)題，本文在等幾何分析框架內(nèi)，結(jié)合離散變量拓?fù)鋬?yōu)化算法展開(kāi)了研究．本文提出的優(yōu)化方法結(jié)合了等幾何分析和離散變量拓?fù)鋬?yōu)化的優(yōu)勢(shì)，一方面允許曲邊單元存在，另一方面單元密度為非黑即白的，因此優(yōu)化結(jié)果在局部某些曲線邊界上清晰光滑．本研究還采用了基于Helmh?ltz方程的隱式過(guò)濾方法，節(jié)省了搜索鄰近單元的計(jì)算量，對(duì)曲邊單元過(guò)濾效果良好．

通過(guò)算例對(duì)比，驗(yàn)證了方法的有效性，展示了方法的優(yōu)勢(shì)．與傳統(tǒng)有限元相比，本文方法優(yōu)化計(jì)算所得的構(gòu)形在整體上更加清晰、在局部上更加光滑，對(duì)后續(xù)設(shè)計(jì)和加工更具指導(dǎo)意義．