連續(xù)體結(jié)構(gòu)拓撲優(yōu)化敏度過濾研究*

張國鋒，徐 雷，李大雙，余方超

(四川大學(xué)機械工程學(xué)院，成都 610065)

0 引言

連續(xù)體結(jié)構(gòu)拓撲優(yōu)化[1]已成為當前拓撲優(yōu)化研究方向的熱點問題之一，其在諸多領(lǐng)域具有普遍應(yīng)用。變密度法[2]作為常用的拓撲優(yōu)化方法之一，因其設(shè)計變量少、程序簡單及應(yīng)用范圍廣等特點被廣泛應(yīng)用，其實質(zhì)是單元密度為0～1的離散變量之間的排列組合的問題。因進行有限元計算分析，導(dǎo)致在整個的優(yōu)化過程中產(chǎn)生如棋盤格現(xiàn)象、網(wǎng)格依賴性等數(shù)值不穩(wěn)定的問題[3]，使其直接制造性不強，制約了變密度法在結(jié)構(gòu)拓撲優(yōu)化領(lǐng)域的發(fā)展。

國內(nèi)外學(xué)者針對上述問題進行了廣泛和深入的研究，Sigmund O等[4]提出半徑過濾方法，有效消除棋盤格問題，但優(yōu)化結(jié)果常出現(xiàn)邊界擴散現(xiàn)象。Haber R B等[5]提出周長控制法，需要人為反復(fù)調(diào)試參數(shù)，且對局部棋盤格現(xiàn)象抑制效果不明顯。Petersson J等[6]提出局部密度斜率控制法，由于引入額外的約束條件，導(dǎo)致運算耗時，求解難度大。朱劍峰等[7]提出一種通過影響因子控制敏度修正程度的方法，盡管對迭代速率，結(jié)構(gòu)剛度均有提升，但仍出現(xiàn)一定的邊界擴散現(xiàn)象。陳垂福等[8]通過對中心及周邊元素采用不同權(quán)重來抑制棋盤格現(xiàn)象，但并不能有效控制最小尺寸。杜義賢等[9]借鑒粒子群優(yōu)化算法中粒子狀態(tài)的更新方法來抑制棋盤格現(xiàn)象，但在優(yōu)化過程中可能出現(xiàn)細小孔洞的情況。龍凱等[10]提出一種考慮密度梯度的方法，但可能出現(xiàn)結(jié)構(gòu)孔洞數(shù)量發(fā)生變化，影響最后的優(yōu)化結(jié)果。

本文在現(xiàn)有研究的基礎(chǔ)上提出了一種改進的敏度過濾方法，該方法通過引入新的卷積因子，結(jié)合原有卷積因子建立三者在一定權(quán)重比分配下的數(shù)值關(guān)系，并采用一種帶有預(yù)設(shè)修正權(quán)值的方法，來弱化邊界擴散的問題。該方法在有效消除棋盤格等問題的同時，有效提升迭代速度，降低柔度收斂值，提升結(jié)構(gòu)剛度。以柔度最小化為優(yōu)化目標，通過選取多個二維平面結(jié)構(gòu)算例驗證了本文方法的可行性。

1 變密度法拓撲優(yōu)化理論及數(shù)學(xué)建模

固體各向同性材料的懲罰法[2](Solid Isotropic Microstructures with Penalization，SIMP)即變密度法是將設(shè)計域離散成一個由X軸及Y軸上的元素集合定義的有限元單元，選擇材料密度作為結(jié)構(gòu)拓撲優(yōu)化的設(shè)計變量，根據(jù)最優(yōu)性準則或數(shù)學(xué)規(guī)劃方法，在參考域內(nèi)重新分配材料，得到最優(yōu)結(jié)構(gòu)拓撲。在大多數(shù)情況下，優(yōu)化目標是在給定的載荷模式和邊界條件下獲得結(jié)構(gòu)的最小重量。

假設(shè)選用各向同性材料作為研究對象，材料的泊松比取值是一個與材料密度等其他參數(shù)無關(guān)的常量，建立相對密度與彈性模量的函數(shù)關(guān)系為：

(1)

式中，ρi為第i個元素相對密度值，p為懲罰因子，E(ρi)為第i個元素的彈性模量，Emin表示孔洞部分中元素的彈性模量，E0為實體部分中元素的彈性模量。為了保證數(shù)值計算的穩(wěn)定性，通常取Emin=E0/1000，并且0

當ρi=1時，單元為實體結(jié)構(gòu)，當ρi=0時，單元為空洞結(jié)構(gòu)。通過控制懲罰因子p的值，實現(xiàn)元素密度向[0, 1]兩端快速收斂，從而得到接近理想狀態(tài)且材料最優(yōu)分布的離散結(jié)構(gòu)，從而得到最優(yōu)的材料分布。

變密度法假定每個單元的剛度矩陣依賴于懲罰因子p相對密度的改變，建立材料相對密度與材料彈性模量間顯示非線性的函數(shù)關(guān)系，如圖1所示。

圖1 變密度法密度插值函數(shù)模型

隨著懲罰因子p的增加，結(jié)構(gòu)剛度逐漸受到懲罰，對給定體積的材料，其中間密度值向兩端收斂，在單元網(wǎng)格的約束下拓撲優(yōu)化單元將重新分布，懲罰因子p一般是通過連續(xù)法從下界遞增到上界，每一步收斂后p遞增，以避免過早收斂到局部最小值。

在給定的體積約束條件下，選取最小柔度作為目標函數(shù)?；谧兠芏确▋?yōu)化問題可表述為：

find：ρ={ρ1,ρ2,···,ρn}T∈R

(2)

(3)

subject to：

(4)

式中，C(ρ)為給定拓撲的柔度，U是全局位移矢量，F(xiàn)是全局負載矢量，K是結(jié)構(gòu)全局剛度矩陣，N為元素個數(shù)，k0是單位楊氏模量的初始元素剛度矩陣，V(ρ)是優(yōu)化后結(jié)構(gòu)體積，V0是設(shè)計域的C初始體積，f是預(yù)先設(shè)定的體積分數(shù)，ρmin是一個包含最低允許相對密度的向量。

優(yōu)化準則法[11](Optimality Criteria，OC)在求解過程中具有收斂迅速，求解方式簡單，容易實現(xiàn)等特點。對于一定體積比約束下的柔度最小化問題，宜采用OC算法。OC算法的主要實現(xiàn)方式是利用目標函數(shù)與約束條件建立拉格朗日方程，將約束問題的約束與目標函數(shù)結(jié)合，成為無約束問題。通過對拉格朗日函數(shù)分析計算，可確定設(shè)計變量的優(yōu)化迭代準則可表示為：

(5)

式中，λ為當柔度取極值時全設(shè)計域中元素的應(yīng)變能密度。

2 敏度過濾方法

2.1 Sigmund敏度過濾方法

目前在變密度法拓撲優(yōu)化中普遍存在著棋盤格、網(wǎng)格依賴性現(xiàn)象。棋盤格現(xiàn)象，就是在不對敏度信息進行處理，使得拓撲優(yōu)化結(jié)果呈現(xiàn)單元之間的單點連接，形成單元密度為0、1交替排布的情形，導(dǎo)致不具備制造性，這種現(xiàn)象的存在是采用有限元求解無法克服的。網(wǎng)格依賴性和網(wǎng)格劃分大小具有密不可分的關(guān)系，不同的劃分結(jié)果對優(yōu)化結(jié)果產(chǎn)生較大影響，網(wǎng)格劃分越大，拓撲優(yōu)化結(jié)果越容易出現(xiàn)細長桿等微小結(jié)構(gòu)，導(dǎo)致優(yōu)化結(jié)果不具備良好的可制造性。

敏度過濾能夠有效解決數(shù)值不穩(wěn)定的問題，結(jié)構(gòu)的靈敏度分析在進行拓撲優(yōu)化時非常重要的，對模型的收斂判斷具有決定性的作用。結(jié)構(gòu)敏度可表示為：

(6)

目前常用的敏度過濾技術(shù)是Sigmund提出的敏度過濾方法，是一種局部約束方法，其本質(zhì)是通過人為設(shè)定一過濾半徑，在此范圍內(nèi)通過引入線性卷子因子將中心單元與其余單元之間的距離進行加權(quán)平均計算，修正目標函數(shù)的靈敏度，進而構(gòu)建這一范圍內(nèi)所有元素的敏度均值，更新敏度后續(xù)的迭代處理，從而解決棋盤格問題。Sigmund敏度過濾方法可表示為：

(7)

式中，Hei為卷積因子，rmin為最小過濾半徑，Δ(e,i)為元素i到單元元素的中心距離，由于設(shè)計變量的取值為[0, 1]，取ρmin=10-3，防止造成計算上的奇異性。

Sigmund敏度過濾方法存在如下問題：當所設(shè)定的過濾半徑較大時，拓撲優(yōu)化邊界會出現(xiàn)過度磨平的情況，這也導(dǎo)致優(yōu)化結(jié)果很容易出現(xiàn)邊界擴散等問題。

2.2 一種改進的敏度過濾方法

為解決傳統(tǒng)Sigmund敏度過濾方法中容易出現(xiàn)邊界擴散現(xiàn)象的問題，采用一種改進的敏度過濾方法，該方法在保留原有卷積因子的條件下，引入兩種非線性卷積因子Hin及Hat，其公式表示為：

(8)

(9)

式中，rmin及Δ(e,i)與Sigmund敏度過濾方法中卷積因子定義一致，λ為模型預(yù)先設(shè)定的網(wǎng)格最大數(shù)，即λ=max{x,y}。

由于Hin、Hat二者均為非線性卷積因子，能夠在根本上保證中心單元的敏度權(quán)值不會因過濾半徑的改變而改變，二者均為二階非線性函數(shù)，可以有效擴展卷積的擬合區(qū)域，使得在過濾半徑領(lǐng)域內(nèi)中心單元的權(quán)重值遠大于其余單元的權(quán)重值，確保中心單元敏度在迭代過程中不被過多的平均處理，保證最終拓撲優(yōu)化結(jié)果盡可能避免過度磨平的情況出現(xiàn)。

(10)

式中，α、β及γ分別為卷積因子Hi、Hin及Hat的修正權(quán)值，滿足α+β+γ=1，其控制各卷積因子對目標函數(shù)敏度值的影響程度。

采用一種帶有預(yù)設(shè)修正權(quán)值的方法，弱化邊界擴散的問題，可表示為：

(11)

式中，預(yù)設(shè)修正權(quán)值k=10-3，經(jīng)數(shù)值計算驗證取閾值η=0.7為宜，|ρi-ρe|為兩離散單元密度差值。

通過該方法的處理，能夠使拓撲優(yōu)化時具有黑白較為分明結(jié)果，與此同時在結(jié)構(gòu)內(nèi)部保證良好的均勻效果。

綜上可得改進的敏度過濾方法可表示為：

(12)

2.3 敏度過濾方法參數(shù)分析

為研究修正權(quán)值對拓撲優(yōu)化結(jié)果的影響，采用算例1所示的二維應(yīng)力結(jié)構(gòu)拓撲優(yōu)化確定本文方法的參數(shù)，為更準確地確定本文方法的參數(shù)，設(shè)定迭代終止準則的設(shè)計變量變化率ε=0.01，數(shù)值實驗結(jié)果見表1所示。

表1 不同參數(shù)計算結(jié)果

由表1數(shù)據(jù)可知，實驗1與實驗2所得結(jié)果顯示迭代次數(shù)較多，且柔度值較大；實驗6與實驗7所得結(jié)果顯示迭代次數(shù)較少，且柔度值較小，但是拓撲優(yōu)化結(jié)果容易出現(xiàn)微小孔洞。所以當α、β及γ滿足α=0.4～0.6、β=0.2～0.3、γ=0.2～0.3且滿足式α+β+γ=1時，該方法具有較好的算法合理性和迭代效率。

該方法可以有效控制棋盤格現(xiàn)象，解決邊界擴散問題，獲得邊界清晰的優(yōu)化結(jié)果，并大幅度提高了優(yōu)化速度，有效提升結(jié)構(gòu)剛度。

3 實驗分析驗證

采用多個拓撲優(yōu)化典型算例來驗證本改進的敏度過濾方法(以下簡稱本文方法)在不同條件下的有效性及可行性。算例通過MTLAB-2016a編程實現(xiàn)。在算例運行計算中，采用平面四結(jié)點雙線性正四邊形單元離散結(jié)構(gòu)，實體材料的彈性模量E=1，泊松比ν=1，選用最小過濾半徑rmin=1.5，設(shè)定迭代終止準則的設(shè)計變量變化率ε=0.01。

算例1：如圖2所示二維平面應(yīng)力結(jié)構(gòu)，該結(jié)構(gòu)設(shè)計區(qū)域為60 mm×40 mm，網(wǎng)格劃分為60×40，結(jié)構(gòu)左端采用全平面固定約束，結(jié)構(gòu)右端中間節(jié)點處受到載荷F=1的豎直向下載荷作用。

圖2 算例一模型示意圖

通過對比無敏度過濾方法、Sigmund敏度過濾方法、反三角函數(shù)因子方法[7]及指數(shù)因子方法[7]拓撲優(yōu)化結(jié)果，分析及驗證本文方法是否具有可行性，采用體積比為0.5的約束條件，懲罰因子p=3，得到多種處理方法下的拓撲優(yōu)化結(jié)果如表2所示。

表2 算例1不同處理方法拓撲優(yōu)化結(jié)果

分析表2優(yōu)化結(jié)果，無敏度過濾方法拓撲優(yōu)化結(jié)果呈現(xiàn)出復(fù)雜的棋盤格現(xiàn)象，采用后4種方法均能很好地抑制棋盤格的產(chǎn)生，結(jié)構(gòu)拓撲圖形清晰。Sigmund敏度過濾方法、反三角函數(shù)因子方法及指數(shù)因子方法均出現(xiàn)灰度單元的現(xiàn)象，采用本文方法得到的拓撲優(yōu)化結(jié)果在一定程度上避免了邊界存在灰度單元的現(xiàn)象。在相同的預(yù)設(shè)條件下，反三角函數(shù)因子方法、指數(shù)因子方法及本文方法在拓撲優(yōu)化時的迭代次數(shù)均在40左右，而采用本文方法時，優(yōu)化結(jié)構(gòu)具有最小的柔度值和良好的優(yōu)化效果。

算例2：如圖3所示經(jīng)典Michell結(jié)構(gòu)，該結(jié)構(gòu)設(shè)計區(qū)域為80 mm×40 mm，左下角節(jié)點處采用固定約束，右下角節(jié)點處采用鉸鏈約束，在圖示結(jié)構(gòu)底部中間節(jié)點處受到載荷F=1的豎直向下載荷作用。

圖3 算例二模型示意圖

分別對算例2采用Sigmund 敏度過濾方法、反三角函數(shù)因子方法、指數(shù)因子方法及本文方法，通過選用不同參數(shù)，分析及驗證后處理方法在不同網(wǎng)格計算密度、不同懲罰因子時是否具有可行性，采用體積比為0.5的約束條件，優(yōu)化結(jié)果如表3所示。

表3 不同優(yōu)化方法拓撲優(yōu)化結(jié)果

分析表3優(yōu)化結(jié)果，本文方法能夠很好地抑制數(shù)值不穩(wěn)定現(xiàn)象，在不同網(wǎng)格劃分、懲罰因子的條件下，本文方法均具有良好的表現(xiàn)，對比其他三種方法，本文方法在拓撲優(yōu)化時具有較少的迭代次數(shù)，與此同時優(yōu)化結(jié)構(gòu)的柔度收斂值最小，在一定程度上，避免了邊界存在中間單元的現(xiàn)象，有效驗證了本文方法的可行性。

圖4 算例三模型示意圖

算例3：如圖4所示二維平面應(yīng)力結(jié)構(gòu)，該結(jié)構(gòu)設(shè)計區(qū)域為60 mm×60 mm，網(wǎng)格劃分為60×60，結(jié)構(gòu)左端采用全平面固定約束，在圖示結(jié)構(gòu)右上角受到載荷F=1的豎直向上載荷作用，右下角受到載荷F=1的豎直向下載荷作用。

通過對比Sigmund敏度過濾方法、反三角函數(shù)因子方法及指數(shù)因子方法拓撲優(yōu)化結(jié)果，分析及驗證本文方法在多重載荷情況下是否具有可行性，采用體積比為0.4的約束條件，懲罰因子p=3，優(yōu)化結(jié)果如表4所示。

表4 算例3不同處理方法拓撲優(yōu)化結(jié)果

續(xù)表

分析表4優(yōu)化結(jié)果，本文方法在多重載荷情況下同樣能夠很好地抑制數(shù)值不穩(wěn)定現(xiàn)象，對比其他三種方法，本文方法在拓撲優(yōu)化時迭代次數(shù)穩(wěn)定，且優(yōu)化結(jié)構(gòu)的柔度收斂值最小，有效驗證了本文方法的在多負載工況下的可行性。

4 結(jié)論

變密度法作為處理拓撲優(yōu)化問題的主流方法之一，具有設(shè)計變量少、效率高等優(yōu)點。但傳統(tǒng)的變密度法在優(yōu)化過程中常出現(xiàn)如棋盤格現(xiàn)象等數(shù)值不穩(wěn)定現(xiàn)象，使得優(yōu)化模型提取較為困難，不具備良好的可制造性。本文在現(xiàn)有研究的基礎(chǔ)上，提出了一種改進的敏度過濾方法，實驗結(jié)果表明，該方法能有效消除棋盤格等問題，且有效提升迭代速度，降低柔度收斂值，提升結(jié)構(gòu)剛度。但是，本文方法在拓撲優(yōu)化結(jié)果的離散度方面還需進一步改進，這也是今后研究工作的主要方向。