基于MATLAB和Optistruct的C形夾拓撲優(yōu)化

1、華中科技大學研究生課程考試答題本謹于MATLAB和Optistruct的C形夾拓撲優(yōu)化學院（部）：機械科學與工程學院課程名稱： 工程優(yōu)化設(shè)計學生姓名：范利洪班 級：機碩11 07班學M201170602號：2012年01月1 0日目錄第1章選題背景介紹及問題描述。錯誤!未定義書簽。1 .1選題背景及意義。錯誤!未定義書簽。1.1.1 工程背景及基本原理。錯誤!未定義書簽。1.1.2 本文研究意義。錯誤!未定義書簽。1.2 研究現(xiàn)狀錯誤!未定義書簽。1.2. 1理論研究現(xiàn)狀錯誤!未定義書簽。1.3. 2應(yīng)用研究現(xiàn)狀。錯誤!未定義書簽。1.3 該研的意義4第2章SIMP變密度法理論基礎(chǔ)錯誤!未定義

2、書簽。2. 1SIM P密度剛度插值法理論基礎(chǔ)錯誤!未定義書簽。2.1 拓撲優(yōu)化的數(shù)學模型。錯誤!未定義書簽。2.2 優(yōu)化準則的基本理論。錯誤!未定義書簽。 第3章優(yōu)化設(shè)計的數(shù)學模型及解決方案。錯誤!未定義書簽。3. 1問題描述錯誤!未定義書簽。3.1 優(yōu)化問題的數(shù)學模型錯誤!未定義書簽。3.2 模型分析求解方法選擇錯誤!未定義書簽。第四章拓撲優(yōu)化步驟及結(jié)果錯誤!未定義書簽。4. 1基于matlab的變密度拓撲優(yōu)化。錯誤!未定義書簽。4.1 1問題求解的關(guān)鍵技術(shù)及代碼。錯誤!未定義書簽。4. 1.2拓撲優(yōu)化結(jié)果及分析錯誤!未定義書簽。4.2基于Optistruct的C形夾拓撲優(yōu)化。錯誤!未定義

3、書簽。1 .2.1有限元模型的建立。錯誤!未定義書簽。4 . 2 .2基于Opti s truct的拓撲優(yōu)化結(jié)果及分析。錯誤!未定義書簽。第6章 結(jié)論及總結(jié)錯誤!未定義書簽。參考文獻錯誤!未定義書簽。附件：懸臂梁拓撲優(yōu)化mat 1ab程序錯誤!未定義書簽。第1章 選題背景介紹及問題描述1.1 選題背景及意義1.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ū)域

4、內(nèi)尋求最優(yōu)材料的分 布問題。尋求一個最佳的拓撲結(jié)構(gòu)形式有兩種基本的原理:一種是退化原理，另一 種是進化原理。退化原理的基本思想是在優(yōu)化前將結(jié)構(gòu)所有可能桿單元或所有材 料都加上，然后構(gòu)造適當?shù)膬?yōu)化模型,通過一定的優(yōu)化方法逐步刪減那些不必要 的結(jié)構(gòu)元素，直至最終得到一個最優(yōu)化的拓撲結(jié)構(gòu)形式。進化原理的基本思想是 把適者生存的生物進化論思想引入結(jié)構(gòu)拓撲優(yōu)化，它通過模擬適者生存、物競天 擇、優(yōu)勝劣汰等自然機理來獲得最優(yōu)的拓撲結(jié)構(gòu)。1.1.2 本文研究意義目前，結(jié)構(gòu)優(yōu)化大部分集中在尺寸設(shè)計變量（如板厚、桿的剖面積及管梁的 直徑）。拓撲結(jié)構(gòu)優(yōu)化較尺寸優(yōu)化復(fù)雜，但對于有些問題拓撲結(jié)構(gòu)優(yōu)化比尺寸 優(yōu)化有效，c

5、形夾是其中的例子之一。本文討論c形夾的拓撲優(yōu)化問題，圍繞這 一問題，怎樣使結(jié)構(gòu)具有最大剛度的設(shè)計占有相當重要的地位;怎樣優(yōu)化結(jié)構(gòu)的形 狀使材料的分布,更加合理從而達到使結(jié)構(gòu)具有最小柔度的目的是本文要研究的 問題。1-2研究現(xiàn)狀1.2 . 1理論研究現(xiàn)狀結(jié)構(gòu)拓撲優(yōu)化是近20年來從結(jié)構(gòu)優(yōu)化研究中派生出來的新分支，它在計算 結(jié)構(gòu)力學中已經(jīng)被認為是最富挑戰(zhàn)性的一類研究工作。目前有關(guān)結(jié)構(gòu)拓撲優(yōu)化的 工程應(yīng)用研究還很不成熟，在國外處在發(fā)展的初期，尤其在國內(nèi)尚屬于起步階段。 1 904年M i chell在桁架理論中首次提出了拓撲優(yōu)化的概念。自1964年Dorn 等人提出基結(jié)構(gòu)法，將數(shù)值方法引入拓撲優(yōu)化領(lǐng)域

6、，拓撲優(yōu)化研究開始活躍。20世 紀8 0年代初，程耿東和N.O 1 ho f f在彈性板的最優(yōu)厚度分布研究中首次將最 優(yōu)拓撲問題轉(zhuǎn)化為尺寸優(yōu)化問題，他們開創(chuàng)性的工作引起了眾多學者的研究興 趣。1 9 88年Be n d so e和Ki ku c h i發(fā)表的基于均勻化理論的結(jié)構(gòu)拓撲優(yōu)化設(shè) 計，開創(chuàng)了連續(xù)體結(jié)構(gòu)拓撲優(yōu)化設(shè)計研究的新局面。I 9 93年Xi e YM和S t evenGP提出了漸進結(jié)構(gòu)優(yōu)化法。1999年Bendsoe和S i gmund證實了變密度法 物理意義的存在性。2002年羅鷹等提出三角網(wǎng)格進化法，該方法在優(yōu)化過程中實 現(xiàn)了退化和進化的統(tǒng)一，提高了優(yōu)化效率。1.3 . 2應(yīng)用

7、研究現(xiàn)狀在前人提出的重要理論基礎(chǔ)上，后人也將其跟其他現(xiàn)代設(shè)計的方法相結(jié)合， 衍生出了其他一些拓撲結(jié)構(gòu)優(yōu)化方法：如與可靠性相結(jié)合的情況下，MAUTE 等應(yīng)用變密度法并結(jié)合可靠性分析對一微機電系統(tǒng)進行了基于可靠性的拓撲優(yōu) 化設(shè)計,PAPADRA KAKIS等將遺傳算法應(yīng)用于具有可靠性約束的桁架結(jié)構(gòu)拓 撲優(yōu)化設(shè)計中，國內(nèi)學者馬洪波也對基于遺傳算法的結(jié)構(gòu)可靠性優(yōu)化問題進行了 討論。華南理工大學機械工程學院歐陽高飛等對基于水平集方法的結(jié)構(gòu)可靠性拓 撲優(yōu)化進行了研究。1.4 本文研究的意義通過這次的作業(yè)加深對工程優(yōu)化算法的學習和使用，提高對拓撲優(yōu)化的方法 和過程的了解和學習。另外對相關(guān)軟件軟件的應(yīng)用能夠達

8、到一個新的高度。這些 不僅能使我們現(xiàn)在的知識體系得到充實和優(yōu)化，而且也是我們今后人生的財富。第2章 SIMP變密度法理論基礎(chǔ)2.1 SIMP密度剛度插值法理論基礎(chǔ)S I M P模型主要通過引入懲罰因子，在材料的彈性模量和單元相對密度之 間建立起一種顯示的非線性對應(yīng)關(guān)系。它的作用是當設(shè)計變量的值在(0, 1 )之 間時.，對中間密度值進行懲罰,使中間密度值逐漸向0/1兩端聚集，這樣可以使 連續(xù)變量的拓撲優(yōu)化模型能很好地逼近原來01離散變量的優(yōu)化模型。這時中 間密度單元對應(yīng)一個很小的彈性模量，對結(jié)構(gòu)剛度矩陣的影響將變得很小，可以 忽略不計。SIMP材料模型的數(shù)學表達形式：¥但)二爐加+必

9、(£° E巧JJ1) .1)的=之(即+叼3儲(2.2 )其中:p為兩數(shù)學模型中對中間密度材料的懲罰因子。懲罰因子的作用是當 設(shè)計變量的值在(0, 1)之間時，通過逐漸增加的值對設(shè)計變量的中間值進行懲 罰，隨著值的增大，設(shè)計逐漸接近0/1設(shè)計。為有效壓縮中間密度材料，要 求 2 2。EP表示插值以后的彈性模量,石。和Em,n分別為固體和空洞部分材料的 彈性模量，AE = E°-Em,n,y'-5。/ ooo。,(/ = 1,2,表示單元/的設(shè)計 變量。K表示插值以后的剛度矩陣,K/表示笫，個單元固體材料的剛度矩陣。2) 2拓撲優(yōu)化的數(shù)學模型以結(jié)構(gòu)的柔度最小

10、化(或剛度最大化、應(yīng)變能最小化)作為優(yōu)化的目標函數(shù)， 以結(jié)構(gòu)整體的體積約束作為優(yōu)化的約束條件。剛度優(yōu)化的數(shù)學模型表示為：'Minimize :C(X) = UrKUX = $ 也,./ P 齡廣小Subject to J 0 < x而門-xj - 1，KU=F.(2.3)S I MP對應(yīng)的柔度函數(shù)和敏度形式：c(x) = X(£min+ 考 AE)U%j=i(2.4)Cx) = £ px/TUAEj=i(2. 5)其中:以上各式中，KJ表示第j個單元的剛度矩陣。U表示結(jié)構(gòu)的位移向 量；X表示設(shè)計變量，為避免總剛度矩陣奇異，取Xmin=0-00 1。為單元數(shù)目，

11、 。表示結(jié)構(gòu)的柔度,C'表示柔度關(guān)于設(shè)計變量的敏度。2.3優(yōu)化準則的基本理論剛度拓撲優(yōu)化問題，是典型的具有不等式約束的非線性規(guī)劃問題。不等式約 束多元函數(shù)極值的必要條件是Kuhn-Tucker條件，它是采用優(yōu)化準則法求解 非線性優(yōu)化問題的重要理論。引入對設(shè)計變量上下限約束的拉格朗日乘子力、萬 以及對體積約束的拉格朗日乘子A ,構(gòu)造Lagrange函數(shù)L為如下形式： 乙=。(%)+ 八43廠 0) +，(KU-用)+ £用(-七)+ 支片區(qū)72) j=ij=j=i(2.6)對于La grange函數(shù)，當x)時,設(shè)計變量的上下限側(cè)面約束均不起 作用，設(shè)計變量是主動變量。主動變量在

12、迭代過程中作為設(shè)計變量允許發(fā)生改 變,右=片=0;當/加之與時，僅設(shè)計變量下限約束起作用，設(shè)計變量為被動變 量,*20,片=0；當心a、",時，僅設(shè)計變量上限約束起作用，設(shè)計變量也為被動 變量，2>0,2；>0o被動變量在迭代過程中不能變化，只能由側(cè)面約束的邊界 值來確定。第3章優(yōu)化設(shè)計的數(shù)學模型及解決方案31問題描述如圖所示，C形夾在自由端口受到三角形分布力F的作用，要求在保持對原 材料體積一定縮減比的情況下對原實體懸臂梁做結(jié)構(gòu)拓撲優(yōu)化設(shè)計，優(yōu)化目標是 使結(jié)構(gòu)剛度最大。（優(yōu)化的結(jié)果應(yīng)該使原設(shè)計區(qū)域產(chǎn)生孔洞，使結(jié)構(gòu)拓撲發(fā)生變 化。）原實體C形夾為如下圖3. 1所示的C形，

13、尺寸如下圖所示，材料為4 5鋼, 密度p為7.8x/加3，彈性模量E = 2e5MPa，泊松比 =0.3。F = 20N。3.2優(yōu)化問題的數(shù)學模型該問題中，要求同時滿足剛度最大，質(zhì)量最輕，這兩個變量若同時改變，則問 題復(fù)雜度太大，并且可能導(dǎo)致問題不可求解。所以我們采用在確定的質(zhì)量下，來討 論剛度最大的問題。由于對特定的材料，其質(zhì)量和體積有一定的關(guān)系，并且我們 采用去除法的思想來建立模型的，故我們可以采用給優(yōu)化后的體積與優(yōu)化前的體 積比賦確定的值，來達到在給定質(zhì)量條件下滿足剛度最大的問題。其數(shù)學模型如 下：NMinimize :C(X) = UKU二X =( X，"2'，芍 7

14、0=1vA-v<o,>1Subject to: <°V/in« 外 41， 7 = l,2,.,n KU=F.(1)注：其中C(X)為結(jié)構(gòu)變形能,U為結(jié)構(gòu)變形總位移矩陣，K為結(jié)構(gòu)總剛度矩 陣,N為劃分單元總數(shù),% 為單元位移向量,。為單元剛度，(由于劃分單元的時 候，我們采用等分矩形單元，所以每個單元的剛度可用一個常量來處理)匕是拓 撲結(jié)構(gòu)優(yōu)化過程中變化著的體積,0為未經(jīng)過優(yōu)化前懸臂梁的體積。F為結(jié)構(gòu)所 受的三角形載荷。，為懸臂梁的相對密度。3.3模型分析求解方法選擇對該問題是用兩種方法求解。方法一:基于mat lab的變密度拓撲優(yōu)化法; 方法二:是用成熟

15、的有限元拓撲優(yōu)化軟件Optis t ruct進行優(yōu)化?；趍atlab的變密度拓撲優(yōu)化該問題的優(yōu)化方法有很多種，常用的有如下方法：Optima 1 it y Crite r ia (OC) methods,(優(yōu)化準則方法)Sequentia 1 Linear P rogramm i ng (SLP ) meth o ds (序列線性規(guī)劃法) Meth o d of Moving Asymp t o tes (MM A by S van ber g 1987)等 為了簡化問題的復(fù)雜度,此處我們采用stand a rd 0 C-meth。d .方法來 實現(xiàn)。在處理過程中，關(guān)于設(shè)計變量相對密度x每一

16、步的更新，我們采用灰公四 在1 99 5年提出的如下算法來實現(xiàn)。第四章拓撲優(yōu)化步驟及結(jié)果4.1基于ma t la b的變密度拓撲優(yōu)化4.1. 1問題求解的關(guān)鍵技術(shù)及代碼一般而言，由于OC法所使用的單元是矩形，所以0C法很適合求解求解域 為矩形的優(yōu)化問題，而本文選題為一 C形結(jié)構(gòu)，若使用OC法，則需將C形結(jié) 構(gòu)劃分成為三個矩形的集合，因而在整體剛度矩陣的組裝方面應(yīng)該考慮如何進 行，由于本程序是參考經(jīng)典9 9行OC法拓撲優(yōu)化程序改變，對于程序中避免邊界 鋸齒現(xiàn)象所使用的check函數(shù)，該如何進行修改，以及如何進行邊界條件的施加。優(yōu)化問題的初值條件：Nelx=6 0 x方向單元的數(shù)目為60ne 1

17、y=50y方向單元的數(shù)目為50vo 1 f rac = 0.35保留原材料的體積分數(shù)為0 . 35p e nal=3. 0抑制權(quán)值為3 .0 (該取值是資料建議的典型值)rmin=l. 2過濾大小為1. 2 (該取值是資料建議的典型值)關(guān)鍵代碼：funct i o n xnew=OC(n e lx, nely,x, volfrac,de)%定義 OC 優(yōu)化準則函數(shù)11 = 0 ； I 2 = 1000 0 0； mo v e = 0 . 2；while (12-11 > 1 e - 4 )Imid = 0.5 *( 1 2+11);%采用折半查找x n e w = m a x (0. 0

18、 01,max(x m o v e , min(l.jnin( x +mo v e, x.*s q r t (-de. / Imid)；%根據(jù)Xj和敏度de更新X的值%if s u m(sum(xnew) - volfrac*n e 1 x * n e ly > 0 ;%根據(jù)迭代過程 中體積比是否達到預(yù)設(shè)的體積比，判斷迭代是否繼續(xù)進行%11 = Imid;else1 2 = Im i d;endendfor e 1 y = 1 :nelyfor elx = I :( n e 1 x-20)%取C形夾的左半矩形區(qū)域n 1 = (ne 1 y+1) * (el x -1) +e 1 y;n2

19、 = (nely+1)* elx +ely;Ue = U(2*n 1 -l;2*n 1 ; 2*n2-l;2 * n 2; 2 * n2+l ;2 * n2+2; 2*nl+ 1 ;2*n 1+2J);%單元位移向量 %x( 1 6： 3 5,41： n elx ) =0.001;c = c + x( e 1 y ,e 1 x )- penal*U e，* KE* Ue;% 目標函數(shù) dc(el y , e 1 x ) = - p enal* x (eljse 1 x)A( penal-1 )*Ue* KE*Ue;% 敏度值deendend4.1.2拓撲優(yōu)化結(jié)果及分析實驗結(jié)果分析：從該實驗結(jié)果

20、來看，在我們給定的體積保留率的情況下,每經(jīng)過一次拓撲結(jié)構(gòu) 優(yōu)化，該優(yōu)化程序就將C形夾的拓撲結(jié)構(gòu)中強度要求不高處材料的密度減小，直 到所有無用的材料都將被去除為止。我們的拓撲結(jié)構(gòu)優(yōu)化模型是建立在結(jié)構(gòu)變形 能最小、體積去除率自己給定的基礎(chǔ)上進行的，故我們可以根據(jù)實際情況，自行確 定體積去除率。在拓撲優(yōu)化的過程中，我們可以觀察到，我無論體積壓縮率如何變化,c形 夾模型最終都向桁架結(jié)構(gòu)進化。這說明，在結(jié)構(gòu)件中，在自身材料多少相同的條 件下,桁架具有很高的剛度和強度，其實這也就是為什么拓撲結(jié)構(gòu)優(yōu)化首先在桁架 結(jié)構(gòu)領(lǐng)域提出。故工程上，我們常見工程人員采用桁架結(jié)構(gòu)來作為一些工程的支 撐結(jié)構(gòu)，如塔吊等。在驗證

21、不同的體積壓縮率時，在不同的給定體積壓縮率下，算法的有效性也 不同，但在驗證過程中，發(fā)現(xiàn)算法一直會收斂。體積壓縮率小的時候，該算法能 很快終止;體積壓縮率較大的時候，該算法的收斂速度較慢,并且還會出現(xiàn)不同程 度的震蕩，并且體積壓縮率越大，該算法的振動也震蕩。4.2基于Optis truct的C形夾拓撲優(yōu)化4.2.1 有限元模型的建立在Optis true t中建立有限元模型并劃分網(wǎng)格。相關(guān)參數(shù)為：E=2xl05：4=03；密度片7800口/ m3添加邊界調(diào)節(jié)和載荷，得到有限元模型如下圖5. 1所示：圖5.1有限元模型4. 2.2基于O p t i s t met的拓撲優(yōu)化結(jié)果及分析當體積分數(shù)v

22、ol frac =0. 35時，拓撲結(jié)果圖如圖5. 2所示:Rssult E /tmp/CMAX/cxingjial-des.KJdDesign Iteration 30Fram31圖5.2 C形夾拓撲優(yōu)化結(jié)果圖(vol f r a c=0.3 5 ) 當體積分數(shù)v olf r ac= 0 . 5時，拓撲結(jié)果圖如圖5 .3所示：Result E./tmp/CMAX/ cxinqj ia1 _dc5. h3dDesign ： Iteration 30LFrame 31圖5.3 C形夾拓撲優(yōu)化結(jié)果圖(volfrac= 0 .5)拓撲優(yōu)化結(jié)構(gòu)分析：在以變形最小為目標函數(shù)，即結(jié)構(gòu)剛度最大，一定的材 料

23、保留比率volf r a c為約束調(diào)節(jié)，當v o Ifrac取不同的值時，優(yōu)化得到的拓撲結(jié) 構(gòu)基本保持不變，說明該拓撲結(jié)構(gòu)存在一理想狀態(tài)，通過不斷嘗試，得到如下發(fā)現(xiàn): 在本例中當vol frac=O. 35時所得的拓撲結(jié)構(gòu)最為規(guī)律，并且該拓撲結(jié)構(gòu)趨向 與桁架結(jié)構(gòu)。(具體拓撲結(jié)構(gòu)演變過程見附件視頻動畫)第6章結(jié)論及總結(jié)本文中，從兩條獨立的途徑來分別對該問題進行研究。Hyper works中，利用hype r mesh自帶的的模塊對該問題進行了建模，網(wǎng) 格劃分，得到有限元模型，然后將有限元模型導(dǎo)入拓撲優(yōu)化模塊Opt i st r u c t, 添加邊界條件及約束，定義目標函數(shù)和約束條件，對C形夾進

24、行拓撲結(jié)構(gòu)優(yōu)化。 所得到的拓撲結(jié)構(gòu)趨向于桁架，說明桁架結(jié)構(gòu)有很高的剛度。Matlab中，首先對該問題構(gòu)建數(shù)學模型，采用SI MP(OC)算法對數(shù)學模 型進行優(yōu)化求解，由于本文所研究的問題的求解域并不是規(guī)則的矩形，而是由三 部分矩形疊加，所以在使用SIMP法時，對總體剛度矩陣組裝，添加邊界和約束 條件成為了本文的關(guān)鍵技術(shù)，涉及到較多的有限元理論。在mat 1 a b中將每次迭 代所得到的處理結(jié)果進行了圖像動態(tài)顯示，以此來清晰的觀察拓撲結(jié)構(gòu)優(yōu)化的動 態(tài)過程，給人以直觀的印象。在用Optistru c t進行拓撲結(jié)構(gòu)優(yōu)化的時候，我們發(fā)現(xiàn)，當材料去除率為 60%時，其所得到的拓撲結(jié)構(gòu)與我們用Matla

25、b進行拓撲結(jié)構(gòu)優(yōu)化是所得到的結(jié) 果的拓撲結(jié)構(gòu)是一致的。這驗證了我們的matlab數(shù)學模型是對的。參考文獻lO.sigm u nd. A 9 9 lin e t opo 1 og y optim i za t ion cod e writ t en in Ma t 1 a b.Ed u cational ar t i c 1 e , 2 0 0 12羅震.基于變密度法的連續(xù)體結(jié)構(gòu)拓撲優(yōu)化設(shè)計技術(shù)研究.博士學位論文.20 0 53 孫靖民.機械優(yōu)化設(shè)計第三版.機械工業(yè)出版社.200 4白新理等.結(jié)構(gòu)優(yōu)化設(shè)計.黃河水利出版社.2 0 0 8.4郭仁生.機械工程設(shè)計分析和MATLAB應(yīng)用.機械工業(yè)出版

26、社.200 6附件：C形夾拓撲優(yōu)化matl a b程序fu nction top(nc 1 x ,ne 1 y.volf r ac, p e n al. r mi n);%變量初始化n elx= 6 0;nel y = 5 0:vo 1 frac=0.3 5;penal = 3 ;rmi n = 2 ;x(l:nely.1 : nelx)= v olfra c ;1 oop =0;change = L;%并始迭代wh i le cha nge> 0.01loop = loop + 1;xold = x;U=FE( n elx. nely # x.pe nal) ;%調(diào)用子函數(shù)得到整體位移

27、向量%將求解域分成三塊矩形進行目標函數(shù)定義和敏度分析KE = lk;c = 0.;for ely = 1 :ne 1 yfor elx = l:(nelx-20)nl = ( n e ly + 1 ) *(elx - 1 ) +ely ;n2 = (ne 1 y + 1 )* e 1 x +ely;Ue =U(2*nl-l;2*nl;2*n2-l ； 2 *n 2 ; 2*n2+l:2*n2 + 2;2*n 1 + 1 ; 2*nl+2J )；x(16:35 , 4 1:nelx)=0 . 001;C = C +乂化1丫，0區(qū))、口1*110'*1<£*110；%目標函

28、數(shù)dc(ely,elx )= - pen a l*x(elyxlx)A(pena 1 -1)*U e1 * KE * Ue;%敏度endendfor e 1 y = 1:15f orelx = 41 : nelxnl = (nely+1) * ( e 1 x-1 ) +e 1 y ;n2 = (nel y+1)* elx +e 1 y;Ue = U(2*nl4;2 *n 1; 2*n2-l;2*n2; 2*n2+l ； 2*n2 + 2;2*nl + l;2*nl+2 ,1)；x(16:35,41:nelx ) = 0.0 0 1;c = c + x(elyxIx)A p enal * Ue

29、1 *KE* U e;% 目標函數(shù)de ( ely, elx ) = -penal * x(ely z e lx)A(penal-l)*U eK E * Ue ;endendf o r el y = 36: ne 1 yfo r el x = 4 l:nelxnl = ( nely+l)*(e 1 x-l)+ely;n2 = (ncly+l)*elx + ely;Ue = U(2*nl-l;2*nl; 2*n2-l;2*n2; 2*n2 + l;2*n2 + 2; 2*nl + l;2 *nl+2,l );x(16:3 5 , 41 ： nelx)=O . 001;c = c + x ( el

30、y9e 1 x ) A p e n a l*Ue 1 *KE*Ue ; % 目標函數(shù) dc(elyxlx) = -penal * x (e 1 y,el x ) A(pe n al-l)*U e "*KE*U e;end end%敏度過濾，避免拓撲結(jié)構(gòu)邊界鋸齒傳de = chec k(nelx, n elyjm i n r x z d c );%通過OC法更新單元密度x= OC(nelx,n e 1 y,x,v ol f rac.dc);fo r j = (nelx-19) : ( ne 1 x-1 ) %指定非優(yōu)化區(qū)域x(15j)=0.9;endf orj =(ne 1 x-1 9

31、):(ne 1 x-1)x(3 6 J) = 0.9 ;endc h a n g e = max(max(abs( x - xold);disp(It.: 1 sp r intf('%4i 1 oop)1 Obj.:1 spr i n t f (1%10.4fx)1 V oL: * sp r i ntf( 1 % 6 . 3f, s u m(s u m(x ) )/ ( n elx* n ely)* ch.: 1 s p ri n tf(r% 6.3f,c h a n ge ) ) %顯示苦果colorma p(gray); imagesc( - x ) ; axis eq u a 1

32、 ; a x i s tig ht; axis of f ;p a use(l e-6);%繪制密度鹵end%最優(yōu)性判據(jù)更新單元密度值%f unc t ion xnew=OC(nelx,nely,x , volfrac,dc)11 = 0: 12= 10000 0; m o ve = 0.2;whi 1 e (12-11 > le - 4)Imid = 0.5*( 12+1 1 );xnew = max ( 0.0 0 1 .max(x-mo ve,min( 1 .,min(x + mo ve,x .*sq r t( - dc./I mid);1 f s uni ( suni(xnew)

33、 - v olfrac*nelx * nel y > 0;11 = Imid;else12 = Imi d ;endend% 單元密度過濾 %f unctio n de n =check(nelx. n elyjmi n , x , de)d c n = z e r o s(n e ly f ne 1 x );for i = 1 :(nel x - 20)f o r j = l:nelysum = 0.0 ;fo r k = max(i - floor(rmin)J ) :min(i+ floor (rm in), n elx)for 1 = max (j - floor(rmin),

34、1) ： min(j + floor(rmin)jiely)f ac = rmi n-sqrt(i- k)A2 +(j-I) A2);sum = sum+inax( 0 ,fac)；den ( j , i )= den (j J) + m a x(0 r fac) * x(l,k)*dc( 1, k )；endendd c n (j,i) = d c n(j,i)/(x( j j)*s u m);endendfor i = (nelx-19) : nelxfo r j = 1:15sum= 0.0 ;for k= max(i-floo r (rm i n),l) : m i n (i+floo

35、r (rmin),nelx)f orl = max (j - floor ( r min),l):min(j + floor(rmin), nely)fac = rmin-s q r t(i - k)A2+(j-l) A 2)；s u m = s u m+m a x(0 , f a c);dcn(j, i ) = dcn( j J ) + max(0 r fac ) *x(l.k) * dc( 1 f k );endendde n (jj) = dcn(j , i) /(x (j, i ) * sum);endendf ori = (nelx - 19) : ne 1 xfor j = (ne

36、 1 y-14):ne 1 ys u m = 0.0;for k = max(i-flo or( r min)J):m i n(i+floor(rmin ) ,nelx)for 1 = ma x(j-floor( r min), 1): min (j+floo r ( r min ) ,nely )fac = r m i n-sq r t ( i -k ) A 2+ (j-I)A2);sum = s um+max( 0 Jac)；dcn(j, i ) = de n (j ,i) + m a x ( OJac) * x(l , k)*dc( 1 / k);ende nddcn(j, i ) =

37、dcn( j ,i)/(x (j , i )*sum );endend% 有限元分析 %f un c tion U=FE ( nel x # nely. x.pen al)KE =lk;K = spa r se(2* ( n e lx+1 ) *(nely+lh 2 *(nel x +1 ) * (ne 1 y+1);F = s pa r s e(2*(ne 1 y + 1 )*( nelx+1)J) ; U = zeros(2 * ( ne 1 y + 1 )*(nelx+l ) J);fore 1 x = l:(nelx - 20)for el y = 1 :nelyn i = (nely

38、 + l)*(elx-1 )+ely ;n2 = (nely+1)*elx +e 1 y ;edof= 2*nl - 1 ; 2*nl; 2 *n2-l;2*n 2 ; 2*n2+ 1:2* n2+2; 2 * n 1+1;2*nl+2;K ( edof,edof) = K(edof.edof) + x(ely,elx ) Apena 1 *K E ； ende ndfor elx = (nelx-1 9 ) :nel xfor ely =1:15n 1 =(nely+l ) * (e 1 x - l) + ely ;n2= (nel y+1 ) * elx +ely;edof =2*nl-l ; 2*nl; 2 *n2-l ； 2 * n2: 2*n2+1; 2 * n2 + 2; 2*n