朱杰

(西安衛(wèi)星測控中心，陜西 西安 710000）

0 引言

優(yōu)化設計，是指在多種設計方案之中選擇最優(yōu)方案的設計方法。它以最優(yōu)化理論為基礎，利用計算機做輔助設計，根據(jù)設計要求的性能目標從而確定目標函數(shù)，在滿足給定條件的各種約束下，尋找最優(yōu)設計方案。

二十世紀六十年代，Schmit[1]在多載荷條件下的彈性結(jié)構(gòu)優(yōu)化設計研究中，將有限元技術(shù)和數(shù)學規(guī)劃思想結(jié)合起來，形成了全新的結(jié)構(gòu)優(yōu)化思想和理論，對現(xiàn)代優(yōu)化設計起到了創(chuàng)造性的指導作用。尺寸優(yōu)化在這期間是結(jié)構(gòu)優(yōu)化設計的主要研究方向。七十年代，隨著優(yōu)化設計的不斷發(fā)展，從尺寸優(yōu)化到形狀優(yōu)化和拓撲優(yōu)化方向的發(fā)展成為了結(jié)構(gòu)優(yōu)化的研究方向。RDssow&Taylor[2]、Cheng&Olhoff[3]將材料微結(jié)構(gòu)概念引入連續(xù)體結(jié)構(gòu)變厚度板優(yōu)化的研究工作；Bendsoe和Kikuchi[4]在連續(xù)體結(jié)構(gòu)優(yōu)化中，將材料微結(jié)構(gòu)概念和均勻化方法相結(jié)合。在此之后的二十年里，拓撲優(yōu)化方法得到了迅速發(fā)展。Hassani和Hinton[5-6]等總結(jié)了拓撲優(yōu)化中均勻化方法和理論的運用。Bendsoe、Rozvany[7]等提出 SIMP（Solid Isotropic Material with Penalization）人工材料模型避免了均勻化方法復雜繁瑣的計算，進一步推動了拓撲優(yōu)化的發(fā)展。Sigmund和Petersson[8]等針對優(yōu)化中存在和棋盤格、網(wǎng)格依賴性、局部極值等數(shù)值不穩(wěn)定現(xiàn)象提出了多種解決策略。

拓撲優(yōu)化在材料設計、汽車工業(yè)和飛機制造等領域的應用得到了不錯的成果。一些大公司還專門內(nèi)部專用開發(fā)了自己公司使用的拓撲優(yōu)化軟件；部分有限元軟件，如 NASTRAN、ANSYS、Hyperworks、COMSOL也具有相應的拓撲優(yōu)化功能模塊，并且支持二次開發(fā)程序編寫。本文旨在結(jié)合學習到的關于最優(yōu)化的知識和自身研究方向相關的拓撲優(yōu)化相關理論，利用MATLAB編程實現(xiàn)對二維懸臂梁的拓撲結(jié)構(gòu)優(yōu)化設計，從而加深對知識的理解和運用。

1 變密度SIMP方法拓撲優(yōu)化基本原理

1.1 拓撲優(yōu)化思想

拓撲優(yōu)化旨在尋求某種結(jié)構(gòu)布局，使其性能指標能夠在滿足特定約束條件下達到最優(yōu)。即在一個給定的連續(xù)區(qū)域內(nèi)找到結(jié)構(gòu)內(nèi)部非實體區(qū)域最佳配置的位置和數(shù)量，尋求結(jié)構(gòu)中最優(yōu)化的節(jié)點聯(lián)結(jié)方式及構(gòu)建布局，使結(jié)構(gòu)能在滿足位移、應力等約束條件下，外載荷作用通過結(jié)構(gòu)傳遞到支撐位置，同時結(jié)構(gòu)達到最優(yōu)的某種性態(tài)指標[9]。

在對設計域內(nèi)的材料進行合理分配，達到最優(yōu)布局時，考慮到材料有限，所以最后的目的轉(zhuǎn)換為確定空間中哪些設計點為空，而哪些設計點有材料。考慮到計算效率問題，不可能將設計區(qū)域每一點都作為設計變量，采用有限元方法對設計區(qū)域進行離散，將離散后的單元作為設計變量。由于優(yōu)化的最終目的是確定哪些設計單元里面有材料，如果把每個單元最后的狀態(tài)描述為有和無，采用基本的枚舉法組合時，其組合方式為2n種，如果設計變量n較大時，將產(chǎn)生組合爆炸現(xiàn)象。同時對于離散后的優(yōu)化問題，不一定所有離散點都位于不等式約束條件所形成的曲面上，因此，計算過程中得到的可行點只能是接近于約束條件的解[10]。為了解決這個問題，需要在設計單元中加入中間變量，對整個問題進行松弛，擴大所允許的材料性能集合，將整數(shù)參量函數(shù)用一個新的連續(xù)分布的參量函數(shù)來代替，這樣原來的離散變量數(shù)學模型轉(zhuǎn)換為連續(xù)變量的數(shù)學模型，可以確?？尚薪獾拇嬖赱11-12]。

1.2 變密度SIMP材料插值模型方法基本原理

變密度法（Variable Density Method）是連續(xù)體結(jié)構(gòu)拓撲優(yōu)化的比較常用方法，它人為地假設材料的宏觀物理常數(shù)與其密度之間存在非線性關系，將連續(xù)體離散為有限元模型后，以每個單元的密度ρ(x)為設計變量，設定優(yōu)化目標后，在一定的約束條件下（一般為體積約束或質(zhì)量約束），通過一定的優(yōu)化算法更新迭代ρ(x)的值，最終得到優(yōu)化結(jié)果，其優(yōu)化數(shù)學模型為[13]：

其中ρi為每個單元的密度，C(ρ)為優(yōu)化目標方程，valume(m ass)constraint 為體積或質(zhì)量約束，equilibrium equation 為結(jié)構(gòu)平衡方程，ρmin為每個單元的密度下限。

固體各向同性材料懲罰函數(shù)模型SIMP（Solid Isotropic Material with Penalization）是目前拓撲優(yōu)化中最常用的密度插值模型方法之一。在SIMP材料插值模型中，有以下兩個前提和假設：

假設在離散單元內(nèi)部的材料是各向同性的，設計變量定義為離散單元的相對密度ρ(x)在[0-1]區(qū)間內(nèi)變化，當ρ(x)=0時代表此單元無材料，當ρ(x)=1時代表此單元為固體材料。

單元材料屬性隨著單元相對密度的變化而變化，并且是與單元相對密度ρ(x)成指數(shù)變化關系。在此條件上令材料的楊氏彈性模量E(x)隨ρ(x)的變化公式為：

并可推導出：

這里E0為固體材料的楊氏模量，ρ為懲罰系數(shù)，Ke為優(yōu)化中結(jié)構(gòu)總剛度矩陣，K0為固體材料的結(jié)構(gòu)總剛度矩陣，ke為優(yōu)化過程中的單元剛度矩陣，k0為固體材料的單元剛度矩陣。

SIMP材料插值模型的主要思想是引入懲罰因子ρ對中間密度值進行懲罰，使中間密度值向0/1兩端聚集收斂，避免中間密度單元的出現(xiàn)，這樣連續(xù)變量的拓撲優(yōu)化問題才能與實際中的0/1離散變量的整數(shù)規(guī)劃問題更接近，同時在減少中間密度單元后，整個結(jié)構(gòu)的材料分布也會更清晰，更利于工程設計中的參考使用[14-15]。

1.3 優(yōu)化準則算法

優(yōu)化準則算法依據(jù)一定的工程假設建立優(yōu)化準則和相應的迭代公式，然后利用數(shù)值法進行循環(huán)求解。對于求解單約束問題具有收斂速度快，求解結(jié)果穩(wěn)定的有點，但是對于不同的優(yōu)化問題需要重新推導優(yōu)化準則。以針對最小柔度為例，先根據(jù)目標函數(shù)和約束條件構(gòu)成拉格朗日函數(shù)

式中，λ,λ1, λ2,λ3為拉格朗日乘子，λ為標量，其他為向量。當ρ取極值ρ?時應該滿足的Kuhn-Tucker條件為：

因為外載荷不隨密度變化而變化：

當約束條件取上邊界時，即ρmin=ρ，設計變量的下限約束不起作用，即λ3=0。同理，當約束取下邊界時，即，設計變量的上約束不起作用，即λ2=0。當ρmin＜ρ＜ρmax時，上下約束都不起作用，即，此時可以將Kuhn-Tucker條件改寫為：

考慮ρmin＜ρ＜ρmax的情況，進一步地優(yōu)化準則為：

優(yōu)化問題的自變量迭代條件為：

m 為移動極限常數(shù)，為了避免在實際優(yōu)化過程中過多的刪除單元，取值一般在0.1-0.3之間；阻尼因子η同樣是為了優(yōu)化的穩(wěn)定性，取值在0.3-0.5之間。

1.4 敏度過濾

采用有限元法進行拓撲優(yōu)化的數(shù)值計算過程中，由于設計區(qū)域內(nèi)的離散變量存在奇異性問題，會出現(xiàn)數(shù)值計算不穩(wěn)定的情況，通常導致的結(jié)果就是得到的最終解是非全局最優(yōu)解而是局部最優(yōu)解。在日常的優(yōu)化設計中，此類現(xiàn)象的出現(xiàn)經(jīng)常會導致優(yōu)化結(jié)果不可用。在結(jié)果中有時會出現(xiàn)如棋盤格現(xiàn)象、網(wǎng)格依賴性等數(shù)值不穩(wěn)定現(xiàn)象。拓撲優(yōu)化中數(shù)值計算不穩(wěn)定性的消除非常重要，它關系到數(shù)值計算的收斂性和計算結(jié)果的可制造性問題[16]。

為了解決這類數(shù)值不穩(wěn)定現(xiàn)象，引入敏度過濾：

式中dist(e,f)表示單元e中心和單元f中心之間的距離。

2 二維懸臂梁算例的MATLAB實現(xiàn)

2.1 問題描述

如圖3-1所示懸臂梁，左端固支，假定材料為各向同性，為建模和計算簡單化，采用無量綱模式計算，長寬比為3∶1，彈性模量為1，右下部受力大小為1，泊松比為0.3。先要求設計體積為原體積的0.5倍，用變密度SIMP方法進行拓撲優(yōu)化，實現(xiàn)最大剛度。

圖3 -1 懸臂梁結(jié)構(gòu)示意圖Fig. 3-1 Schematic diagram of cantilever beam structure

2.2 優(yōu)化過程

這是一道典型的簡單的資源約束條件下的最大剛度設計。故先建立數(shù)學模型如下：

用MATLAB程序?qū)崿F(xiàn)優(yōu)化過程的流程圖如下圖3-2，全部MATLAB程序在附錄中可以查看

圖3 -2 計算流程圖Fig.3-2 Calculationflowchart

2.3 優(yōu)化結(jié)果

在MATLAB中輸入下面語句：

＞＞topmain（60,20,0.5,3,1.5）

調(diào)用主程序程序，即給定參數(shù)為橫坐標劃分網(wǎng)格數(shù)60，縱坐標20，體積分數(shù)0.5，懲罰系數(shù)3，過濾半徑1.5。程序運行迭代總次數(shù)為94次，分別取第5、10、15、20、25、45、75、94次迭代時的過程圖像來展示迭代的優(yōu)化過程，如圖3-3（a）-（h）。

圖3 -3 迭代過程中的拓撲結(jié)構(gòu)圖Fig. 3-3 Topology diagram in the iterative process

從上面的優(yōu)化過程圖中可以看出，在迭代次數(shù)到了25次以后，基本已經(jīng)形成了最后迭代終止時的拓撲結(jié)構(gòu)，查看迭代過程中的變量，發(fā)現(xiàn)25次之前目標函數(shù)值減小的很快，第1、5、10、15、20、25次迭代時目標函數(shù)值分別為 1007.0221、322.0344、272.6553、232.2493、212.4871、206.2165，而到了25次之后收斂速度開始減慢，第45、75、94次迭代對應的目標函數(shù)值分別為 204.2806、203.4899、203.3061。

3 結(jié)論

文章利用MATLAB軟件結(jié)合學習到的關于最優(yōu)化的知識和自身研究方向相關的拓撲優(yōu)化相關理論，采用變密度SIMP材料插值模型和優(yōu)化準則迭代方法，以整體柔度最小，即剛度最大為目標函數(shù)，實體材料所占體積為設計區(qū)域50%為約束條件，對二維懸臂梁的拓撲結(jié)構(gòu)進行了優(yōu)化設計[17]。優(yōu)化結(jié)果表明，拓撲優(yōu)化能夠有效的在滿足工況條件下減少對材料的使用，且因為加入了有較成熟理論基礎的敏度過濾步驟，所以優(yōu)化后的結(jié)果并沒有出現(xiàn)棋盤格式或者網(wǎng)格依賴性現(xiàn)象。