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

張國(guó)鋒，徐 雷，李大雙，余方超

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

0 引言

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

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

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

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

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

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

(1)

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

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

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

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

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

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

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

(2)

(3)

subject to：

(4)

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

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

(5)

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

2 敏度過(guò)濾方法

2.1 Sigmund敏度過(guò)濾方法

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

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

(6)

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

(7)

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

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

2.2 一種改進(jìn)的敏度過(guò)濾方法

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

(8)

(9)

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

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

(10)

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

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

(11)

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

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

綜上可得改進(jìn)的敏度過(guò)濾方法可表示為：

(12)

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

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

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

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

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

3 實(shí)驗(yàn)分析驗(yàn)證

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

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

圖2 算例一模型示意圖

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

表2 算例1不同處理方法拓?fù)鋬?yōu)化結(jié)果

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

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

圖3 算例二模型示意圖

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

表3 不同優(yōu)化方法拓?fù)鋬?yōu)化結(jié)果

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

圖4 算例三模型示意圖

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

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

表4 算例3不同處理方法拓?fù)鋬?yōu)化結(jié)果

續(xù)表

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

4 結(jié)論

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