面向連續(xù)體結(jié)構(gòu)拓撲優(yōu)化的分區(qū)密度修正敏度過濾方法研究

張國鋒，徐雷，王鑫，李大雙

（四川大學 機械工程學院，成都 610065）

在現(xiàn)代科學、經(jīng)濟和工程應用中，優(yōu)化顯得尤 為重要。得益于計算機科學的迅猛發(fā)展，各種復雜的優(yōu)化問題能夠在有限元環(huán)境中通過數(shù)學模型及算法進行模擬和快速求解。連續(xù)體結(jié)構(gòu)拓撲優(yōu)化[1]已成為當前拓撲優(yōu)化研究方向的熱點問題之一，其在機械工程、土木工程、汽車輕量化等多領域具有普遍應用。目前，常用的拓撲優(yōu)化方法是均勻化方法[2]、變密度法[3]、水平集方法[4]、漸進結(jié)構(gòu)優(yōu)化法[5]和ICM 方法[6]等。其中，變密度法因其設計變量少、程序簡單及應用范圍廣等特點被廣泛應用，其實質(zhì)是單元密度為0 到1 的離散變量之間的排列組合的問題[7]。在優(yōu)化過程中普遍存在如棋盤格現(xiàn)象、網(wǎng)格依賴性等數(shù)值不穩(wěn)定的問題[8]，制約了變密度法在結(jié)構(gòu)拓撲優(yōu)化領域的發(fā)展。

文獻[9]提出一種考慮過濾半徑的敏度過濾方法，有效解決棋盤格、網(wǎng)格依賴性等問題，當所設定的過濾半徑較大時，拓撲優(yōu)化邊界會出現(xiàn)邊界擴散等問題，使其直接制造性不強。為了獲取清晰的拓撲優(yōu)化邊界，文獻[10]在原有敏度過濾方法的基礎上通過忽略密度項進行修正，但該方法存在迭代速度較慢，個別區(qū)域仍存在灰度單元的問題。諸多學者[11-12]通過引入Heaviside 函數(shù)的方法來控制邊界擴散現(xiàn)象，但迭代次數(shù)較多，效率較低。國內(nèi)學者對此也做了大量研究，羅震等[13]提出一種二重敏度過濾技術(shù)有效獲得邊界清晰的拓撲結(jié)構(gòu)，但在過濾半徑較小、單元數(shù)目較多時，優(yōu)化結(jié)果會產(chǎn)生多余的細小結(jié)構(gòu)。廉睿超等[14]提出一種考慮分區(qū)混合權(quán)重的敏度過濾方法，在網(wǎng)格劃分數(shù)目較多時，能有效抑制細小結(jié)構(gòu)的產(chǎn)生，但迭代速度較慢。此外，龍凱等[15]提出一種考慮密度梯度的敏度過濾方法，但個別條件下結(jié)構(gòu)開孔數(shù)目發(fā)生變化，影響最終拓撲結(jié)果。張志飛等[16]提出一種雙重SIMP 方法，但個別工況下仍存在灰度擴散現(xiàn)象。陳垂福等[17]提出一種變過濾半徑的敏度過濾方法，但不能有效控制單元的最小尺寸。

針對Sigmund 敏度過濾方法中邊界灰度問題，本文在現(xiàn)有研究的基礎上提出了一種分區(qū)密度修正的敏度過濾方法，將原敏度過濾區(qū)域進行劃分，采用新的復合卷積因子對不同區(qū)域處理，進一步自動判斷單元敏度絕對差值大小，采用密度修正權(quán)值，有效弱化邊界擴散的問題。該方法能有效獲得邊界清晰的拓撲優(yōu)化結(jié)果，并保證了求解的穩(wěn)定性，且提升迭代速度，降低柔度收斂值，提高結(jié)構(gòu)剛度。本文以柔度最小化為優(yōu)化目標，通過選取多個二維平面結(jié)構(gòu)算例驗證本文方法的可行性。

1 變密度法拓撲優(yōu)化理論

拓撲優(yōu)化是根據(jù)給定的負載及邊界條件對特定的設計空間進行優(yōu)化的一種數(shù)學梯度方法，利用最優(yōu)準則算法求解拓撲優(yōu)化問題，其目標是使系統(tǒng)柔度最小化。在連續(xù)體拓撲優(yōu)化中，首先確定設計域及設計變量，將設計區(qū)域按需求劃分為多個子區(qū)域，對其進行結(jié)構(gòu)分析以滿足預設優(yōu)化條件。如圖1 所示， ΩS為 實體區(qū)域， ΩV為空白區(qū)域，則設計域為Ω=ΩS+ΩV。

圖1 設計域與單元密度

固體各向同性材料的懲罰法（Solid isotropic microstructures with penalization，SIMP）即變密度法主要用于處理中間密度，是目前主流的基于有限元的復合連續(xù)體拓撲優(yōu)化方法。變密度法通常將元素相對密度作為設計變量，將0 到1 的整數(shù)優(yōu)化問題轉(zhuǎn)化為連續(xù)變量優(yōu)化問題。

對于變密度方法，拓撲變量被映射到密度變量中。選用各向同性材料，建立變密度法插值模型的表達式定義為

式中： ρi為單元i 的設計變量；p 為材料的懲罰因子，使設計變量 ρi趨近于0 或1； E(ρi)為第i 個單元的彈性模量； Emin為 孔洞元素的彈性模量； E0為實體元素的彈性模量；K 為插值處理后的剛度矩陣； ki為單元i 的剛度矩陣； ΔE=E0-Emin，為保證數(shù)值計算的穩(wěn)定性，通常 Emin=E0/1000 ， 且0

建立材料相對密度與材料彈性模量的函數(shù)關(guān)系，如圖2 所示。

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

在給定的體積約束條件下，選取最小柔度作為目標函數(shù)，基于變密度法優(yōu)化問題可表述為：

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

優(yōu)化準則法[18]（Optimality criteria，OC）在拓撲優(yōu)化中得到了廣泛的應用，對于單約束拓撲優(yōu)化而言，它是一種相當簡單且易于實現(xiàn)的算法。其設計變量迭代公式可表示為

式中：t 為平移限度， t =0.2；ω 為阻尼系數(shù)，ω=0.5。參數(shù)t 與ω 用來控制迭代穩(wěn)定及快速收斂，Be=為中間變量；ξ 為拉格朗日乘子，一般由二分法確定其值。

OC 算法的主要實現(xiàn)方式是利用目標函數(shù)與約束條件建立拉格朗日方程，將約束條件與目標函數(shù)結(jié)合成為無約束問題。結(jié)合式（3）所列目標函數(shù)與式（4）所列迭代設計變量，則優(yōu)化準則法可表示為

通過分析計算，可確定最終設計變量的優(yōu)化迭代準則可表示為

式中 ?為當柔度取極值時全設計域中元素的應變能密度。

當兩次連續(xù)的目標函數(shù)迭代差值的變化率小于預設收斂率時，優(yōu)化迭代終止，可表示為

2 分區(qū)密度修正的敏度過濾方法

2.1 Sigmund 敏度過濾方法

目前在變密度法拓撲優(yōu)化中普遍存在著棋盤格、網(wǎng)格依賴性現(xiàn)象。敏度過濾能夠有效解決數(shù)值不穩(wěn)定的問題，結(jié)構(gòu)的靈敏度分析在進行拓撲優(yōu)化時尤為重要，對模型的收斂判斷具有決定性的作用。結(jié)構(gòu)敏度可表示為

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

圖3 Sigmund 敏度過濾示意圖

Sigmund 敏度過濾方法可表示為

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

由式（9）可看出，在敏度區(qū)域范圍內(nèi)，任一單元到中心單元的距離越遠，其受到權(quán)重約束的程度越小，因?qū)蓡卧嚯x進行加權(quán)平均計算處理，使得拓撲優(yōu)化邊界產(chǎn)生過度磨平的情況，形成灰度單元帶，出現(xiàn)邊界擴散問題。如圖4 所示，當所設定的的過濾半徑較大，邊界含有過多低密度單元時，這種情況尤為明顯。

圖4 r min=5 時拓撲優(yōu)化結(jié)果

2.2 分區(qū)密度修正的敏度過濾方法

本文提出分區(qū)密度修正的敏度過濾方法來克服Sigmund 敏度過濾方法中邊界擴散的問題。其基本思想是將原敏度過濾區(qū)域劃分為新的子區(qū)域，對不同的區(qū)域采用不同的敏度過濾因子，保證在提高中心位置近距離單元快速收斂迭代的同時，降低外圍區(qū)域邊界處單元的權(quán)重值，進而降低遠距離單元對中心單元敏度的影響程度。同時采用一種帶有預設修正權(quán)值的方法，進一步弱化邊界擴散的問題。

如圖5 所示，將原過濾半徑為的區(qū)域進行劃分為Ⅰ、Ⅱ兩個子區(qū)域。其中Ⅰ區(qū)域的半徑為rm，采用常數(shù)1 作為敏度過濾因子，即采用原Sigmund 敏度過濾因子，保證單元中心在Ⅰ區(qū)域中心單元（紅色）的近距離單元（黃色）具有較快的收斂迭代速度。對單元中心在Ⅱ區(qū)域的遠距離單元（藍色）采用非線性函數(shù)因子處理，可以有效擴展卷積的擬合區(qū)域，以此降低該區(qū)域外圍邊界處單元的權(quán)重值，確保中心單元的敏度在迭代過程中不被過多的平均化，保證最終拓撲優(yōu)化結(jié)果盡可能避免過度磨平的情況出現(xiàn)。

圖5 敏度過濾區(qū)域劃分示意圖

分區(qū)域因子可表示為

采用一種帶有預設密度修正權(quán)值的方法，自動判斷單元敏度絕對差值大小，當差值大于預設閾值時，采用密度修正權(quán)值，進一步弱化邊界擴散的問題。可表示為

式中：q 為預設密度修正權(quán)值， q=10-3； |ρi-ρe|為兩離散單元密度差值 i ∈Ne； η為預設閾值，經(jīng)數(shù)值計算驗證發(fā)現(xiàn) η過大則導致迭代速度過慢， η過小則容易出現(xiàn)細小結(jié)構(gòu)，取η =0.8為宜。

綜上可得分區(qū)密度修正的敏度過濾方法可表示為

為考量拓撲優(yōu)化結(jié)果的收斂程度，Sigmund[10]最先提出引入離散率作為評價指標，其反映了中間單元密度值向0、1 兩端收斂的程度，可表示為

式中：當全部單元密度 ρi為0 或1 時， Md=0取得最小值，表示所有單元均已離散；當全部單元密度ρi=0.5時 ， Md=100取得最大值，表示所有單元均收斂。

為考量拓撲優(yōu)化結(jié)果中灰度單元占比程度，引入灰度率作為評價指標，可表示為

式中： ρj表示密度值介于0、1 之間的灰度單元； v0表示預設的體積分數(shù)。

2.3 方法參數(shù)確定

由圖5 可知，半徑rm的大小直接影響Ⅰ區(qū)域所含單元的個數(shù)，進而影響敏度過濾的效果。rm值過大時，因Ⅰ區(qū)域包含過多單元，并未有效弱化單元距離加權(quán)平均的效果，其趨向于文獻[10]所得結(jié)果；rm值過小時，因Ⅰ區(qū)域包含較少單元，迭代過程中極易出現(xiàn)出現(xiàn)數(shù)值不穩(wěn)定問題。由式（10）知rm取值與過濾半徑rmin具有一定關(guān)系，因此，為研究二者關(guān)系并確定影響因子β 的值，對算例1 所示的二維應力結(jié)構(gòu)進行拓撲優(yōu)化，取懲罰因子p=2，假定修正參數(shù)k=2、λ=0.5，得到rm與rmin在迭代次數(shù)、柔度值、離散率及灰度率的分析結(jié)果，如圖6 所示。

圖6 不同β 值多指標分析示意圖

由圖6 可知，當β=1 時，具有較好的迭代速率，但因離散率及灰度率過高，不宜采用；當β=5 時，在4 個指標中均有良好表現(xiàn)，但優(yōu)化結(jié)果容易出現(xiàn)數(shù)值不穩(wěn)定的問題，不宜采用。綜合多種因素考慮，取β 的值為3 附近時，具有良好的算法合理性及優(yōu)化結(jié)果。

當β=3 時rm=rmin-，分別取rmin=3 與rmin=5，探究分區(qū)域因子 He′i的 修正參數(shù)k 與λ 隨 Δ(e,i)的變化規(guī)律。假定k 一定，λ的變化曲線如圖7 所示；假定λ 一定，k 的變化曲線如圖8 所示；k 與λ 均變化的曲線如圖9 所示。

圖7 分區(qū)域因子修正參數(shù)k 一定λ 變化曲線

圖8 分區(qū)域因子修正參數(shù)λ 一定k 變化曲線

如圖7 變化曲線可知，當k、 Δ(e,i)一定時，λ 越小， He′i的值越大，曲線越陡峭；如圖8 變化曲線可知，當λ、 Δ (e,i)一 定時，k 越越小， He′i的值越小，曲線越平緩。修正參數(shù)k 與λ 過大或過小均不利于獲取較優(yōu)的拓撲優(yōu)化結(jié)果。如圖9 變化曲線及式（10）可知，針對不同過濾半徑rmin下，當 Δ(e,i)=1+rmin時He′i=0， 不同的修正參數(shù)得到的曲線所得到的 He′i的值在 Δ(e,i)=rmin處均大于0，有效避免了原過濾半徑rmin在某種程度上縮小，保證在rmin范圍內(nèi)各單元均參與分析計算，使得原過濾半徑rmin對優(yōu)化結(jié)構(gòu)的最小尺寸控制力，并保留了Sigmund 敏度過濾方法求解計算時的穩(wěn)定性，確保了拓撲優(yōu)化結(jié)果的可靠性，并采用一種帶有預設密度權(quán)值的方法，進一步的弱化了邊界擴散問題。

圖9 分區(qū)域因子修正參數(shù)λ 與k 均變化曲線

為進一步研究修正參數(shù)k 與λ 對拓撲優(yōu)化結(jié)果的影響，采用控制變量法仍對算例1 所示的結(jié)構(gòu)進行拓撲優(yōu)化來確定其值，取懲罰因子 p=3，過濾半徑rmin=3， 影響因子 β=3，數(shù)值實驗結(jié)果見表1 所示。

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

由表1 可知，不同參數(shù)組合所得結(jié)果的迭代次數(shù)、柔度值、離散率及灰度率有較大差異。綜合考慮多種因素，設定參數(shù) k=2 、 λ=0.5時，該方法具有較好的算法合理性和迭代效率。參數(shù)設定并不唯一，可根據(jù)具體研究對象略加調(diào)整。

3 實驗驗證

采用多個拓撲優(yōu)化典型算例來驗證分區(qū)密度修正的敏度過濾方法（以下簡稱本文方法）在不同條件下的有效性及可行性，算例通過MTLAB-2016a 編程實現(xiàn)。在本文算例運行計算中，采用平面四結(jié)點雙線性正四邊形單元離散結(jié)構(gòu)，實體材料的彈性模量 E =1， 泊松比ν =1， 懲罰因子 p=3[19]。

算例1 如圖10 二維平面應力結(jié)構(gòu)，該結(jié)構(gòu)設計區(qū)域為120 mm×30 mm，左下角節(jié)點處采用固定約束，右下角節(jié)點處采用鉸鏈約束，在圖示結(jié)構(gòu)上部中間節(jié)點處受到載荷F=1 的豎直向下載荷作用。該模型屬于對稱結(jié)構(gòu)，為減小計算量，故采用模型一半結(jié)構(gòu)進行分析。網(wǎng)格計算密度取60×30，過濾半徑取rmin=3。

圖10 算例1 模型示意圖

通過對比變密度法、文獻[10, 14-16]方法及本文方法的拓撲優(yōu)化結(jié)果，分析及驗證本文方法是否具有可行性。采用體積比為0.5 的約束條件，得到多種處理方法下的拓撲優(yōu)化結(jié)果如圖11 所示，不同指標對比結(jié)果如表2 所示。

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

圖11 算例1 不同優(yōu)化方法結(jié)果對比

對比上述實驗結(jié)果可看出，除變密度法外，文獻[10，14-16]方法及本文方法均能夠獲得邊界清晰的拓撲優(yōu)化結(jié)構(gòu)，有效抑制邊界擴散現(xiàn)象。然而，文獻[10]方法迭代速率較為緩慢，文獻[14]方法優(yōu)化結(jié)果邊界仍然存在個別數(shù)值不穩(wěn)定現(xiàn)象，文獻[15]方法不能保證其結(jié)構(gòu)完整性，文獻[16]方法優(yōu)化結(jié)果仍存在部分灰度單元。對比上述方法，本文方法在有效獲得清晰拓撲結(jié)構(gòu)的同時，具有較快的收斂速度和較低的柔度值，對網(wǎng)格離散的控制較好，能有效控制結(jié)構(gòu)，確保其完整性。

算例2 如圖12 所示二維平面應力結(jié)構(gòu)，該結(jié)構(gòu)設計區(qū)域為120 mm×80 mm，結(jié)構(gòu)左端采用全平面固定約束，結(jié)構(gòu)右端中間節(jié)點處受到載荷F=1 的豎直向下載荷作用。

圖12 算例2 模型示意圖

通過選用不同參數(shù)，分析及驗證本文方法在不同網(wǎng)格計算密度、不同體積比時是否具有可行性，采用體積比為0.5 的約束條件。拓撲優(yōu)化結(jié)果如圖13 所示，不同指標對比結(jié)果如表3 所示。

圖13 算例2 拓撲優(yōu)化優(yōu)化結(jié)果

表3 算例2 優(yōu)化結(jié)果

由上述實驗b、c、d 在不同過濾半徑條件下的優(yōu)化結(jié)果對比可知，本文方法在不同指標下均具有良好的結(jié)果，在過濾半徑較大時，有效抑制邊界擴散問題。由實驗b、e、f 在不同網(wǎng)格單元劃分條件下的優(yōu)化結(jié)果可知，本文方法在網(wǎng)格劃分數(shù)目較多時，仍具有良好的結(jié)構(gòu)性能。

為進一步驗證本文方法是否具有穩(wěn)定性，通過對比變密度法與本文方法，將實驗a 與b 的柔度值進行數(shù)據(jù)處理。為清晰顯示比較結(jié)果，從第五次迭代開始記錄數(shù)據(jù)，結(jié)果如圖14 所示。

圖14 算例2 不同方法迭代曲線示意圖

由圖14 可知，本文方法較之變密度法具有更快的收斂速度，同時在結(jié)構(gòu)柔度方面具有良好的穩(wěn)定性。

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

圖15 算例3 模型示意圖

通過對比變密度法與本文方法拓撲優(yōu)化結(jié)果，進一步分析及驗證本文方法在多重載荷情況下是否具有可行性，采用體積比為0.4 的約束條件，過濾半徑取rmin= 2。拓撲優(yōu)化結(jié)果如圖16 所示，不同指標對比結(jié)果如表4 所示。

圖16 算例3 拓撲優(yōu)化優(yōu)化結(jié)果

表4 算例3 優(yōu)化結(jié)果

由上述實驗結(jié)果可知，本文方法在多重載荷情況下同樣具有良好的可行性。其拓撲優(yōu)化迭代次數(shù)穩(wěn)定，優(yōu)化結(jié)構(gòu)的柔度收斂值較小，能有效控制離散率及灰度率指標。

算例4 圖17 所示為二維平面應力結(jié)構(gòu)，該結(jié)構(gòu)設計區(qū)域為160 mm×80 mm，設計區(qū)域中兩孔洞圓直徑均為40 mm，左孔洞圓中心距左端面為45 mm，距底面為40 mm，右孔洞圓中心距右端面為45 mm，距底面為40 mm，網(wǎng)格劃分為160×80。結(jié)構(gòu)左下角節(jié)點處采用固定約束，右下角節(jié)點處采用鉸鏈約束，在圖示結(jié)構(gòu)頂部中間節(jié)點處受到載荷F = 1 的豎直向下載荷作用。

圖17 算例4 模型示意圖

通過對比變密度法與本文方法拓撲優(yōu)化結(jié)果，進一步分析及驗證本文方法在含有多個孔洞結(jié)構(gòu)的情況下是否具有可行性，采用體積比為0.5 的約束條件，過濾半徑取rmin= 2。拓撲優(yōu)化結(jié)果如圖18所示，不同指標對比結(jié)果如表5 所示。

圖18 算例4 拓撲優(yōu)化優(yōu)化結(jié)果

表5 算例4 優(yōu)化結(jié)果

由上述實驗結(jié)果可知，本文方法在含有多個孔洞結(jié)構(gòu)的情況下同樣具有良好的可行性。其拓撲優(yōu)化迭代次數(shù)穩(wěn)定，優(yōu)化結(jié)構(gòu)的柔度收斂值較小，能有效控制離散率及灰度率指標。對比圖18 優(yōu)化結(jié)果，較之變密度法，本文方法在一定程度上能夠有效解決網(wǎng)格依賴性的問題。

4 結(jié)論

變密度法在拓撲優(yōu)化過程中常出現(xiàn)棋盤格現(xiàn)象等數(shù)值不穩(wěn)定問題，使優(yōu)化結(jié)果常具有灰度單元，致使優(yōu)化模型提取較為困難。雖然Sigmund 敏度過濾方法能夠有效減少數(shù)值不穩(wěn)定問題，但該方法容易導致邊界擴散，不具備良好的可制造性。為解決這一問題，本文在現(xiàn)有研究的基礎上提出了一種分區(qū)密度修正的敏度過濾方法。實驗結(jié)果表明，該方法能有效避免邊界擴散等問題，迭代速度穩(wěn)定，離散率及灰度率指標較低，同時能夠降低柔度收斂值，提升結(jié)構(gòu)剛度，保證結(jié)構(gòu)完整性。當前，該方法僅適用于二維平面結(jié)構(gòu)的拓撲優(yōu)化，在今后的研究中，需進一步探究其在三維空間結(jié)構(gòu)中的應用途徑，進一步提升方法的適用性及有效性。