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

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

在現(xiàn)代科學、經(jīng)濟和工程應用中，優(yōu)化顯得尤 為重要。得益于計算機科學的迅猛發(fā)展，各種復雜的優(yōu)化問題能夠在有限元環(huán)境中通過數(shù)學模型及算法進行模擬和快速求解。連續(xù)體結構拓撲優(yōu)化[1]已成為當前拓撲優(yōu)化研究方向的熱點問題之一，其在機械工程、土木工程、汽車輕量化等多領域具有普遍應用。目前，常用的拓撲優(yōu)化方法是均勻化方法[2]、變密度法[3]、水平集方法[4]、漸進結構優(yōu)化法[5]和ICM 方法[6]等。其中，變密度法因其設計變量少、程序簡單及應用范圍廣等特點被廣泛應用，其實質(zhì)是單元密度為0 到1 的離散變量之間的排列組合的問題[7]。在優(yōu)化過程中普遍存在如棋盤格現(xiàn)象、網(wǎng)格依賴性等數(shù)值不穩(wěn)定的問題[8]，制約了變密度法在結構拓撲優(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ù)目較多時，優(yōu)化結果會產(chǎn)生多余的細小結構。廉睿超等[14]提出一種考慮分區(qū)混合權重的敏度過濾方法，在網(wǎng)格劃分數(shù)目較多時，能有效抑制細小結構的產(chǎn)生，但迭代速度較慢。此外，龍凱等[15]提出一種考慮密度梯度的敏度過濾方法，但個別條件下結構開孔數(shù)目發(fā)生變化，影響最終拓撲結果。張志飛等[16]提出一種雙重SIMP 方法，但個別工況下仍存在灰度擴散現(xiàn)象。陳垂福等[17]提出一種變過濾半徑的敏度過濾方法，但不能有效控制單元的最小尺寸。

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

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

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

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

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

式中：ρ 為設計變量；C(ρ)為給定拓撲的柔度；U 是單元節(jié)點位移矢量；F 是元素節(jié)點負載向量；K 是全局剛度矩陣；N 為元素個數(shù)；k0是單位楊氏模量的元素剛度矩陣；V(ρ)是優(yō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ù)結合成為無約束問題。結合式（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)定的問題，結構的靈敏度分析在進行拓撲優(yōu)化時尤為重要，對模型的收斂判斷具有決定性的作用。結構敏度可表示為

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

圖3 Sigmund 敏度過濾示意圖

Sigmund 敏度過濾方法可表示為

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.3 方法參數(shù)確定

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

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

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

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

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

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

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

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

3 實驗驗證

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

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

圖10 算例1 模型示意圖

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

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

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

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

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

圖12 算例2 模型示意圖

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

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

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

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

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

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

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

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

圖15 算例3 模型示意圖

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

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

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

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

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

圖17 算例4 模型示意圖

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

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

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

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

4 結論

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