基于SIMP插值模型的棋盤格現(xiàn)象濾波控制算法研究

張 宇，羅陸鋒

（天津職業(yè)技術師范大學 機械工程學院，天津 300222）

0 引言

棋盤格是連續(xù)體結構拓撲優(yōu)化中較為常見的一種數(shù)值不穩(wěn)定現(xiàn)象[1,2]。所謂棋盤格，是指在結構優(yōu)化過程中某些區(qū)域出現(xiàn)的單元材質(zhì)密度周期性分布的一種現(xiàn)象。棋盤格的出現(xiàn)與分析單元的選擇有關，與單元的類型、單元的劃分、單元彈性模量與單元密度之間的關系以及優(yōu)化算法本身等都有關系。除了棋盤格以外，在拓撲優(yōu)化的數(shù)值計算中一般還存在著網(wǎng)格依賴性、局部極值、多孔材料等問題。棋盤格和網(wǎng)格依賴兩種現(xiàn)象一般同時出現(xiàn)在優(yōu)化結果中，能夠有效去除棋盤格的方法通常也能有效克服網(wǎng)格依賴現(xiàn)象[3]。

拓撲優(yōu)化中數(shù)值計算不穩(wěn)定性的消除非常重要，它關系到數(shù)值計算的收斂性和計算結果的可制造性問題。解決棋盤格除了采用高階有限單元代替低階單元方法外，一般還有周長約束、濾波法等。從算法本身來說，周長約束等于在原有的模型中添加了約束類型，可能造成收斂性和數(shù)值計算問題。濾波法是圖像處理技術中的常用方法，其基本思想是將原函數(shù)與濾波函數(shù)進行卷積運算[3]，對原函數(shù)進行規(guī)整化處理，提高函數(shù)光滑性。該思想應用到拓撲優(yōu)化中的具體表現(xiàn)，就是通過調(diào)整每次迭代后單元的靈敏度從而避免棋盤格的出現(xiàn)。在這種方法中，某一單元的靈敏度依賴于其自身以及濾波半徑范圍內(nèi)相鄰單元靈敏度的加權平均值。其優(yōu)點是不需要在優(yōu)化問題中加入額外約束，另

外容易實施，收斂性好，計算穩(wěn)定。這種方法對消除棋盤格問題非常有效，也能一并解決網(wǎng)格依賴性問題。

1 基于SIMP插值模型的濾波算法

1.1 靈敏度分析

SIMP[1,4,5]方法的思想和前提是將離散單元內(nèi)部的材料屬性定義為常數(shù)，設計變量定義為離散單元的相對密度，盡量減少結構中間密度單元的數(shù)目，使結構單元密度盡可能為0或1[6]，從而用連續(xù)優(yōu)化設計方法來近似離散優(yōu)化設計[7]。

根據(jù)變密度方法的思想，單元材料為各向同性材料，假設單元內(nèi)材料的宏觀性質(zhì)為各向同性，其泊松比與實體材料的泊松比相等。單元剛度矩陣與單元材料的等效性質(zhì)相關[8]。為了控制整個結構位移，設目標函數(shù)為：

其中，C( ρ)為結構柔度，F(xiàn)和U分別為節(jié)點載荷陣和位移陣，K為剛度矩陣。如果我們以最大位移為約束條件，根據(jù)變密度方法的思想，當對連續(xù)體結構進行有限元離散以后，假設位移最大的節(jié)點的位移受到了限制，則約束條件如下：

式中u*為根據(jù)有限元分析的結果所給定的合理的位移上限，ρi代表單元密度，ε為人為給定的密度下限，i=j1, j2,…, jk指優(yōu)化后密度保持不變的單元號。

有限元平衡方程為：

在插值前和插值后的材料彈性模量之間引入關系式：

其中，EP代表插值以后的彈性模量，ρiP代表迭代后的單元密度，E0和Emin分別代表固體和去除部分的彈性模量，彈性模量E( ρ)與實體彈性模量E0的比值同相對密度近似成指數(shù)關系，對上式進行簡化：則目標函數(shù)可以修改為：對材料的平衡方程進行微分：該方程與目標函數(shù)相結合，可得：

根據(jù)以上各式，容易得到結構總體柔度的靈敏度等效形式：

1.2 濾波算法實現(xiàn)

濾波法本來是圖像處理技術中的常用方法，其數(shù)學表達式為：

式中，F(xiàn)代表原函數(shù)，G代表濾波函數(shù)，該表達式的基本思想是將原函數(shù)與濾波函數(shù)進行卷積，對原函數(shù)進行規(guī)整化處理，提高函數(shù)光滑性。

將該思想引用到連續(xù)體拓撲優(yōu)化中，針對本文變密度法模型，利用有限元法對設計區(qū)域進行離散化，每次迭代計算后對靈敏度場進行濾波，以濾波后的靈敏度場作為下次迭代的初始值，直至收斂。設j單元在某一次迭代后的靈敏度數(shù)值為，

現(xiàn)計算j單元濾波處理后的靈敏度值。在濾波范圍內(nèi)每個單元具有一個權重因子hi，hi的大小與該單元和當前計算單元的中心距離成反比，計算式如下：

式中，rf表示濾波特征半徑，disk(i, j )表示單元i和當前計算單元j的距離。當disk(i, j )＞rf時，hi取值恒為0。

則單元濾波后的靈敏度值可以表示為：

在本文所采用的拓撲優(yōu)化算法中，每一次迭代后，有一部分單元會被刪除，也就是說這一部分已被刪除的單元不應該再在濾波過程中起作用，為了忽略這部分單元的影響，將模型中加入懲罰系數(shù)p，則濾波函數(shù)可以修改為：

當單元存在時，pi取值為1，否則為零。利用這種方法在每一次迭代后對設計區(qū)域中所有設計單元進行濾波，形成濾波后的靈敏度場，作為下一次迭代的初始值，如此循環(huán)直到得到最優(yōu)拓撲形式。

從該算法的濾波半徑設置以及濾波公式可以看出，在濾波半徑范圍之內(nèi)的單元對單元靈敏度更新有影響，并且當濾波半徑趨近于0時，修改后的靈敏度值趨于不變，當濾波半徑趨近于無限大時，各單元靈敏度相等。通過對每次迭代后全局靈敏度進行修改，可以有效克服數(shù)值計算中的棋盤格現(xiàn)象和網(wǎng)格依賴。

2 數(shù)值模擬分析

添加了濾波算法后的基于SIMP的變秘密度法拓撲優(yōu)化方法的算法流程圖如圖1所示。

采用數(shù)值模擬方法加以驗證。分析模型采用單邊（左側(cè)）固支懸臂梁，右邊中點加載集中載荷。梁的彈性模量為210×109，泊松比為0.3，厚度為2mm，右側(cè)中點加載向下集中載荷為1000N，位移約束為-0.18×10-3m。有限元劃分網(wǎng)格24×16個。

圖1 加入濾波算法后的變密度拓撲優(yōu)化流程

圖2和圖3分別給出了未加入濾波算法的結構拓撲優(yōu)化結果和加入了濾波算法之后的新的拓撲優(yōu)化結果。表1是加入濾波算法前后的各項數(shù)據(jù)指標。

圖2 未加入濾波算法前的拓撲結構

圖3 加入濾波算法后的拓撲結構

從圖2到圖3的結構拓撲優(yōu)化結果可以看出，在加入濾波算法后，棋盤格現(xiàn)象消除；由表1可以看出，濾波后所得結構滿足約束要求，結構最終質(zhì)量略有變動但變動不大。棋盤格結構中存在的應力集中現(xiàn)象也得到了一定的控制，整體應力分布更加均勻。

表1 加入濾波算法前后的各項指標

3 結束語

基于工程中常用的SIMP變密度函數(shù)模型，進行了靈敏度分析，提出了以結構最小柔度為目標函數(shù)、位移為約束的變密度插值模型為基礎的濾波思想，并推導了具體的計算方法。通過算例分析，明顯消除了棋盤格現(xiàn)象，同時對網(wǎng)格依賴性有所控制。證明了這種濾波方法是切實有效的。

在濾波算法算例中，設置的濾波半徑為1.8，也就是說采取的是二階濾波，對于濾波算法來講，選取合適的濾波半徑也是非常重要的。濾波半徑過小，則靈敏度變化不大，不能起到很好的消除棋盤格式的效果；反之，則會導致靈敏度變化過于劇烈，所得結構完全脫離原來拓撲模型。同時，濾波法難以對全局數(shù)值變量進行控制，而且沒有嚴格的理論依據(jù)，有時會導致邊界擴散化，產(chǎn)生模糊的邊界。這也是濾波方法的不足之一，是在具體實施過程中需要密切注意的問題。

[1] 張麗, 饒華球, 昌俊康. 對連續(xù)體結構拓撲優(yōu)化的一點認識[J]. 制造業(yè)自動化. 2010, 32(2): 93-96.

[2] 左孔天. 連續(xù)體結構拓撲優(yōu)化理論與應用研究[D]. 武漢:華中科技大學, 2004.

[3] 郭中澤, 張衛(wèi)紅, 陳裕澤. 結構拓撲優(yōu)化設計綜述[J], 機械設計, 2007, 24(8): 1-4.

[4] Mlejnek H. P, Schirrmacher R..An engineer's approach to optimal material distribution and shape finding[J]. Comp.Meth. Appl. Mech. Engrg, 1993, 106: 1-26.

[5] M le jned H p, S chinm ascher R. An engineers approach to optimal material distribution and shape finding[J]. Comput Method Appl Mech Engrg 1993, 106(1-2): 1-26.

[6] 湯穎穎. 基于變密度法的連續(xù)體拓撲優(yōu)化設計[D], 西安:長安大學, 2008.

[7] 鄒春江, 左孔天, 向宇. 基于SIMP方法微電容加速度計結構固有頻率拓撲優(yōu)化[J], 科學技術與工程,2011,11(29): 7088-7091.

[8] 牟淑志, 杜春江. 基于ANSYS二次開發(fā)的結構拓撲優(yōu)化[J], 計算機應用與軟件, 2010, 27(2): 228-230.