一種用于密度法的SR插值模型

劉曉晨，孫 宇 ，敬石開，郄龍飛

（1.南京理工大學(xué) 機(jī)械工程學(xué)院，南京 210094；2.北京理工大學(xué) 機(jī)械與車輛學(xué)院，北京 100029）

0 引言

拓?fù)鋬?yōu)化是一種在給定設(shè)計(jì)區(qū)域內(nèi)，根據(jù)邊界及載荷條件，在滿足特定限制要求的情況下實(shí)現(xiàn)結(jié)構(gòu)性能最優(yōu)的先進(jìn)設(shè)計(jì)方法，其屬于結(jié)構(gòu)優(yōu)化設(shè)計(jì)的一種[1]。自1988年Bendsoe和Kikuchi[2]提出均勻化方法以來，拓?fù)鋬?yōu)化得到快速發(fā)展，目前已發(fā)展出多種類型的優(yōu)化方法，如密度法[3]、水平集法[4]、漸進(jìn)結(jié)構(gòu)法[5]、組件優(yōu)化法[6]、拓?fù)鋵?dǎo)數(shù)法[7]等。

密度法是常用的一種拓?fù)鋬?yōu)化方法，其基本思想是人為的引入一種假想的取值為{0,1}的密度變量，其中0代表空相材料，由白色表示，1代表實(shí)體相材料，由黑色表示。建立拓?fù)鋬?yōu)化模型時(shí)，以材料密度為設(shè)計(jì)變量，通過單元密度值來控制單元的增刪，實(shí)現(xiàn)拓?fù)浣Y(jié)構(gòu)的改變。

因設(shè)計(jì)變量只能取值0或1，上述拓?fù)鋬?yōu)化問題是一個(gè)典型的離散變量拓?fù)鋬?yōu)化問題。離散變量拓?fù)鋬?yōu)化問題是指設(shè)計(jì)變量在優(yōu)化過程中只能取某些離散值，而不是在一個(gè)區(qū)間內(nèi)連續(xù)變化[8]。因設(shè)計(jì)變量不連續(xù)，導(dǎo)致目標(biāo)函數(shù)和約束條件不可微，因此離散變量拓?fù)鋬?yōu)化問題在求解時(shí)只能借助組合數(shù)學(xué)方法，而不能使用解析數(shù)學(xué)方法[9]。但拓?fù)鋬?yōu)化的典型特征就是設(shè)計(jì)變量規(guī)模大，密度法中設(shè)計(jì)變量數(shù)量與有限元網(wǎng)格劃分?jǐn)?shù)量相等，均勻化方法中設(shè)計(jì)變量數(shù)量是有限元網(wǎng)格劃分?jǐn)?shù)量的3倍。這帶來的問題就是“組合爆炸”，某些情況即使借助計(jì)算機(jī)也難以求解。

為了避免離散變量拓?fù)鋬?yōu)化中的上述問題，通常對(duì)設(shè)計(jì)變量進(jìn)行松弛，將設(shè)計(jì)變量取值范圍由{0,1}放松到[0,1]，允許其在0到1之間連續(xù)取值。此舉有效緩解了優(yōu)化求解過程中的組合爆炸問題，但是將密度值松弛到[0,1]意味著在優(yōu)化過程中引入了中間密度變量（即處于0和1之間的密度變量），中間密度變量沒有物理意義，并且使得優(yōu)化結(jié)果中產(chǎn)生大量的灰度單元（介于黑與白之間的灰色單元）。

為了抑制中間密度單元的產(chǎn)生，密度法中通常采用懲罰策略對(duì)中間密度變量進(jìn)行處理，迫使其向0或者1兩端靠攏，從而得到清晰的0-1解，此類方法又稱為材料插值模型。密度法中最有代表性的材料插值模型是SIMP模型[10,11]和RAMP模型[12]。SIMP模型通過冪指數(shù)的形式對(duì)中間密度單元進(jìn)行懲罰，RAMP模型通過有理數(shù)的形式對(duì)中間密度單元進(jìn)行懲罰，其與SIMP模型的主要區(qū)別是設(shè)計(jì)變量在0處的敏度不為0，這有助于抑制低密度值單元的產(chǎn)生，并提高優(yōu)化過程穩(wěn)定性。但RAMP模型采用有理數(shù)的形式進(jìn)行懲罰，其優(yōu)化效率遠(yuǎn)低于SIMP模型。

總的來說，SIMP模型具有較高的優(yōu)化效率，RAMP模型具有較高的穩(wěn)定性。但二者仍處于割裂狀態(tài)，為了充分利用二者的優(yōu)勢，可以構(gòu)建統(tǒng)一的懲罰模型。統(tǒng)一模型一般采用加權(quán)的形式，這種處理方式在多相材料中已有出現(xiàn)[13]：對(duì)于沒有空白相的兩相材料結(jié)構(gòu)，可以使用插值，該插值對(duì)于每個(gè)材料相采用Hashin-Shtrikman上限和下限的加權(quán)平均值方式進(jìn)行處理；對(duì)于帶有空白相的兩相材料結(jié)構(gòu)，則使用SIMP模型和Hashin-Shtrikman模型進(jìn)行混合插值處理[14]。

受此啟發(fā)，本文首先提出一種SIMP模型和RAMP模型加權(quán)的插值模型，即SR模型。SR模型通過加權(quán)因子將SIMP和RAMP向結(jié)合，通過調(diào)整加權(quán)因子的取值，SR模型可以轉(zhuǎn)換為SIMP模型或者RAMP模型；其次，在SR模型框架下，構(gòu)建了以機(jī)械柔順度為目標(biāo)函數(shù)，體積約束為約束條件的拓?fù)鋬?yōu)化模型；最后，研究了適用于該優(yōu)化模型的優(yōu)化準(zhǔn)則求解算法，以及SR模型在不同加權(quán)因子取值策略下的懲罰效果。

1 SR模型構(gòu)建基礎(chǔ)

SIMP模型和RAMP模型均屬于隱式懲罰法，其直接對(duì)設(shè)計(jì)變量進(jìn)行懲罰。SIMP模型采用冪指數(shù)的形式對(duì)中間密度單元進(jìn)行懲罰，單元i的楊氏模量Ei可表示為：

其中，E0表示固體相的楊氏模量，Emin表示空白相的楊氏模量，其被賦予一個(gè)非常小的值以防止優(yōu)化過程中剛度矩陣變得奇異，xi為設(shè)計(jì)變量密度，p為懲罰因子。圖1展示了懲罰因子p在不同取值下的SIMP插值模型的變化曲線。

圖1 SIMP插值模型曲線圖

圖1中X坐標(biāo)為設(shè)計(jì)變量密度，Y坐標(biāo)為單元楊氏模量。如圖所示，當(dāng)懲罰因子p取值為1時(shí)，SIMP模型為一條直線，此時(shí)SIMP模型對(duì)中間密度單元沒有懲罰力度，目標(biāo)函數(shù)為凸函數(shù)，優(yōu)化模型有唯一的全局最優(yōu)解；隨著懲罰因子p取值增加，SIMP模型對(duì)中間密度單元的懲罰力度逐漸加強(qiáng)，目標(biāo)函數(shù)由凸函數(shù)轉(zhuǎn)為非凸函數(shù)。從圖1中可以看出，當(dāng)懲罰因子p取值大于20時(shí)，密度值處于0.8以下的單元基本被歸為0。

RAMP模型采用有理數(shù)的形式對(duì)中間密度單元進(jìn)行懲罰，單元i的楊氏模量Ei可表示為：

其中，q為懲罰因子。圖2展示了懲罰因子取值相同時(shí)，SIMP模型和RAMP模型懲罰力度的對(duì)比情況。

圖2 SIMP模型和RAMP模型對(duì)比圖

圖2中X坐標(biāo)為設(shè)計(jì)變量密度，Y坐標(biāo)為單元楊氏模量。SIMP模型懲罰曲線由實(shí)線表示，RAMP模型懲罰曲線由虛線表示。區(qū)別于SIMP模型，當(dāng)懲罰因子取值為1時(shí)，RAMP模型對(duì)中間密度單元有懲罰力度，對(duì)于其他懲罰因子取值，RAMP模型的懲罰力度均低于SIMP模型。

2 SR模型構(gòu)建

假設(shè)最優(yōu)拓?fù)錁?gòu)型由?s表示，由密度法的定義可知，?s是設(shè)計(jì)域?的子域，在?s上設(shè)計(jì)變量數(shù)值為1，在?/?s區(qū)域設(shè)計(jì)變量值為0。如果用E1和E2分別表示兩種材料的楊氏模量，則第i個(gè)單元的楊氏模量需要滿足：

其中，ΔE=E1-E2。式(3)中，假設(shè)材料1的泊松比與材料2的泊松比相同。在密度法中，固體相被視為材料E1，空白相被視為材料E2，則插值模型需滿足：

基于此，SR模型中楊氏模量和設(shè)計(jì)變量的關(guān)系式為：

其中，r為懲罰因子，θ為加權(quán)因子。加權(quán)因子θ決定SR模型中SIMP模型和RAMP模型的權(quán)重：當(dāng)θ=0時(shí)，SR模型轉(zhuǎn)化為RAMP模型；當(dāng)θ=1時(shí)，SR模型轉(zhuǎn)化為SIMP模型。

圖3展示了懲罰因子r取值為5時(shí)的SR模型的插值曲面圖。

圖3 SR插值模型曲面圖

圖3中，X坐標(biāo)為設(shè)計(jì)變量，Y坐標(biāo)為加權(quán)因子，Z坐標(biāo)為楊氏模量，懲罰因子r取值為5。當(dāng)懲罰因子取不同數(shù)值時(shí)，區(qū)別于SIMP模型和RAMP模型的曲線圖，SR模型為一系列不同的曲面。如上所述，當(dāng)懲罰因子取值為0時(shí)，SR模型轉(zhuǎn)化為RAMP模型；當(dāng)懲罰因子取值為0時(shí)，SR模型轉(zhuǎn)化為SIMP模型。在圖3中，θ=0時(shí)SR模型的曲線與RAMP模型中q=5的曲線相同；θ=1時(shí)SR模型的曲線與SIMP模型中p=5的曲線相同。因此，r=5的SR模型曲面為q=5的RAMP模型過渡到p=5的SIMP模型所形成。

3 基于SR模型的拓?fù)鋬?yōu)化模型

對(duì)于以機(jī)械柔順度為目標(biāo)函數(shù)，體積約束為約束條件，采用SR插值模型的拓?fù)鋬?yōu)化模型可表示為：

其中，c=UTKU為目標(biāo)函數(shù)機(jī)械柔順度，K，U和F為別表示全局剛度矩陣、全局位移矩陣和全局載荷矩陣，ui表示單元i的位移，ki為單元?jiǎng)偠染仃?，k0表示單位楊氏模量的單元?jiǎng)偠染仃嚕珿0為體積約束，V表示優(yōu)化后設(shè)計(jì)域的體積，V0表示原始設(shè)計(jì)域的體積，vi表示單元i的體積，n表示有限元網(wǎng)格劃分的數(shù)量，xmin為一很小的正數(shù)，為了避免優(yōu)化過程整體剛度矩陣變得奇異。

4 SR拓?fù)鋬?yōu)化模型求解

根據(jù)OC算法的原理[15]，用于求解SR拓?fù)鋬?yōu)化模型的優(yōu)化準(zhǔn)則算法可以通過優(yōu)化模型中目標(biāo)函數(shù)和約束條件所構(gòu)建的拉格朗日函數(shù)進(jìn)行推導(dǎo)。公式(6)對(duì)應(yīng)拓?fù)鋬?yōu)化問題的拉格朗日函數(shù)為：

當(dāng)xi=xmin時(shí)，設(shè)計(jì)變量上限無效，此時(shí)；當(dāng)xi=xmax時(shí)，設(shè)計(jì)變量下限無效，此時(shí)；當(dāng)時(shí)，設(shè)計(jì)變量的上下限均無效，此時(shí)?；诖?，式(8)可轉(zhuǎn)變?yōu)椋?/p>

將機(jī)械柔順度c=UTKU代入式(9)，考慮等式項(xiàng)，可得：

在SR模型中，單元?jiǎng)偠染仃噆i和單位楊氏模量剛度矩陣k0關(guān)系為：

將其代入式(10)，并結(jié)合整體剛度矩陣的對(duì)稱性，可得：

OC算法的優(yōu)化設(shè)計(jì)準(zhǔn)則表示為：

將設(shè)計(jì)變量的上下限代入式(14)中，可得OC算法的迭代優(yōu)化公式：

其中，引入阻尼系數(shù)η以保證優(yōu)化過程的穩(wěn)定性和收斂性。一般地，當(dāng)滿足以下條件時(shí)，OC算法停止迭代：

其中，ε為任意小的正數(shù)。

5 數(shù)值算例

本節(jié)以MBB梁為對(duì)象測試SR模型在不同懲罰因子取值策略下的優(yōu)化效果。采用敏度過濾器進(jìn)行過濾。MBB梁的設(shè)計(jì)域、邊界條件及載荷情況如圖4所示（由于對(duì)稱性，圖中僅展示一半MBB梁），其左端面為對(duì)稱結(jié)構(gòu)，左上角作用一個(gè)大小為1，垂直向下的力F，右下角處為水平支撐，設(shè)計(jì)域的長寬比例為3:1。參數(shù)設(shè)置如下：網(wǎng)格劃分120×40，體積分?jǐn)?shù)取值為0.5，懲罰因子取值為4，過濾半徑取值為4。

圖4 MBB梁的設(shè)計(jì)域、邊界條件及載荷情況

為了評(píng)價(jià)優(yōu)化結(jié)果的0/1離散率，本節(jié)引入指標(biāo)Mnd指標(biāo)[16]，其表達(dá)式為：

圖5展示的是SIMP模型和RAMP模型的優(yōu)化結(jié)果。

圖5 MBB梁的優(yōu)化結(jié)果

圖5(a)為SIMP模型優(yōu)化結(jié)果，其右上部分和左下部分支撐處出現(xiàn)明顯的灰度單元聚集區(qū)域。RAMP模型的優(yōu)化結(jié)果較為清晰。兩者的優(yōu)化結(jié)果中均未出現(xiàn)棋盤格和網(wǎng)格依賴問題。

圖6展示的SR模型θ取值為0.25，0.5和0.75的優(yōu)化結(jié)果。

圖6 MBB梁的優(yōu)化結(jié)果

由圖6可以看出，SR模型的優(yōu)化結(jié)果中在上下兩支撐處未出現(xiàn)明顯的灰度單元聚集區(qū)域，并且在非支撐區(qū)域也存在較少中間密度單元，優(yōu)化結(jié)果清晰，沒有產(chǎn)生網(wǎng)格依賴和棋盤格現(xiàn)象。上述測試的優(yōu)化數(shù)據(jù)如表1所示。

表1 MBB梁優(yōu)化數(shù)據(jù)

SIMP模型采用冪指數(shù)的懲罰形式，相同參數(shù)下，其收斂速度最快，迭代次數(shù)最少。當(dāng)θ取值為0.25時(shí)，SR模型體現(xiàn)出SIMP模型的特征，具有較快的收斂速度，且目標(biāo)函數(shù)值和Mnd值均低于SIMP優(yōu)化結(jié)果。隨著加權(quán)因子取值增大，SR模型迭代次數(shù)不斷增加。其原因是SR模型逐漸由SIMP模型過渡到RAMP模型，其懲罰效率逐漸降低。三種加權(quán)因子取值下SR模型優(yōu)的迭代次數(shù)均低于RAMP模型的迭代次數(shù)。

目標(biāo)函數(shù)值也體現(xiàn)出同樣的規(guī)律：SR模型θ=0.25優(yōu)化結(jié)果的目標(biāo)函數(shù)值最低，目標(biāo)函數(shù)值隨著加權(quán)因子取值增大而增大。需要指出的是，θ=0.25時(shí)SR模型目標(biāo)函數(shù)值低于SIMP模型，θ=0.75時(shí)目標(biāo)函數(shù)值高于RAMP模型，這是SR模型與SIMP和RAMP模型的差異。

0-1離散率方面，SIMP模型因其優(yōu)化結(jié)果中上下兩支撐處產(chǎn)生灰度單元聚集區(qū)域，其Mnd值較高。SR模型和RAMP模型的Mnd值明顯低于SIMP模型Mnd值。加權(quán)因子取越大值時(shí)，SR模型優(yōu)化結(jié)果的Mnd值越低。

6 結(jié)論

密度法中通過對(duì)設(shè)計(jì)變量進(jìn)行松弛來避免離散變量拓?fù)鋬?yōu)化問題求解時(shí)的組合爆炸問題，為了抑制因松弛操作而引入的中間密度單元，密度法中通常引入材料插值模型，對(duì)中間密度單元進(jìn)行懲罰，使其向0或1兩端靠攏。

本文首先提出一種SIMP模型和RAMP模型的加權(quán)材料插值模型，即SR模型。該模型通過引入加權(quán)因子將SIMP與RAMP相結(jié)合，通過調(diào)整加權(quán)因子的取值，SR模型可轉(zhuǎn)化為SIMP模型或者RAMP模型?；诖?，本文分析了SR模型的懲罰曲面特征，并構(gòu)建了以機(jī)械柔順度為目標(biāo)函數(shù)，體積約束為約束條件，采用SR插值模型的拓?fù)鋬?yōu)化模型；其次，研究了適用于該模型的優(yōu)化準(zhǔn)則求解算法。以MBB梁為例，驗(yàn)證SIMP模型、RAMP模型以及SR模型的優(yōu)化特征，研究加權(quán)因子對(duì)優(yōu)化結(jié)果的影響關(guān)系。

實(shí)驗(yàn)表明，SR模型可獲得較為清晰的優(yōu)化構(gòu)型，優(yōu)化結(jié)果中無棋盤格和網(wǎng)格依賴現(xiàn)象，并且沒產(chǎn)生明顯的中間密度單元聚集區(qū)域。通過調(diào)整加權(quán)因子取值，SR模型表現(xiàn)出SIMP模型收斂速度快以及RAMP穩(wěn)定性高的特點(diǎn)；通過分配較大的懲罰因子，SR模型可獲得低Mnd值的優(yōu)化結(jié)果，并且優(yōu)化過程沒有出現(xiàn)震蕩。