基于MATLAB和ANSYS的懸臂梁拓撲優(yōu)化

1、現(xiàn)代設(shè)計理論和方法大作業(yè)-基于MATLA和ANSYS勺懸臂梁拓撲優(yōu)化指導(dǎo)老師：志岡H 項目組長：薛亞波 項目成員：機自 66 學(xué)院：機械工程學(xué)院基于MATLA和ANSYS勺懸臂梁拓撲優(yōu)化、計劃和任務(wù)安排表1.1學(xué)習(xí)時間安排表:任務(wù)時間查閱整理資料5月03日一5月08日學(xué)習(xí)相關(guān)知識和論文5月09日一5月15日設(shè)計實驗方案5月16日一5月17日學(xué)習(xí)優(yōu)化理論5月18日一5月20日撰寫論文5月21日一5月22日集體交流學(xué)習(xí)心得5月23日一5月24日1.2任務(wù)分配表學(xué)號任務(wù)06011136齊寰宇查閱資料、答辯PPT制作06011137王曉慶查閱資料、答辯PPT制作06011139冰查閱資料、研究問題分析

2、06011141康辰龍查閱資料、ansys拓撲模 塊研究06011143梁偉云查閱資料、項目背景研究06011145盧崗查閱資料ansys拓撲模塊 研究06011153清偉查閱資料ansys拓撲模塊 研究06011154薛亞波查閱資料、后期工作總結(jié)06011155嚴宏偉查閱資料ansys拓撲模塊 研究06011156志亮查閱資料、項目背景研究06011157龍查閱資料、答辯PPT制作06012020喬百杰查閱資料、Matlab結(jié)構(gòu)優(yōu) 化06073012鄧磊查閱資料、項目背景研究06182047立峰查閱資料、研究問題分析06055164杜博查閱資料、研究問題分析06011146牛玉輝查閱資料、研

3、究問題分析06011147樸群星查閱資料、項目背景研究、項目背景介紹及問題描述2.1項目背景及意義：2.1.1工程背景及基本原理：通常把結(jié)構(gòu)優(yōu)化按設(shè)計變量的類型劃分成三 個層次：結(jié)構(gòu)尺寸優(yōu)化、形狀優(yōu)化和拓撲優(yōu)化。尺寸優(yōu)化和形狀優(yōu)化已得到充分 的發(fā)展，但它們存在著不能變更拓撲結(jié)構(gòu)的缺陷。在這樣的背景下，人們開始研究拓撲優(yōu)化。拓撲優(yōu)化的基本思想是將尋求結(jié)構(gòu)的最優(yōu)拓撲問題轉(zhuǎn)化為在給定的 設(shè)計區(qū)域?qū)で笞顑?yōu)材料的分布問題。尋求一個最佳的拓撲結(jié)構(gòu)形式有兩種基本的 原理：一種是退化原理，另一種是進化原理。退化原理的基本思想是在優(yōu)化前將 結(jié)構(gòu)所有可能桿單元或所有材料都加上，然后構(gòu)造適當?shù)膬?yōu)化模型，通過一定的優(yōu)

4、化方法逐步刪減那些不必要的結(jié)構(gòu)元素，直至最終得到一個最優(yōu)化的拓撲結(jié)構(gòu) 形式。進化原理的基本思想是把適者生存的生物進化論思想引入結(jié)構(gòu)拓撲優(yōu)化， 它通過模擬適者生存、物競天擇、優(yōu)勝劣汰等自然機理來獲得最優(yōu)的拓撲結(jié)構(gòu)。2.1.2本文研究意義：目前，結(jié)構(gòu)優(yōu)化大部分集中在尺寸設(shè)計變量 （如板 厚、桿的剖面積及管梁的直徑）。拓撲結(jié)構(gòu)優(yōu)化較尺寸優(yōu)化復(fù)雜，但對于有些問 題拓撲結(jié)構(gòu)優(yōu)化比尺寸優(yōu)化有效，懸臂梁是其中的例子之一。本文討論懸臂梁的 拓撲優(yōu)化問題，圍繞這一問題，怎樣使結(jié)構(gòu)具有最大剛度的設(shè)計占有相當重要的 地位；怎樣優(yōu)化結(jié)構(gòu)的形狀使材料的分布，更加合理從而達到使結(jié)構(gòu)具有最大剛 度的目的是本文要研究的問題。

5、2.2研究現(xiàn)狀2.2.1理論研究現(xiàn)狀：結(jié)構(gòu)拓撲優(yōu)化是近20年來從結(jié)構(gòu)優(yōu)化研究中派生出來 的新分支，它在計算結(jié)構(gòu)力學(xué)中已經(jīng)被認為是最富挑戰(zhàn)性的一類研究工作。目前有關(guān)結(jié)構(gòu)拓撲優(yōu)化的工程應(yīng)用研究還很不成熟，在國外處在發(fā)展的初期，尤其在 國尚屬于起步階段。1904年Michell在桁架理論中首次提出了拓撲優(yōu)化的概念。 自1964年 Dorn等人提出基結(jié)構(gòu)法，將數(shù)值方法引入拓撲優(yōu)化領(lǐng)域，拓撲優(yōu)化研究 開始活躍。20世紀80年代初，程耿東和N.OIhoff在彈性板的最優(yōu)厚度分布研究中 首次將最優(yōu)拓撲問題轉(zhuǎn)化為尺寸優(yōu)化問題，他們開創(chuàng)性的工作引起了眾多學(xué)者的 研究興趣。1988年Bendsoe和Kikuchi

6、發(fā)表的基于均勻化理論的結(jié)構(gòu)拓撲優(yōu)化設(shè) 計，開創(chuàng)了連續(xù)體結(jié)構(gòu)拓撲優(yōu)化設(shè)計研究的新局面。1993年XieYM和Steve nGP提出了漸進結(jié)構(gòu)優(yōu)化法。1999年 Ben dsoe和Sigm unc證實了變密度法物理意義的存在 性。2002年羅鷹等提出三角網(wǎng)格進化法，該方法在優(yōu)化過程中實現(xiàn)了退化和進化 的統(tǒng)一，提高了優(yōu)化效率。2.2.2應(yīng)用研究現(xiàn)狀：在前人提出的重要理論基礎(chǔ)上，后人也將其跟其他現(xiàn)代設(shè) 計的方法相結(jié)合，衍生出了其他一些拓撲結(jié)構(gòu)優(yōu)化方法：如與可靠性相結(jié)合的情況下，MAUT等應(yīng)用變密度法并結(jié)合可靠性分析對一微機電系統(tǒng)進行了基于可靠 性的拓撲優(yōu)化設(shè)計，PAPADRAKAKIS將遺傳算法應(yīng)用于

7、具有可靠性約束的桁架結(jié) 構(gòu)拓撲優(yōu)化設(shè)計中，國學(xué)者馬洪波也對基于遺傳算法的結(jié)構(gòu)可靠性優(yōu)化問題進行 了討論。華南理工大學(xué)機械工程學(xué)院歐陽高飛等對基于水平集方法的結(jié)構(gòu)可靠性 拓撲優(yōu)化進行了研究。2.3研究目標：2.3.1設(shè)計目標：設(shè)計一懸臂梁在自由端受到豎直向下集中力P的作用，要求保持懸臂梁長度尺寸不變。2.3.2優(yōu)化目標：在滿足以上條件情況下對原實體懸臂梁做結(jié)構(gòu)拓撲優(yōu)化設(shè) 計，使其結(jié)構(gòu)剛度最大、重量最輕。（優(yōu)化的結(jié)果應(yīng)該使原設(shè)計區(qū)域產(chǎn)生孔洞， 使結(jié)構(gòu)拓撲發(fā)生變化。）優(yōu)化目標分析：要使梁滿足以上條件情況必須同時兼顧重量和 剛度問題。重量太輕不能滿足剛度問題。剛度太大又必須足夠的材 料以滿足其要求。將

8、兩個因素進行耦合分析，我們可以得出的結(jié)論 是必須通過結(jié)構(gòu)層次來進行優(yōu)化。連續(xù)體結(jié)構(gòu)拓撲優(yōu)化較成熟的方 法主要有均勻化方法、變密度方法和漸進結(jié)構(gòu)優(yōu)化法，通過綜合比 較我們最終選用變密度法。2.3.3實現(xiàn)方法：本文在對題目進行分析后，決定從兩個獨立的方向來分別 研究，以此來互相驗證，保證結(jié)果的可靠性。實驗中要使用ansys和matlab兩種軟件分別獨立進行優(yōu)化分析，然后將兩種結(jié)果來進行比較。2.3.4實踐目標：通過這次的作業(yè)我們期望通過搜索資料，團體討論，分組作業(yè)的方式，以期望達到對拓撲優(yōu)化的方法和過程有個具體的了解和學(xué)習(xí)。另外對分析軟件的應(yīng)用能夠達到一個新的高度。這些不僅能使我們現(xiàn)在的知識體系得

9、 到充實和優(yōu)化，而且也是我們今后人生的財富2.4研究容2.4.1問題描述如圖所示，懸臂梁在自由端受到豎直向下集中力 P的作用，要求在保持懸臂 梁長度尺寸不變的情況下對原實體懸臂梁做結(jié)構(gòu)拓撲優(yōu)化設(shè)計，優(yōu)化目標是使結(jié)構(gòu)剛度最大、重量最輕。（優(yōu)化的結(jié)果應(yīng)該使原設(shè)計區(qū)域產(chǎn)生孔洞，使結(jié)構(gòu)拓撲 發(fā)生變化。）原實體懸臂梁橫截面為矩形，矩形寬度為 0.2m,懸臂梁的材料為45鋼，密 度p為7.8 x 10-6kg/mm3,彈性模量E=2X 105MPa許用切應(yīng)力t =60MPa豎 直向下集中力為P 1.5 104n。圖 懸臂梁受集中力作用(單位： m242問題轉(zhuǎn)化本研究中，要求在保持懸臂梁在結(jié)構(gòu)尺寸不變的情況

10、下， 使懸臂梁的剛度最 大、質(zhì)量最輕。一般懸臂梁都為三維立體結(jié)構(gòu)，載荷是在豎直方向上加載的，所以在垂直紙 面方向上，并且在懸臂梁的對稱平面上加載，根據(jù)材料力學(xué)知識可知，所加載在 和只對加載區(qū)周圍較小圍產(chǎn)生影響， 其余部分與加載形式無關(guān)。其受力與如圖所 示平面的手里相同，故可以將問題抽象為平面問題來處理。當占有區(qū)域S的一個結(jié)構(gòu)上作用有強度為P ( S )的載荷時,使結(jié)構(gòu)具有最大 剛度的問題和使結(jié)構(gòu)所受外力功 W具有最小值的問題是等價的。結(jié)構(gòu)所受的外力 功W與結(jié)構(gòu)變形能C之間的關(guān)系為W = 2C,所以我們可以將該問題處理為求在一 定的約束條件下，該懸臂梁的的最小變形能最小的問題。243數(shù)學(xué)模型：該

11、問題中，要求同時滿足剛度最大，質(zhì)量最輕，這兩個變量若同時改變， 則問題復(fù)雜度太大，并親切可能導(dǎo)致問題不可求解。所以我們采用在確定的質(zhì)量 下，來討論剛度最大的問題。由于對特定的材料，其質(zhì)量和體積有一定的關(guān)系， 并且我們采用去除法的思想來建立模型的， 故我們可以采用給優(yōu)化后的體積與優(yōu) 化前的體積比賦確定的值，來達到在給定質(zhì)量條件下滿足剛度最大的問題。 其數(shù) 學(xué)模型如下：Nmin : c(x) = UtKU 二uj k0 ueC=1subject to: f:KU= F:0 Xm 迫 xl(1)注：其中C (x)為結(jié)構(gòu)變形能，U為結(jié)構(gòu)變形總位移矩陣，K為結(jié)構(gòu) 總剛度矩陣，N為劃分單元總數(shù)，Ue為單元

12、位移向量，ko為單元剛度， (由于劃分單元的時候，我們采用等分舉行單元，所以個單元的剛度可用 一個常量來處理)V(x)是拓撲結(jié)構(gòu)優(yōu)化過程中變化著的體積，Vo為未經(jīng)過優(yōu)化前懸臂梁的體積。F為結(jié)構(gòu)所受的總載荷。x為懸臂梁的相對密度。243模型分析求解：該問題的優(yōu)化方法有很多種，常用的有如下方法：Optimality Criteria(OC) methods, (優(yōu)化準則方法)Seque ntial Linear Programmi ng (SLP) methods(序列線性規(guī)劃法)Method of Movi ng Asymptotes (MMA bySva nberg 1987) 等 為了簡化問

13、題的復(fù)雜度，此處我們采用 standard OC-method.方法來實 現(xiàn)。在處理過程中，關(guān)于設(shè)計變量相對密度x每一步的更新，我們采用在1995 年提出的如下算法來實現(xiàn)：礙榨=:max(厲金亡m):if %石2 一 m),I if m) xeB niin( 1, re + m),:+ m)i if mm(l5 xe + m T注：其中m是移動限制量，即相對密度變化的最大步長，更新條件中，加上最小邊界條件判斷，是防止相對密度變?yōu)?,使問題出 現(xiàn)奇異；加上相對密度最大邊界條件判斷，是為了防止出現(xiàn)相對密度大于1,是問題是去物理意義。是數(shù)值阻尼系數(shù)。Be可以由以下條件來確定:其中 為拉格朗日乘數(shù)，由

14、雙向切片算法來確定目標函數(shù)對單元相對密度的變化率，即單元敏感度如下確定：deb丁二-PMP uf k0 % (4)為了確保該問題的拓撲結(jié)構(gòu)優(yōu)化解的存在，在處理技術(shù)上，我們得對該 問題的求解過程加些限制。此處我們采用過濾技術(shù)來進行處理。我們采用網(wǎng)格獨立性濾波器來對單元敏感度的計算來進行改進，其改進 算法如下：de1Nf=i工丑宀卻/=1丁其中Hf為卷積因子(即權(quán)重因子)，其由以下表達式來確定：/C JV|dist(e, /) AnalysisObject Name Shape OptimizationStateFully DefinedDefinitionPhysics TypeStructur

15、al An alysis TypeShape Optimizati onOptionsReference Temp22. Q C在上截面施加 15000N的力Model Shape Optimization LoadsObject NameFixed SupportForceStateFully DefinedScopeScoping MethodGeometry SelectionGeometry1 FaceDefinitionTypeFixed SupportForceSuppressedNoDefine ByVectorMagnitude15000 NDirectionDefined32

16、1.4 求解：32141設(shè)置材料屬性:StructuralYoun gs Modulus 氏模量2.e+011 PaPoiss ons Ratio0.3Den sity 密度7850. kg/m 3Thermal Expa nsion1.2e- 005 1/ CTen sile Yield Stre ngth2.5e+008 PaCompressive Yield Stre ngth2.5e+008 PaTen sile Ultimate Stre ngth4.6e+008 PaCompressive Ultimate Stren gth 0. PaThermalThermal Con duc

17、tivity60.5 W/m CSpecific Heat434. J/kg CElectromag neticsRelative Permeability10000Resistivity1.7e- 007 Ohm-m32142 定義期望目標值的減少量，即期望減少的質(zhì)量或體積分數(shù):ScopeGeometryAll BodiesDefin iti onTargetReduct ion50. %32143 拓撲結(jié)果如下:32144 其他兩種結(jié)果:定義質(zhì)量減少率為 80%ScopeGeometryAll BodiesDefin iti onTargetReducti on80. %拓撲結(jié)果:60%勺

18、拓撲結(jié)果如下定義質(zhì)量減少率為 60%ScopeGeometryAll BodiesDefin iti onTarget Reducti on60. %321.4 試驗結(jié)果分析及結(jié)論：從拓撲結(jié)構(gòu)我們可以看到，在不同的材料去除率下，總的趨勢都 是邊界面梁，其中包含著等強度設(shè)計的思想，這種拓撲結(jié)果較為符合 現(xiàn)實情況，在工程上具有很強的實用性。在不同的材料去除率下，我們發(fā)現(xiàn)結(jié)果具有不同的拓撲結(jié)構(gòu)。3.2.2拓撲結(jié)構(gòu)優(yōu)化：注：以下結(jié)果是在材料去除率為40%勺情況下進行的優(yōu)化。拓撲結(jié)構(gòu)優(yōu)化前，該結(jié)構(gòu)的各項參數(shù)如下：Model Shape Optimization Solution ResultsObjec

19、t NameShape FinderStateSolvedScopeGeometryAll BodiesDefinitionTarget Reduction60. %ResultsOriginal Mass2512. kgOptimized Mass1188.9 kgMarginal Mass35.046 kg3.2.2.1 優(yōu)化分析：在上面Ansys圖中，我們可以發(fā)現(xiàn)，改圖樸結(jié)構(gòu)十分不規(guī)則 工程實際中，這樣的結(jié)構(gòu)在加工上很難實現(xiàn)?；诖?，我們決定對 上面的結(jié)果進行優(yōu)化，得到較為規(guī)則的拓撲結(jié)構(gòu)，以便工程上的加 工。3.2.2.1有限元分析：優(yōu)化前材料去除率為34%BoxLen gth X1.6

20、 mLen gth Y1. mLen gth Z0.2 mPropertiesVolume0.21212 m 3Mass1665.1 kgOrigi nal Mass2512 kgCAD Parametersds_d520ds_s150ds_D1320ds_ss250ds_dd990總位移變形圖應(yīng)變圖應(yīng)力圖優(yōu)化后的分析結(jié)果,Object NameTotal Deformati onEquivale nt Elastic StrainEquivale nt StressMini mum0. m5.4554e-008 m/m10911 PaMaximum5.1147e-006 m3.119e-00

21、6 m/m6.2381e+005 Pa3.221 優(yōu)化設(shè)計：根據(jù)拓撲結(jié)果，我們設(shè)置了 5個參數(shù)作為優(yōu)化設(shè)計變量，其各參數(shù) 取值如下：Design ParametersTABLE 1Coritinuoqs input 卩占rariEt廿對efinitiongIDNameLowei BoundUpper HoundUsable ValuesInitial Valuepeds_d260.780.ContinuciLis520.reds_s75.225.Continuous160.P10ds_D1160430.Continuous320.phdsss125376.Continuous250.P12ds

22、 dd4951485.Continuaus990.我們采用多目標優(yōu)化，將質(zhì)量和結(jié)構(gòu)最大變形作為目標函數(shù)，使其取值最小。在 將多目標轉(zhuǎn)化為單目標時，我們對質(zhì)量和最大變形分別賦予了不同的權(quán)重，具體設(shè)置如 下表：Pararnetei GalsClick iuw/=： ri this labl& to 菇sign(1 四!gg曲Ie lu f碼pun* uarjrr 軌sr. Defining Ldigtl 旳luw is 口pt or 11.NnmifiTtirgieirbesired VmtneIm|inrt4iiiceIraJeOffSDlid Ma強-Mlnirnum P壇ibbIxrwerO

23、nToial Deformation Maxrnuim-Minimum PossibleHigherOn優(yōu)化結(jié)果：PdinoierI* Cilflilid.ilt! A Cdfiilhhlle FC Cftndiikim Cds_d501.0M W564.2de 171 B4143.72*w129.b6dslDl兀 .07235923H.72248 25比龍 A 213.26MJRi290 25ds_dcSuliJ MassTotal Deformation Maximum833.331665.5 kg4.7237e-r06 m1186.91B593 ky 一5.326e-C0S m *8C6

24、.91E63.5 kg5.8108e-006 m通過對比以上三種結(jié)果，我們發(fā)現(xiàn)結(jié)果A較為理想。Rating definitions:The most positive rating of means that the design is excellent in terms of satisfy ing the goal defi ned for the parameter.The most n egative rati ng of means that the desig n is poor in terms of satisfy ing the goal.A rati ng of is n

25、 eutral with respect to the goal.響應(yīng)面：3.222 優(yōu)化實驗結(jié)果:3.221 重要參數(shù)設(shè)計：Nelx=32 x 方向單元的數(shù)目為 32nely=20 y 方向單元的數(shù)目為 20注：該懸臂梁的長寬比為16: 10f/volfrac=0.4 pen al=3.0保留原材料的體積分數(shù)為0.4抑制權(quán)值為3.0 （該取值是資料建議的典型值）=1.2 過濾大小為1.2 （該取值是資料建議的典型值）3.2.2.3試驗結(jié)果分析從該實驗結(jié)果來看，在我們給定的體積保留率的情況下， 每經(jīng)過一次拓撲結(jié)構(gòu)優(yōu)化，該優(yōu)化程序就將懸臂梁的拓撲結(jié)構(gòu) 中強度要求不高處材料的密度減小，直到所有無用

26、的材料都將 被去除為止。我們的拓撲結(jié)構(gòu)優(yōu)化模型是建立在結(jié)構(gòu)變形能最 小、體積去除率自己給定的基礎(chǔ)上進行的，故我們可以根據(jù)實 際情況，自行確定體積去除率。在拓撲優(yōu)化的過程中，我們可以觀察到，我無論體積壓縮 率如何變化，懸臂梁模型最終都向桁架結(jié)構(gòu)進化。這說明，在 結(jié)構(gòu)件中，在自身材料多少相同的條件下，桁架具有很高的剛 度和強度，其實這也就是為什么拓撲結(jié)構(gòu)優(yōu)化首先在桁架結(jié)構(gòu) 領(lǐng)域提出。故工程上，我們常見工程人員采用桁架結(jié)構(gòu)來作為 一些工程的支撐結(jié)構(gòu)，如塔吊等。在實驗中，在驗證不同的體積壓縮率時，我們發(fā)現(xiàn)，在不 同的給定體積壓縮率下，算法的有效性也不同，但在驗證過程 中，發(fā)現(xiàn)算法一直會收斂。體積壓縮率

27、小的時候，該算法能很 快終止；體積壓縮率較大的時候，該算法的收斂速度較慢，并 且還會出現(xiàn)不同程度的震蕩，并且體積壓縮率越大，該算法的 振動也震蕩。3.3實驗結(jié)論及總結(jié)：實驗中，為了充分發(fā)揮我們小組人員的優(yōu)勢，我們從兩條獨立的途徑來 分別對該問題進行研究。An sys 中，我們利用其自帶的的模塊對該問題進行了建模，拓撲結(jié)構(gòu)優(yōu)化。 為了工程上的加工方便，我們在拓撲結(jié)構(gòu)優(yōu)化的基礎(chǔ)上，對其形狀進行了規(guī)則化， 與此同時，進行再次優(yōu)化。Matlab 中，我們將建立的數(shù)學(xué)模型用數(shù)值模擬的形式進行了處理，并將所 得到的處理結(jié)果進行了圖像動態(tài)顯示，以此來清晰的觀察拓撲結(jié)構(gòu)優(yōu)化的動態(tài)過 程，給人以直觀的印象。在用

28、Ansys進行拓撲結(jié)構(gòu)優(yōu)化的時候，我們發(fā)現(xiàn)，當材料去除率為60%寸,其所得到的拓撲結(jié)構(gòu)與我們用 Matlab進行拓撲結(jié)構(gòu)優(yōu)化是所得到的結(jié)果的拓撲 結(jié)構(gòu)是一致的。這驗證了我們的數(shù)學(xué)模型是對的。但是我們的模型也存在著一定的問題，即在材料去除率很小的時候，預(yù)加 載荷未變的情況下，該模型的解均存在。出現(xiàn)這一情況的原因是，我們在計算懸 臂梁的位移矩陣，一直認為其應(yīng)力未超過破壞強度，整個過程均在彈性變形圍。 故在實際應(yīng)用的時候，需要我們結(jié)合一定的工程經(jīng)驗來確定個重要參數(shù)的取值。四、參考文獻及附錄：4.1主要參考文獻：1. Sigmund, O. 1994: Design of material struc

29、tures using topologyoptimization. Ph.D. Thesis, Department of Solid Mechanics,Tech ni cal Uni versity of Denmark2. Sigm und, O. 1997: On the desig n of complia nt mecha ni smsusing topology optimization. Mech. Struct. Mach. 25, 495- 5263. ANSYS/WORKBE設(shè)計、仿真與優(yōu)化兵何正嘉雪峰等清華大學(xué)2008年8月第一版4. .topopt.dtu.dk .5

30、. 基于結(jié)構(gòu)最大剛度的形狀優(yōu)化方法歐陽高飛建民機械工程學(xué)報2008年10月4.2相關(guān)附件：附件1:懸臂梁拓撲結(jié)構(gòu)優(yōu)化Matlab程序：% A 99 LINE TOPOLOGY OPTIMIZATION CODE BY OLESIGMUND, OCTOBER 1999 %fun ctio n top( nelx, nely,volfrac,pe nal,rmi n);%主程序，即調(diào)用程序。其中的參數(shù)需要用戶提供。依次為x方向單元數(shù)、y方向單元數(shù)、期望優(yōu)化的體 積或質(zhì)量百分比、抑制權(quán)值、過濾因子% INITIALIZE 初始化x(1: nely,1: nelx) = volfrac; %在整個求解

31、域均分材料loop = 0;cha nge = 1.; % 新舊改變之差% START ITERATIONwhile cha nge 0.01 %第一層，設(shè)置圖像顯示次數(shù)，即優(yōu)化的次數(shù)loop = loop + 1;xold = x;% FE-ANALYSIS 有限元U=FE( nelx, nely,x,pe nal); %調(diào)用有限元子程序，返回位移矢量% OBJECTIVE FUNCTION AND SENSITIVITY ANALYSISKE = lk; %調(diào)用單元剛度陣c = 0.;for ely = 1:nely %16-24行一次循環(huán)。外層循環(huán)，y方向循環(huán)nely次數(shù)for elx

32、= 1:nelx %層循環(huán)，x方向循環(huán)次數(shù)為nelxn1 = (n ely+1)*(elx-1)+ely; %在全局坐標系中，左上方的節(jié)點數(shù)目n2 = (n ely+1)* elx +ely;%在全局坐標系中，右上方的節(jié)點數(shù)目n1、n2作用是從全局位移矢量提取局部或單元位移矢量Ue =U(2* n1-1;2* n1;2*n 2-1;2* n2;2*n 2+1;2* n2+2;2*n 1+1;2*n1+2,1);% 返回單元位移c = c + x(ely,elx)Ape nal*Ue*KE*Ue;dc(ely,elx) = -pe nal*x(ely,elx)A(pe nal-1)*Ue*KE*

33、Ue;endend% FILTERING OF SENSITIVITIES 靈敏度檢驗dc = check( nelx ,n ely,rmi n,x,dc); %調(diào)用網(wǎng)格濾波子函數(shù)% DESIGN UPDATE BY THE OPTIMALITY CRITERIA METH 優(yōu)化x = OC( nelx, nely,x,volfrac,dc); %調(diào)用優(yōu)化子函數(shù)% PRINT RESULTS拓撲過程顯示cha nge = max(max(abs(x-xold);disp( 循環(huán)次數(shù) It.: sprintf(%4i,loop)優(yōu)化目標值 Obj.:spri ntf(%10.4f ,c).拓撲質(zhì)

34、量分數(shù) Vol.: sprintf(%6.3f,sum(sum(x)/(nelx*nely).終止條件ch.: sprintf(%6.3f,change ) %優(yōu)化過程指標顯示% PLOT DENSITIES 密度分布colormap(gray); imagesc(-x); axis equal; axis tight; axis off;pause(1e-6); %圖像顯示設(shè)置 end% MESH-INDEPENDENCY FILTER % function de n=check (n elx ,n ely,rm in, x,dc)dcn=zeros( nely,n elx);for i =

35、 1: nelxfor j = 1: nelysum=0.0;for k = max(i-ro un d(rmi n),1):mi n(i+ro un d(rmi n), nelx)for l = max(j-ro un d(rmi n),1):mi n(j+ro un d(rmi n), n ely)fac = rmin-sqrt(i-k)A2+(j-l)A2);sum = sum+max(0,fac);dcn (j,i) = dcn(j,i) + max(0,fac)*x(l,k)*dc(l,k);endenddcn (j,i) = dcn (j,i)/(x(j,i)*sum);enden

36、d% FE-ANALYSIS %fun ctio n U=FE( nelx ,n ely,x,pe nal)KE=lk;K = sparse(2*(nelx+1)*(nely+1), 2*(nelx+1)*(nely+1);F=sparse(2*( nely+1)* (n elx+1),1);U=sparse(2*( nely+1)*( nelx+1),1);for ely = 1:n ely %全局剛度矩陣的形成for elx = 1:n elxn1 = (n ely+1)*(elx-1)+ely;n2 = (n ely+1)* elx +ely;edof = 2*n1-1; 2*n1; 2

37、*n2-1; 2*n2; 2*n2+1;2*n2+2;2*n1+1; 2*n 1+2;K(edof,edof) = K(edof,edof) +x(ely,elx)Ape nal*KE;endend% DEFINE LOADSAND SUPPORTS(HALF MBB載EA的施加和邊界條件的定 義%F(2,1) = -1; % 在單元力在左上方%fixeddofs = un io n(1:2:2*( nely+1),2*( nelx+1)*( nely+1);F(2*( nelx)*( nely+1)+2,1) = -1;fixeddofs = 1:2*( nely+1);alldofs =

38、1:2*( nely+1)*( nelx+1);freedofs = setdiff(alldofs,fixeddofs);% SOLVINU(freedofs,:) = K(freedofs,freedofs) F(freedofs,:); %添加約束條件，限制自由度U(fixeddofs,:)= 0; %freedofs沒有約束的自由度，fixeddofs 受到約束的自由度% ELEMENT STIFFNESS MATRIX %fun ctio n KE=lkE = 2.;nu = 0.3;k=1/2-nu/6 1/8+nu/8 -1/4-nu/12 -1/8+3*nu/8 .-1/4+n

39、u/12 -1/8-nu/8 nu/6 1/8-3*nu/8;KE = E/(1- nu A2)*k(1) k(2) k(3) k k(5) k(6) k(7) k(8)k(2) k(1) k(8) k(7) k(6) k(5) k k(3)k(3) k(8) k(1) k(6) k(7) k k(5) k(2)k k(7) k(6) k(1) k(8) k(3) k(2) k(5)k(5) k(6) k(7) k(8) k(1) k(2) k(3) kk(6) k(5) k k(3) k(2) k(1) k(8) k(7)k(7) k k(5) k(2) k(3) k(8) k(1) k(6

40、)k(8) k(3) k(2) k(5) k k(7) k(6) k(1);% OPTIMALITY CRITERIA UPDATE %基于標準優(yōu)化方法的優(yōu)化程序fun ctio n x new=OC( nelx, nely,x,volfrac,dc) %返回最新的變量值，不斷更新變量l1 = 0; l2 = 100000; move = 0.2;while (12-11 1e-4)lmid = 0.5*(12+11);xnew=max(0.001,max(x-move, min (1.,m in (x+move,x.*sqrt(-dc./lmid);if sum(sum(x new) - v

41、olfrac* nelx* nely 0; %bi-sect ioning algorithm11 = lmid;else12 = lmid;endend% A 99 LINE TOPOLOGY OPTIMIZATION CODE BY OLESIGMUND, OCTOBER 1999 %function top(nelx,nely,volfrac,penal,rmin);% INITIALIZEX(1 me_yAne_x) H vo_faG-OOP H o八change 丄八% START _TERAT_ONwhi-e change V 0.2-OOP H -OOP +八xo-d H7cHF

42、E(ne_pena-*ue-*KE*updc(e-ykx) H pena-*x(e-ykx)(penaL)*uQ.*KE*up endend% F-LTER-NG OF SENS_T_v_Tms-del H checkmeFne-yJminxdc)八% DESON UPDATE BY THE OPT_M A 匚 TY CR_TER_A METK6OS s H oc(ne_ne-yxNO-facQ.c=% PR_NT RESULT也芋K就越引change H max(max(abs*xo-d)八disp(r _r- sprin戈 4-_00p) - obj:.5|35%10.4o) :-stw幷

43、s午w幷_mid ho.5*(_2+= =xnewHmax(oo2-max*movemn(rmin(x+movpx.*sqrt(dcj-mid)八if sum(sum(x new) - volfrac* nelx* nely 0;11 = Imid;else12 = lmid;endend% MESH-INDEPENDENCY FILTER % function de n=check (n elx ,n ely,rm in, x,dc)dcn=zeros( nely,n elx);for i = 1: nelxfor j = 1: nelysum=0.0;for k = max(i-ro un d(rmi n),1):mi n(i+ro un d(rmi n), nelx)for l = max(j-ro