基于等幾何分析的離散變量拓撲優化方法研究

胡 清 元， 李 小 漁， 梁 緣

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

0 引 言

結構拓撲優化旨在根據給定的約束條件，尋求最優的材料分布，發揮材料最高性價比[1]．基于傳統有限元分析的拓撲優化方法發展成熟，涌現出了眾多的優秀方法，例如基于單元懲罰密度的SIMP方法[2-3]、基于演化結構的ESO方法[4]、基于雙向漸進結構優化的BESO方法[5]、基于水平集的level-set優化方法[6-7]，以及基于顯示邊界描述的MMC/MMV方法[8-9]等．這些方法仍在不斷發展，也被用于指導實際工程問題[10]，取得了不同程度的成效．

然而，就目前廣泛被采用的基于有限元SIMP方法的拓撲優化來說，仍存在兩個亟待解決的問題．首先，受限于傳統有限元單元的以直代曲特性，拓撲優化所得結果邊界往往不夠光滑，需要網格足夠細密或經可視化后處理才能得到光滑邊界；其次，SIMP方法通常會產生灰度單元，導致該處材料布局結果模糊．這些問題導致優化結果難以指導后續設計加工．

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

本文基于等幾何分析，采用離散變量拓撲優化方法，針對平面問題展開理論算法研究．基于序列近似整數規劃和正則松弛算法的離散變量拓撲優化算法由Liang等[15]提出，該方法將單元密度嚴格限制為0或1，因此可以得到完全黑白的材料布局．將離散變量拓撲優化方法引入等幾何分析框架，一方面可以利用等幾何分析曲邊單元的優勢，另一方面可以優化出清晰的結構邊界，使得優化結果對后續設計加工更具指導性．

1 等幾何分析簡介

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

從一維單變量B樣條函數出發，B樣條函數可由預先給定的節點矢量Ξ生成：

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

(1)

節點下標i=1，2，…，n+p+1，最后一個下標編號給出了所生成的B樣條函數的個數n和階次p．接下來，利用Cox-de Boor遞推公式

(2)

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

由于B樣條函數無法描述圓錐曲線，故對基函數進行有理化處理，給定權重wi，即可得到NURBS函數：

(3)

當權重都為1時，NURBS函數退化為B樣條函數．搭配相應控制點pi，可生成曲線

(4)

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

生成二維曲面需要擴充維度，為此，給定另一組節點矢量Ψ和權重，得到NURBS函數Rj，q(η)，其中η為該維度的方向參數，j為函數個數，q表示函數階次．搭配相應控制點pi，j，得到張量積曲面

(5)

與幾何模型的插值類似，等幾何分析采用相同的NURBS基函數近似位移場：

(6)

(a)NURBS基函數

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

等幾何分析后續步驟與等參有限元相同．需要注意的是，與傳統的等參有限元相比，等幾何分析直接采用原始的幾何模型進行分析，無須劃分網格，也不產生離散誤差．但實際上等幾何分析仍然依賴網格，該網格在幾何建模的過程中由節點矢量Ξ×Ψ自然地生成．

2 離散變量拓撲優化方法

考慮在給定材料用量的約束條件下，求結構最小柔順性的問題．將結構離散為M個單元，令ρ為單元密度，問題的優化列式為

(7)

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

E=E0ρP

(8)

其中P=3是懲罰參數．目標函數的靈敏度為

(9)

根據Liang等[15]的方法，將式(7)轉化為求解如下的整數線性規劃問題：

(10)

由于整數線性規劃的約束條件以及目標函數的組合復雜性，很難獲得全局最優解．通過轉化整數松弛變量約束以及引入拉格朗日乘子λ和σ，得到如下的對偶關系：

(11)

同時有

(12)

以及

(13)

這里β=2×104，是攝動參數．給定λ初始值，通過求解式(13)得到對偶變量σi，結合對偶關系式(11)得到原設計變量ρi，如果更新后的設計變量ρi滿足收斂準則，即認為達到了近似解，跳出迭代；否則通過式(12)更新λ，繼續迭代．

3 算法實現細節

本研究所述算法的代碼實現基于開源等幾何分析框架IGAFEM[16]和離散變量拓撲優化開源代碼DVTOPCRA[17]．在將兩者有機融合的基礎上，需要注意以下細節：

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

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

此外，為避免出現棋盤格現象，需對靈敏度進行過濾．傳統有限元中的靈敏度過濾方法高度依賴結構化網格，需要對鄰近單元進行搜索．由于等幾何分析中存在曲邊單元，本文采用基于Helmh?ltz 方程的隱式過濾方法[18]，該方法無須搜索鄰近單元，易于處理復雜網格．過濾公式為

(14)

(15)

為標量問題的等幾何剛度陣，Re為基函數插值矩陣，?Re為基函數梯度矩陣，

(16)

(17)

最終得到過濾后的單元靈敏度列陣

(18)

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

4 數值算例

4.1 懸臂梁

圖3 懸臂梁Fig.3 Cantilever beam

本文方法所得優化結果見圖4，迭代過程中的結構柔順性變化和體分比v變化見圖5．目標函數值為3 210 N·m．圖6展示了優化后結構的應力分布，可見應力分布均勻，整體處于較低水平．

圖4 懸臂梁優化結果Fig.4 Optimization result of the cantilever beam

圖5 懸臂梁迭代數據Fig.5 Iteration data of the cantilever beam

圖6 懸臂梁應力分布Fig.6 Stress distribution of the cantilever beam

為了驗證本文算法的有效性，在ABAQUS中采用了細密的有限元網格，圖7中展示了ABAQUS的優化結果，目標函數值為3 394 N·m．由于連續體結構拓撲優化依賴于網格劃分，而ABAQUS采用的有限元網格更加細密，因此本文方法所得結果與ABAQUS所得結果在構形上有所不同．離散變量下拓撲優化問題的數學本質是偏微分方程約束的非線性整數規劃，因此不同的離散方法也會得到不同的局部最優解．圖4與圖7兩種構形的目標函數值相近，則一定程度上可以說明兩種優化結果構形的受力特性相似．

圖7 懸臂梁ABAQUS優化結果Fig.7 Optimization result of the cantilever beam by ABAQUS

采用基于傳統有限元的SIMP方法計算，網格劃分與本文方法(圖4)相同，優化結果見圖8．可見SIMP方法優化結果邊界產生大量灰度單元，邊界模糊，而本文方法(圖4)邊界清晰，體現了本文方法的優勢．

圖8 懸臂梁SIMP法優化結果Fig.8 Optimization result of the cantilever beam by SIMP method

4.2 四分之一圓環

圖9 四分之一圓環Fig.9 One-quarter annulus

采用本文方法得到的優化結果見圖10，目標函數和材料用量迭代曲線見圖11，目標函數值為28.6 N·m．應力分布圖12表明，優化后的構形應力分布較為平均．圖13為ABAQUS的優化結果，目標函數值為30.5 N·m．本文方法所得優化結果與ABAQUS結果相似，兩種構形目標函數相近，驗證了方法的有效性．

圖10 四分之一圓環優化結果Fig.10 Optimization result of the one-quarter annulus

圖11 四分之一圓環迭代數據Fig.11 Iteration data of the one-quarter annulus

圖12 四分之一圓環應力分布Fig.12 Stress distribution of the one-quarter annulus

圖13 四分之一圓環ABAQUS優化結果Fig.13 Optimization result of the one-quarter annulus by ABAQUS

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

(a)本文

5 結 語

針對基于傳統有限元SIMP拓撲優化方法存在的問題，本文在等幾何分析框架內，結合離散變量拓撲優化算法展開了研究．本文提出的優化方法結合了等幾何分析和離散變量拓撲優化的優勢，一方面允許曲邊單元存在，另一方面單元密度為非黑即白的，因此優化結果在局部某些曲線邊界上清晰光滑．本研究還采用了基于Helmh?ltz方程的隱式過濾方法，節省了搜索鄰近單元的計算量，對曲邊單元過濾效果良好．

通過算例對比，驗證了方法的有效性，展示了方法的優勢．與傳統有限元相比，本文方法優化計算所得的構形在整體上更加清晰、在局部上更加光滑，對后續設計和加工更具指導意義．