基于改進水平集的六邊形網格拓撲優化

張健林,李海艷,周健松,陳慶杰

(廣東工業大學機電工程學院,廣州 510006)

0 引言

拓撲優化是在給定約束條件下,確定設計域內材料的最優分配[1]。水平集方法是拓撲優化的一種主流方法,它通過將設計域嵌入高一維的水平集函數的零級界面,通過水平集函數的整體變化,引起設計域的拓撲變化,從而達到拓撲優化的目的[2]。水平集方法具有最終拓撲結構清晰,邊界光滑的特點,因此廣泛應用于工程實際中[3-5]。

傳統水平集的演化是在界面的水平方向進行,需要求解偏微分方程[6],且用到了基于矩形網格的有限差分法,上風方案[7],在處理多邊形網格上存在著困難。LIU等[8]提出了一種豎直方向演化的水平集方法,以基于節點的速度驅動水平集函數演化,可以很方便地應用于不同網格。

KUMAR[9]提出了一種基于六邊形網格的拓撲優化方法,通過將有限元中的矩形網格換為六邊形網格,可以避免單元之間的非奇異連接,將點連接轉變為線連接,緩解棋盤格現象。該種方法是基于變密度法[10]實現的,自然而然地繼承了變密度法的一些數值缺陷,在最終拓撲結構中出現了大面積的中間密度[11],也即是灰度。

對于基于六邊形網格的拓撲優化方法,目前只有用變密度法實現,還沒有發現有使用水平集方法實現,針對文獻[9]方法最終拓撲出現大面積灰度的問題,引入LIU等[8]改進的水平集方法,結合六邊形網格,提出一種新的方法,解決受體積約束的柔度最小化問題,并用該方法與文獻[9]方法進行實驗對比,驗證所提方法的有效性。

1 水平集函數表達

水平集函數φ的定義如下:

(1)

式中:X為設計域內的點,D為設計域,Ω為材料域,?Ω為材料域的邊界。

執行有限元分析之前,首先需要將水平集函數轉換為物理模型,使用精確的Heaviside函數將水平集函數映射到結構模型上。

(2)

基于H(φ)進行等效的密度變換。

(3)

式中:ρe為材料域內的單元密度。

從式(3)得到了單元密度之后,那么楊氏模量Ee就可以通過以下式子計算:

Ee(ρe)=Emin+(ρe)p(E0-Emin)

(4)

式中:Emin是用于防止剛度奇異的微小正數,p為懲罰系數,E0為實體元素的剛度。

根據文獻[8],水平集的演化最簡單的形式為:

(5)

式中:t為引入的虛擬時間步長,水平集的演化依賴于節點的靈敏度,節點靈敏度VNS的定義如下,它為包含該節點的元素靈敏度之和的平均值。

(6)

式中:ne為有著公共節點的元素總數,R為要優化的目標函數。

2 拓撲優化問題求解

受體積約束的柔度最小化問題表示為:

(7)

式中:C為結構的柔度,ρe為單元密度,Ne為單元總數,ue為單元位移,ke為單元剛度矩陣,K為全局剛度,U為全局位移,F為外載荷,ve為單元體積,Vf為目標體積分數。

拉格朗日函數L可表達為:

(8)

式中:Λ為朗格朗日乘子。

拉格朗日函數對于密度進行微分,可得:

(9)

式中:ve可以看作為1,由KKT條件,可得:

(10)

從式(6)得到,節點靈敏度是由元素靈敏度擴展而來,故式(10)可改寫為:

VNS+Λ=0

(11)

根據文獻[8],可得修正后的水平集演化方程:

(12)

式中:拉格朗日乘子Λ是通過二分法尋找到的。

柔度C對于單元密度ρe微分,便得到了元素靈敏度:

(13)

3 實現細節

3.1 六邊形網格

如圖1所示,黑色字體為單元編號,灰色字體為節點編號,以3×3網格規模為示例,從左下角開始,采用從左到右,從下到上的規律進行編號。

圖1 六邊形網格編號

對于六邊形網格,節點有3種情況,用于式(6)節點靈敏度的計算。

(1)3個單元共用一個節點。

(2)2個單元共用一個節點。

(3)不與其他單元共用的節點。

3.2 過濾器

在本文方法中,靈敏度濾波是唯一的濾波方案,為了實驗對比,故對于文獻[9]也采用相同的靈敏度濾波方案進行處理。

元素j的過濾靈敏度是通過在元素j的鄰域上使用濾波半徑rmin定義的元素靈敏度的加權平均得到的,過濾后的靈敏度場為:

(14)

式中:Nm是元素m的集合,其中元素m的中心到元素j中心距離dist(ρj,ρm)小于濾波半徑rmin,Hjm是一個權重因子,由下式給出:

Hjm=max(0,rmin-dist(ρj,ρm))

(15)

權重的定義使得靠近元素j的元素與遠離元素j的元素相比,對元素j的過濾靈敏度的貢獻更大。

3.3 灰度單元

由于式(3)的緣故,水平集函數切割邊界的單元,所提方法在且只在結構邊界處存在灰度單元,大大減少了灰度單元的數量,隨著網格的細分,灰度單元在最終拓撲結構幾乎可以忽略不計。

3.4 收斂準則

當優化迭代過程中柔度變化量不超過一個微小的閾值ε,或是到達定義的最大迭代次數itermax,則認為迭代已經收斂。

(16)

式中:ci為第i次迭代的柔度,ci-1為第i-1次迭代的柔度。本文柔度變化量的閾值ε取為1×10-4,最大迭代次數取為200。

4 實驗算例

采用半MBB梁與米歇爾結構兩組算例來驗證本文所提方法的有效性,在這兩個算例中,采用一樣的參數數值,E0=1,Emin=10-8,p=3。

4.1 半MBB梁

如圖2所示,半MBB梁設計域的長寬比為3∶1,目標體積選為0.5,單位力F施加在設計域的左上角,對比了文獻[9]與本文方法在網格分別取為75×25,165×55,225×75時的最終拓撲結構及柔度,最終拓撲結構如表1所示,柔度分析如表2所示。

表1 半MBB梁拓撲優化

表2 半MBB梁柔度分析

圖2 半MBB梁

表1對比了文獻[9]方法與本文方法得到的最終拓撲結構,可以看出本文方法得到了與文獻[9]相似的拓撲構型,且比文獻[9]方法更加清晰。

表2對比分析了文獻[9]方法與本文方法最終獲得的半MBB梁柔度,從結果上來看,本文方法獲得了比文獻[9]方法更低的柔度,且隨著網格的細分,柔度下降率的數值更大。

4.2 米歇爾結構

如圖3所示,米歇爾結構設計域的長寬比為1∶1,目標體積選為0.2,單位力F施加在設計域的左下角,對比了文獻[9]與本文方法在網格分別取為45×45,85×85,125×125時的最終拓撲結構及柔度,最終拓撲結構如表3所示,柔度分析如表4所示。

表3 米歇爾結構拓撲優化

表4 米歇爾結構柔度分析

圖3 米歇爾結構

表3對比了文獻[9]方法與本文方法得到的最終拓撲結構,可以看出本文方法得到了與文獻[9]相似的拓撲構型,且比文獻[9]方法更加清晰。

表4對比分析了文獻[9]方法與本文方法最終獲得的米歇爾結構柔度,從結果上來看,本文方法獲得了比文獻[9]方法更低的柔度,且隨著網格的細分,柔度下降率的數值更大。

5 結論

針對基于變密度法的六邊形網格拓撲優化中,獲得的結構存在大量中間密度的問題,引入改進的水平集方法,進而提出一種新的方法,并通過實驗算例進行驗證,通過分析實驗算例結果,可得出以下結論:

(1)本文所提方法能獲得柔度更小,結構清晰,邊界光滑的拓撲優化結構,為實際工程應用中的生產加工提供了便捷,有效的方法。

(2)隨著網格規模的增大,本文所提方法能進一步降低柔度,抑制灰度單元。

本文將該方法運用到了二維拓撲優化,并取得了較好的結果,未來的工作將致力于三維拓撲優化。