印亞榮　傘冰冰

(河海大學 土木與交通學院， 南京　210098)

?

基于擴展有限元法的連續(xù)體結構拓撲優(yōu)化

印亞榮傘冰冰

(河海大學 土木與交通學院， 南京210098)

摘要：利用擴展有限元法能夠在結構內部出現(xiàn)缺陷時無需重新劃分網格、簡化有限元分析計算的特點，將其與拓撲優(yōu)化相結合，計算在變密度法的SIMP (Solid isotropic microstructures with penalization)模型下的連續(xù)體結構拓撲優(yōu)化．建立結構在體積約束下的結構拓撲優(yōu)化模型，將其與結合普通有限元的拓撲優(yōu)化進行對比分析，對普通結構分析比較結果的一致性表明其對于結構拓撲優(yōu)化問題的可用性．對有孔洞約束下的平面結構和殼體結構的分析結果表明其對于有缺陷結構拓撲優(yōu)化問題的網格劃分更加簡單，最終拓撲圖形不會產生由孔洞約束而產生的尖端和應力不均勻現(xiàn)象．

關鍵詞：拓撲優(yōu)化；擴展有限元；變密度法

結構拓撲優(yōu)化是結構優(yōu)化領域研究的熱點問題，從研究Michael桁架的離散體結構優(yōu)化到Bendson和Kikuchi提出連續(xù)體拓撲優(yōu)化的均勻化方法[1]以來，連續(xù)體結構拓撲優(yōu)化的方法不斷出現(xiàn)，如漸近結構優(yōu)化法[2](ESO)、變密度法[3-4]、變厚度法等．Sigmund等對基于變密度法的建模方法進行深入研究，提出了正交各向同性微結構材料懲罰模型法[5](SIMP)，G. L. N. Rozany詳細闡述了該模型的主要優(yōu)勢和背景資料．

目前，變密度法成為應用最廣泛的拓撲優(yōu)化方法之一，廣泛應用于多種約束和響應的連續(xù)體拓撲優(yōu)化問題．因其理論難度較低，數學模型容易實現(xiàn)、計算效率高，現(xiàn)已被廣泛應用．但是Diaz和Sigmund[7]指出，拓撲優(yōu)化采用有限元方法對設計區(qū)域進行離散化時，有時網格的不連續(xù)的排列反而使材料具有更好的“虛擬”剛度．這就造成了數值不穩(wěn)定現(xiàn)象的發(fā)生，例如棋盤格式、網格依賴性等問題．

美國西北大學Belytschko教授等人首先提出擴展有限元思想[8]，并將其應用到拓撲優(yōu)化中[9]，求解結構中有多種材料存在的優(yōu)化問題．擴展有限元法是以有限元為基本框架，以單位分解為基礎，主要針對不連續(xù)問題進行研究．其網格劃分與結構內部的幾何、物理界面無關，克服了在高應力和變形集中區(qū)進行高密度網格劃分所帶來的困難，從而可以克服由于網格問題引起的數值不穩(wěn)定現(xiàn)象．因此本文利用其擴展有限元法對網格劃分的優(yōu)勢，與拓撲優(yōu)化相結合來解決結構中出現(xiàn)的數值不穩(wěn)定現(xiàn)象．

由于一般結構的數值不穩(wěn)定現(xiàn)象發(fā)生的不確定性，對于內部有固定孔洞的結構，內部的形式是由孔洞決定的，利用有限元法時肯定會出現(xiàn)局部網格劃分不均勻，造成計算量的增大和優(yōu)化后應力集中的情況．因此本文分析內部具有孔洞的結構，利用ABAQUS軟件進行分析優(yōu)化，對平面結構不直接在結構開孔洞而是在孔洞外側建立擴展有限元性質的缺陷繼而進行拓撲優(yōu)化，結果證明可以讓結構在不提高計算量和無需二次網格過濾的基礎上，得到同樣具有孔洞且無應力集中的結構．

1SIMP拓撲優(yōu)化理論

變密度法是從均勻化方法發(fā)展來的，它把拓撲優(yōu)化問題轉化為材料的最優(yōu)分布問題，其基本思想是人為假定一種實際工程中并不存在的密度可變材料單元，從而材料的物理屬性可以以材料單元密度函數的形式表達．以區(qū)間[0，1]內的密度值為設計變量，直接定義一個經驗公式來表達密度與彈性模量間假定的函數關系．繼而引入懲罰因子，對設計變量在(0，1)之間的中間密度值進行懲罰，使中間密度值逐漸向0/1兩端聚集．這時除了密度靠近1的單元保留下來，其他單元對結構剛度矩陣影響將變得很小，可以忽略不計．由于變密度法設計變量的數量較少，容易通過編程實現(xiàn)，計算效率也比較高，因此得到了廣泛的應用．

固體各向同性材料懲罰模型是一種常用的密度-剛度插值模型，從對中間密度的懲罰效果來看，SIMP插值模型的處理效果比RAMP插值模型略好．SIGMUND等對變密度法材料插值模型進行了深入研究[10]，從理論上研究了各種不同變密度法的材料插值模型，證明了中間密度物理意義的存在．

1.1優(yōu)化的數學模型

基于變密度理論的拓撲優(yōu)化是以單元相對密度為設計變量，拓撲優(yōu)化數學模型為：求x={x1，x2，…，xN}T

(1)

圖1　SIMP插值模型

p為懲罰系數，p>1會降低中間密度材料的產生，p越大懲罰效果越強，但p太大會引起結構剛度矩陣的病態(tài)，以及迭代過度懲罰得不到合理的結果，這里取p=3．

1.2優(yōu)化算法

優(yōu)化模型為帶不等式約束的非線性規(guī)劃問題，特點為設計變量多，每個迭代步都要進行結構分析，求解拓撲優(yōu)化的主流算法可以分為數學規(guī)劃法[11]、啟發(fā)式算法和優(yōu)化準則法[12]3類．

優(yōu)化準則法是把最優(yōu)解需要滿足的Kuhn-Tucker條件[13]作為最優(yōu)結構應滿足的準則，通過引入Lagrange函數，引入Lagrange乘子，將有約束的非線性優(yōu)化問題轉化為無約束的優(yōu)化問題．優(yōu)化準則法下的變密度法拓撲優(yōu)化就是通過改變密度使得結構目標函數的改變和約束條件的變化相平衡，逐步迭代使得結構滿足Kuhn-Tucker所要求的駐值條件．優(yōu)化準則法要求重分析的次數一般跟變量的數目沒有多大關系，收斂速度快，迭代次數少，最主要是和結構的復雜程度和大小無關，所以對中型和大型結構的優(yōu)化設計有重要的實際意義．

2擴展有限元法

傳統(tǒng)有限元對于位移場的描述是基于單元的，每個單元內部的位移場函數總是通過形函數和單元結點位移來表達的：

(2)

擴展有限元法是單位分解方法[14-15]的特例，利用有限元形函數作為單位分解函數，引入非連續(xù)位移模式來描述非連續(xù)性位移場．擴展有限元法將模型分為：(Ⅰ)在忽略內邊界情況下對區(qū)域進行網格剖分；(Ⅱ)在單元形狀函數中增加與內邊界有關的附加函數，改進有限元逼近空間，既得擴展有限元法下全域的位移逼近函數為：

(3)

式中，Ni、Nj為形狀函數矩陣，ψj為改進函數．ui表示常規(guī)自由度；Ne表示單元改進結點數；Nn表示不存在不連續(xù)界面的單元結點數；aj表示單元改進結點上的附加自由度．

對裂紋問題所涉及的兩個改進函數描述如下：廣義Heaviside函數H(x)：對于裂紋表面采用Heaviside函數H(x)，在裂紋上方H(x)取1，裂紋下方H(x)取-1；裂尖函數ψa(x)：為了模擬裂紋尖端，改善斷裂計算中裂尖場的精度，各種問題下裂尖改進函數不相同．

擴展有限元法在模擬裂紋擴展問題的應用最為成熟，可以通過在含有裂紋的結點處增加額外自由度來反映裂紋的存在．相對有限元方法來說，擴展有限元法顯著提高了描述復雜位移場的能力，對于處理靈活的非連續(xù)邊界問題，避免了網格重新劃分的工作．對于拓撲優(yōu)化方法中，利用擴展有限元法復雜內邊界的網格劃分問題，其優(yōu)勢是去除了可能發(fā)生的拓撲圖形凹凸或尖端的問題，省去了后期對不穩(wěn)定現(xiàn)象的二次優(yōu)化，較好地提高了拓撲優(yōu)化的效率．

3算例分析

算例1：如圖2所示的0.16m×0.1m×0.006m平面矩形懸臂梁結構，材料彈性模量為10GPa，泊松比為0.3，密度為10t/m3，右邊界中點施加一集中荷載F=9kN．采用SIMP方法，以結構的最小化應變能為目標函數，體積比小于40%為約束條件，結合擴展有限元法的拓撲優(yōu)化分析中．以最大主應力失效準則作為擴展起始的依據，最大主應力為300MPa，擴展演化選取基于能量的、線性軟化的、混合模式的指數演化規(guī)律，粘性系數為0.001．利用ATOM模塊計算在減少的體積下，材料的合理排布可以得到的最優(yōu)剛度解，即得到最優(yōu)拓撲結構．本例經過51次迭代結構重分析和優(yōu)化求解后達到收斂狀態(tài)，圖3為拓撲優(yōu)化結果的mises應力云圖．

圖2　拓撲模型　　圖3　拓撲結果的Mises應力云圖

拓撲優(yōu)化中應變能的迭代變化曲線和體積比迭代變化曲線如圖4和圖5所示，該拓撲優(yōu)化過程是在體積比小于40%的基礎上將應變能優(yōu)化到最小值的過程，結構首先滿足約束條件體積比的要求，進而進行應變能最小化的進程．表1列出了兩次優(yōu)化結果的應變能、體積比以及最大mises應力的數值結果．

圖4　應變能迭代變化曲線　　圖5　體積比迭代變化曲線

參數應變能/J最大mises應力/MPaFEM33.4136344.722XFEM33.4089394.428

由于變密度法計算的體積比是結構中有效單元的體積比，與實際刪減單元后的體積比不同，所以與實際拓撲迭代過程體積變化不一致．在ATOM模塊中進行拓撲優(yōu)化的計算，結果顯示其最優(yōu)拓撲形式和運用普通有限元法計算的最優(yōu)拓撲形式幾乎相同，一方面說明了結合擴展有限元法的拓撲優(yōu)化的計算結果正確可用，另一方面表明其為未來結合擴展有限元的拓撲優(yōu)化進一步的研究提供了一個分析的基礎空間．

算例2：由于在實際工程結構中的拓撲優(yōu)化有額外的實際要求，比如有某一些構件需預留孔洞以滿足特定位置的構造要求，就需要結構拓撲優(yōu)化后的特定位置留有孔洞，對于普通有限元方法而言，可以在算例1結構開孔洞后再進行拓撲優(yōu)化．利用擴展有限元法對于處理缺陷情況的優(yōu)勢，在其孔洞兩側建立擴展有限元性質的缺陷，由于所需要的孔洞包含在建立的缺陷內部，因此缺陷經過優(yōu)化過程中的擴展就可以形成所需要的孔洞，如圖6是擴展有限元方法下的結構的網格劃分情況和最終拓撲圖形．該結構經過33次的迭代過程達到收斂狀態(tài)，圖7顯示出了其拓撲優(yōu)化迭代過程圖形．

圖6　XFEM方法下的起始和最終拓撲圖形

圖7　XFEM拓撲優(yōu)化迭代過程

基于普通有限元的拓撲優(yōu)化經過38次迭代達到收斂，分析結果圖形和應變能變化曲線如圖8和圖9所示．

圖8　拓撲優(yōu)化起始和最終圖形

圖9　應變能變化曲線

由于孔洞對于結構網格的劃分的影響，使得拓撲優(yōu)化結果的圖形在預開洞口的位置會出現(xiàn)尖端，為了使結果更加合理就必須利用網格過濾法對其進行二次優(yōu)化．基于普通有限元方法的拓撲優(yōu)化結果顯示，在最后變形圖中會出現(xiàn)結構的1桿件結構較細、應力值較小的情況，造成了最終拓撲圖形不符合所要求的結構形式，不利于在實際工程中的應用．

算例3：如圖10所示為雙曲扁殼模型，尺寸為1 m×1 m，厚度0.01 m，四角固定，材料的彈性模量為E0=2×1011Pa，泊松比為0.3，集中力F=100 N．結構在中間需開有半徑為100 mm的圓形孔洞，網格劃分如圖11所示．建立SIMP模型，以最小化應變能為目標函數，體積比小于0.4作為約束條件，選取懲罰因子為p=3．

圖10　雙曲扁曲殼模型　　　　圖11　結構網格劃分

基于擴展有限元方法下的拓撲優(yōu)化過程中，經過45次迭代過程達到收斂狀態(tài)．普通有限元方法下的拓撲優(yōu)化經過55次迭代達到收斂狀態(tài)，其拓撲優(yōu)化結果和應變能的變化曲線如圖12和13所示．

圖12　拓撲優(yōu)化最終圖形

圖13　應變能變化曲線

兩者結果做對比可以看出最終結構在應變能相差不大的情況下，擴展有限元方法收斂速度更快，結果更加簡單，其最終結果在滿足孔洞要求的同時，其內部邊界又更利于操作加工．由于結構受力對稱，所以最終優(yōu)化結構在上下方向處于對稱狀態(tài)，但是普通有限元方法下的結構由于單元網格的劃分不能與擴展有限元方法下和缺陷的擴展無關，影響了后期孔洞的形成，在左右方向上并不能和想象的一樣處于對稱狀態(tài)．

4結論

1)本文通過將擴展有限元法與拓撲優(yōu)化方法的結合，利用擴展有限元方法對于網格劃分處理的優(yōu)點，結合SIMP變密度法，對同一模型分別進行普通有限元和擴展有限元方法下的拓撲優(yōu)化分析，結果表明了基于擴展有限元法的拓撲優(yōu)化對于結構的優(yōu)化合理且正確．

2)本文對于需在特定位置留構造孔洞的平面結構和殼體結構分別進行計算比較，結果顯示其結果更加合理，在實際中具有較好的實用性．對于平面結構和殼體結構將結構內邊界脫離了孔洞的約束，在滿足空洞要求的前提下結構的應變能與普通有限元方法下的結果一致，也避免了由于孔洞而可能會出現(xiàn)的結構形式問題和應力問題．

3)本文的計算結果從整體模型來看與前者差距不是很大，但是從應用角度而言其優(yōu)越性就能體現(xiàn)出來，當計算到體量較大、模型較為復雜的時候，計算效率的優(yōu)勢以及模型局部合理性的優(yōu)勢就會體現(xiàn)出來，更好地為實際工程服務．

參考文獻：

[1]Bendsoe M P, Kikuchi N. A Homogenization Method for Shape and Topology Optimization[J]． Computer Methods in Applied Mechanics and Engineering, 1991, 93: 291-381．

[2]Xie Y M, Steven G P. A Simple Evolutionary Procedure for Structural Pptimization[J]． Computers & Structures, 1993, 49(5): 885-896．

[3]Sigmund O. Design of Material Structures Using Topology Optimization: [PhD thesis][M]. Denmark: Department of Solid Mechanics, Technical University of Denmark, 1994．

[4]Sigmund O, Bendsoe M P. Material Interpolations in Topology Optimization[J]． Archive of Applied Mechanics, 1999, 69: 635-654．

[5]Rietz A. Sufficiency of a Finite Exponent in SIMP (Power Law) Methods[J]． Structural and Multidiscipline Optimization, 2001, 21: 159-163．

[6]Stolpe M, Svanberg K. An Alternative Interpolation Scheme for Minimum Compliance Topology Optimization[J]． Structural and Multidiscipline Optimization, 2001, 22: 116-124．

[7]Diaz A R, Sigmund O. Checkerboard Patterns in Layout Optimization[J]． Structural and Multidisciplinary Optimization, 1995, 10: 40-45．

[8]Moes N, Dolbow J, Belytschko T. A Finite Element Method for Crack Growth Without Remeshing[J]． International Journal for Numerical Method in Engineering, 1999, 46: 131-150．

[9]Belytschko T, Xiao S P, Parimi C. Topology Optimization with Implicit Functions and Regularization[J]． International Journal For Numerical Methods In Engineering, 57:1177 - 1196, 2003．

[10] Bendsoe M P, Sigmund O. Material Interpolation Schemes in Topology Optimization[J]． Archive of Applied Mechanics, 1999, 69: 635-654．

[11] Sehmit L A, Mallett R H. Structural Synthesis and Development[J]． Journal of Structural Engineering-American Society of Civil Engineers, 1963, 89(3): 269-299．

[12] Zhou M. An Efficient DCOC Algorithm Based on High-quality Approximation for Problems Including Eigenvalue Constraints[J]． Computer Methods in Applied Mechanics and Engineering, 1995, 128: 3837-394．

[13] Rozvany G I N, Zhou M. The COC Algorithm, part I: Cross-section Optimization Sizing[J]． Computer Methods in Applied Mechanics and Engineering, 1991, 89: 281-308．

[14] Melenk J M, Babuska I. The Partition of Unity Finite Element Method: Basic Theory and Applications [J]． Computer Methods in Applied Mechanics and Engineering, 1996, 139: 289-314．

[15] 李錄賢，王鐵軍. 擴展有限元法(XFEM)及其應用[J]．力學進展， 2005(1)：5-20．

[責任編輯周文凱]

收稿日期：2015-10-20

基金項目：國家自然科學基金資助(項目批準號：51578211)

通信作者：印亞榮(1991-)，男，碩士研究生，研究方向為拓撲優(yōu)化．E-mail：yinyarong@163.com

DOI：10.13393/j.cnki.issn.1672-948X.2016.01.012

中圖分類號：TU311.41

文獻標識碼：A

文章編號：1672-948X(2016)01-0057-05

Structural Topology Optimization of Continuum Using Extended Finite Element Method

Yin YarongSan Bingbing

(College of Civil & Transportation Engineering, Hohai Univ., Nanjing 210098, China)

AbstractThe extended finite element method can solve the problem when the mesh is splited by crack without remeshing, which can simplify the finite element analysis. We analyze the structural topology optimization of continuum using the Solid isotropic microstructures with penalization, taking full advantage of extended finite element method. The model under the constraint of volume was established using alterable density method for taking the topology optimization with generalized finite element method and that with extended finite element method into comparative analysis. The result supports that the topology optimization using extended finite element method is available and efficient. Then we analyze the original planar and shell structure constraint with a fixed hole, which show topology optimization with extended finite element method can simplify the meshing and eliminate the possibility of producing the tip or inhomogeneous stresses.

Keywordstopology optimization；extended finite element method；alterable density method