研究生課程考試答題本基于MATLAB和Optistruct的C形夾拓?fù)鋬?yōu)化

1、華中科技大學(xué)研究生課程考試答題本基于matlab和optistruct的c形夾拓?fù)鋬?yōu)化學(xué) 院（部）： 機(jī)械科學(xué)與工程學(xué)院 課 程 名 稱： 工程優(yōu)化設(shè)計 學(xué) 生 姓 名： 班 級： 機(jī)碩1107 班 學(xué) 號： 2012年01月10日目錄第1章 選題背景介紹及問題描述31.1選題背景及意義31.1.1工程背景及基本原理31.1.2 本文研究意義31.2研究現(xiàn)狀31.2.1 理論研究現(xiàn)狀31.2.2 應(yīng)用研究現(xiàn)狀41.3 該研的意義4第2章 simp變密度法理論基礎(chǔ)52.1 simp密度剛度插值法理論基礎(chǔ)52.2 拓?fù)鋬?yōu)化的數(shù)學(xué)模型52.3 優(yōu)化準(zhǔn)則的基本理論6第3章 優(yōu)化設(shè)計的數(shù)學(xué)模型及解決方案

2、73.1問題描述73.2 優(yōu)化問題的數(shù)學(xué)模型73.3 模型分析求解方法選擇8第四章 拓?fù)鋬?yōu)化步驟及結(jié)果94.1 基于matlab的變密度拓?fù)鋬?yōu)化94.1.1 問題求解的關(guān)鍵技術(shù)及代碼94.1.2 拓?fù)鋬?yōu)化結(jié)果及分析104.2 基于optistruct的c形夾拓?fù)鋬?yōu)化114.2.1有限元模型的建立114.2.2 基于optistruct的拓?fù)鋬?yōu)化結(jié)果及分析12第6章 結(jié)論及總結(jié)14參考文獻(xiàn)15附件：懸臂梁拓?fù)鋬?yōu)化matlab程序15第1章 選題背景介紹及問題描述1.1選題背景及意義1.1.1工程背景及基本原理通常把結(jié)構(gòu)優(yōu)化按設(shè)計變量的類型劃分成三個層次：結(jié)構(gòu)尺寸優(yōu)化、形狀優(yōu)化和拓?fù)鋬?yōu)化。尺寸優(yōu)化

3、和形狀優(yōu)化已得到充分的發(fā)展，但它們存在著不能變更拓?fù)浣Y(jié)構(gòu)的缺陷。在這樣的背景下，人們開始研究拓?fù)鋬?yōu)化。拓?fù)鋬?yōu)化的基本思想是將尋求結(jié)構(gòu)的最優(yōu)拓?fù)鋯栴}轉(zhuǎn)化為在給定的設(shè)計區(qū)域內(nèi)尋求最優(yōu)材料的分布問題。尋求一個最佳的拓?fù)浣Y(jié)構(gòu)形式有兩種基本的原理：一種是退化原理，另一種是進(jìn)化原理。退化原理的基本思想是在優(yōu)化前將結(jié)構(gòu)所有可能桿單元或所有材料都加上，然后構(gòu)造適當(dāng)?shù)膬?yōu)化模型，通過一定的優(yōu)化方法逐步刪減那些不必要的結(jié)構(gòu)元素，直至最終得到一個最優(yōu)化的拓?fù)浣Y(jié)構(gòu)形式。進(jìn)化原理的基本思想是把適者生存的生物進(jìn)化論思想引入結(jié)構(gòu)拓?fù)鋬?yōu)化，它通過模擬適者生存、物競天擇、優(yōu)勝劣汰等自然機(jī)理來獲得最優(yōu)的拓?fù)浣Y(jié)構(gòu)。 1.1.2 本文

4、研究意義目前，結(jié)構(gòu)優(yōu)化大部分集中在尺寸設(shè)計變量 (如板厚、桿的剖面積及管梁的直徑 )。拓?fù)浣Y(jié)構(gòu)優(yōu)化較尺寸優(yōu)化復(fù)雜 ,但對于有些問題拓?fù)浣Y(jié)構(gòu)優(yōu)化比尺寸優(yōu)化有效 ,c形夾是其中的例子之一。本文討論c形夾的拓?fù)鋬?yōu)化問題 ,圍繞這一問題,怎樣使結(jié)構(gòu)具有最大剛度的設(shè)計占有相當(dāng)重要的地位;怎樣優(yōu)化結(jié)構(gòu)的形狀使材料的分布，更加合理從而達(dá)到使結(jié)構(gòu)具有最小柔度的目的是本文要研究的問題。1.2研究現(xiàn)狀1.2.1 理論研究現(xiàn)狀結(jié)構(gòu)拓?fù)鋬?yōu)化是近20年來從結(jié)構(gòu)優(yōu)化研究中派生出來的新分支，它在計算結(jié)構(gòu)力學(xué)中已經(jīng)被認(rèn)為是最富挑戰(zhàn)性的一類研究工作。目前有關(guān)結(jié)構(gòu)拓?fù)鋬?yōu)化的工程應(yīng)用研究還很不成熟，在國外處在發(fā)展的初期，尤其在國內(nèi)

5、尚屬于起步階段。1904年michell在桁架理論中首次提出了拓?fù)鋬?yōu)化的概念。自1964年dorn等人提出基結(jié)構(gòu)法，將數(shù)值方法引入拓?fù)鋬?yōu)化領(lǐng)域，拓?fù)鋬?yōu)化研究開始活躍。20世紀(jì)80年代初，程耿東和n.olhoff在彈性板的最優(yōu)厚度分布研究中首次將最優(yōu)拓?fù)鋯栴}轉(zhuǎn)化為尺寸優(yōu)化問題，他們開創(chuàng)性的工作引起了眾多學(xué)者的研究興趣。1988年bendsoe和kikuchi發(fā)表的基于均勻化理論的結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計，開創(chuàng)了連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計研究的新局面。1993年xieym和stevengp提出了漸進(jìn)結(jié)構(gòu)優(yōu)化法。1999年bendsoe和sigmund證實(shí)了變密度法物理意義的存在性。2002年羅鷹等提出三角網(wǎng)格

6、進(jìn)化法，該方法在優(yōu)化過程中實(shí)現(xiàn)了退化和進(jìn)化的統(tǒng)一，提高了優(yōu)化效率。1.2.2 應(yīng)用研究現(xiàn)狀在前人提出的重要理論基礎(chǔ)上，后人也將其跟其他現(xiàn)代設(shè)計的方法相結(jié)合，衍生出了其他一些拓?fù)浣Y(jié)構(gòu)優(yōu)化方法：如與可靠性相結(jié)合的情況下，maute等應(yīng)用變密度法并結(jié)合可靠性分析對一微機(jī)電系統(tǒng)進(jìn)行了基于可靠性的拓?fù)鋬?yōu)化設(shè)計，papadrakakis等將遺傳算法應(yīng)用于具有可靠性約束的桁架結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計中，國內(nèi)學(xué)者馬洪波也對基于遺傳算法的結(jié)構(gòu)可靠性優(yōu)化問題進(jìn)行了討論。華南理工大學(xué)機(jī)械工程學(xué)院歐陽高飛等對基于水平集方法的結(jié)構(gòu)可靠性拓?fù)鋬?yōu)化進(jìn)行了研究。1.3 本文研究的意義通過這次的作業(yè)加深對工程優(yōu)化算法的學(xué)習(xí)和使用，提高

7、對拓?fù)鋬?yōu)化的方法和過程的了解和學(xué)習(xí)。另外對相關(guān)軟件軟件的應(yīng)用能夠達(dá)到一個新的高度。這些不僅能使我們現(xiàn)在的知識體系得到充實(shí)和優(yōu)化，而且也是我們今后人生的財富。第2章 simp變密度法理論基礎(chǔ)2.1 simp密度剛度插值法理論基礎(chǔ)simp模型主要通過引入懲罰因子，在材料的彈性模量和單元相對密度之間建立起一種顯示的非線性對應(yīng)關(guān)系。它的作用是當(dāng)設(shè)計變量的值在(0,1)之間時，對中間密度值進(jìn)行懲罰，使中間密度值逐漸向0/1兩端聚集，這樣可以使連續(xù)變量的拓?fù)鋬?yōu)化模型能很好地逼近原來0-1離散變量的優(yōu)化模型。這時中間密度單元對應(yīng)一個很小的彈性模量，對結(jié)構(gòu)剛度矩陣的影響將變得很小，可以忽略不計。simp材料模

8、型的數(shù)學(xué)表達(dá)形式： （2.1） （2.2）其中：為兩數(shù)學(xué)模型中對中間密度材料的懲罰因子。懲罰因子的作用是當(dāng)設(shè)計變量的值在(0,1)之間時，通過逐漸增加的值對設(shè)計變量的中間值進(jìn)行懲罰，隨著值的增大，設(shè)計逐漸接近0/1設(shè)計。為有效壓縮中間密度材料，要求。表示插值以后的彈性模量，和分別為固體和空洞部分材料的彈性模量，/1000。表示單元的設(shè)計變量。表示插值以后的剛度矩陣，表示第個單元固體材料的剛度矩陣。2.2 拓?fù)鋬?yōu)化的數(shù)學(xué)模型以結(jié)構(gòu)的柔度最小化(或剛度最大化、應(yīng)變能最小化)作為優(yōu)化的目標(biāo)函數(shù)，以結(jié)構(gòu)整體的體積約束作為優(yōu)化的約束條件。剛度優(yōu)化的數(shù)學(xué)模型表示為： （2.3）simp對應(yīng)的柔度函數(shù)和敏度

9、形式： （2.4） （2.5）其中：以上各式中，表示第個單元的剛度矩陣。表示結(jié)構(gòu)的位移向量；表示設(shè)計變量，為避免總剛度矩陣奇異，取=0.001。為單元數(shù)目，表示結(jié)構(gòu)的柔度，表示柔度關(guān)于設(shè)計變量的敏度。2.3 優(yōu)化準(zhǔn)則的基本理論剛度拓?fù)鋬?yōu)化問題，是典型的具有不等式約束的非線性規(guī)劃問題。不等式約束多元函數(shù)極值的必要條件是kuhn-tucker條件，它是采用優(yōu)化準(zhǔn)則法求解非線性優(yōu)化問題的重要理論。引入對設(shè)計變量上下限約束的拉格朗日乘子、以及對體積約束的拉格朗日乘子，構(gòu)造lagrange函數(shù)為如下形式：(2.6)對于lagrange函數(shù)，當(dāng)時，設(shè)計變量的上下限側(cè)面約束均不起作用，設(shè)計變量是主動變量。主

10、動變量在迭代過程中作為設(shè)計變量允許發(fā)生改變，；當(dāng)時，僅設(shè)計變量下限約束起作用，設(shè)計變量為被動變量，；當(dāng)時，僅設(shè)計變量上限約束起作用，設(shè)計變量也為被動變量，。被動變量在迭代過程中不能變化，只能由側(cè)面約束的邊界值來確定2。第3章 優(yōu)化設(shè)計的數(shù)學(xué)模型及解決方案3.1問題描述如圖所示，c形夾在自由端口受到三角形分布力f的作用，要求在保持對原材料體積一定縮減比的情況下對原實(shí)體懸臂梁做結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計，優(yōu)化目標(biāo)是使結(jié)構(gòu)剛度最大。（優(yōu)化的結(jié)果應(yīng)該使原設(shè)計區(qū)域產(chǎn)生孔洞，使結(jié)構(gòu)拓?fù)浒l(fā)生變化。）原實(shí)體c形夾為如下圖3.1所示的c形，尺寸如下圖所示，材料為45鋼，密度為，彈性模量，泊松比。圖3.1 c形夾的尺寸和受

11、力（單位：mm）3.2 優(yōu)化問題的數(shù)學(xué)模型該問題中，要求同時滿足剛度最大，質(zhì)量最輕，這兩個變量若同時改變，則問題復(fù)雜度太大，并且可能導(dǎo)致問題不可求解。所以我們采用在確定的質(zhì)量下，來討論剛度最大的問題。由于對特定的材料，其質(zhì)量和體積有一定的關(guān)系，并且我們采用去除法的思想來建立模型的，故我們可以采用給優(yōu)化后的體積與優(yōu)化前的體積比賦確定的值，來達(dá)到在給定質(zhì)量條件下滿足剛度最大的問題。其數(shù)學(xué)模型如下： (1)注：其中為結(jié)構(gòu)變形能，u為結(jié)構(gòu)變形總位移矩陣，k為結(jié)構(gòu)總剛度矩陣，n為劃分單元總數(shù)， 為單元位移向量，為單元剛度，（由于劃分單元的時候，我們采用等分矩形單元，所以每個單元的剛度可用一個常量來處理）

12、是拓?fù)浣Y(jié)構(gòu)優(yōu)化過程中變化著的體積，為未經(jīng)過優(yōu)化前懸臂梁的體積。f為結(jié)構(gòu)所受的三角形載荷。為懸臂梁的相對密度。3.3 模型分析求解方法選擇對該問題是用兩種方法求解。方法一：基于matlab的變密度拓?fù)鋬?yōu)化法；方法二：是用成熟的有限元拓?fù)鋬?yōu)化軟件optistruct進(jìn)行優(yōu)化?；趍atlab的變密度拓?fù)鋬?yōu)化該問題的優(yōu)化方法有很多種，常用的有如下方法：optimality criteria(oc) methods,（優(yōu)化準(zhǔn)則方法）sequential linear programming (slp) methods（序列線性規(guī)劃法）method of moving asymptotes (mma b

13、ysvanberg 1987)等為了簡化問題的復(fù)雜度，此處我們采用standard oc-method.方法來實(shí)現(xiàn)。在處理過程中，關(guān)于設(shè)計變量相對密度x每一步的更新，我們采用在1995年提出的如下算法來實(shí)現(xiàn)。第四章 拓?fù)鋬?yōu)化步驟及結(jié)果4.1 基于matlab的變密度拓?fù)鋬?yōu)化4.1.1 問題求解的關(guān)鍵技術(shù)及代碼一般而言，由于oc法所使用的單元是矩形，所以oc法很適合求解求解域?yàn)榫匦蔚膬?yōu)化問題，而本文選題為一c形結(jié)構(gòu)，若使用oc法，則需將c形結(jié)構(gòu)劃分成為三個矩形的集合，因而在整體剛度矩陣的組裝方面應(yīng)該考慮如何進(jìn)行，由于本程序是參考經(jīng)典99行oc法拓?fù)鋬?yōu)化程序改變，對于程序中避免邊界鋸齒現(xiàn)象所使用的

14、check函數(shù)，該如何進(jìn)行修改，以及如何進(jìn)行邊界條件的施加。優(yōu)化問題的初值條件：nelx=60 x方向單元的數(shù)目為60nely=50 y方向單元的數(shù)目為50volfrac=0.35 保留原材料的體積分?jǐn)?shù)為0.35penal=3.0 抑制權(quán)值為3.0（該取值是資料建議的典型值）rmin=1.2 過濾大小為1.2（該取值是資料建議的典型值）關(guān)鍵代碼：function xnew=oc(nelx,nely,x,volfrac,dc) %定義oc優(yōu)化準(zhǔn)則函數(shù)l1 = 0; l2 = 100000; move = 0.2;while (l2-l1 1e-4) lmid = 0.5*(l2+l1); %采用

15、折半查找 xnew = max(0.001,max(x-move,min(1.,min(x+move,x.*sqrt(-dc./lmid);%根據(jù)和敏度dc更新x的值% if sum(sum(xnew) - volfrac*nelx*nely 0;%根據(jù)迭代過程中體積比是否達(dá)到預(yù)設(shè)的體積比，判斷迭代是否繼續(xù)進(jìn)行% l1 = lmid; else l2 = lmid; endend for ely = 1:nely for elx = 1:(nelx-20) %取c形夾的左半矩形區(qū)域 n1 = (nely+1)*(elx-1)+ely; n2 = (nely+1)* elx +ely; ue =

16、 u(2*n1-1;2*n1; 2*n2-1;2*n2; 2*n2+1;2*n2+2; 2*n1+1;2*n1+2,1);%單元位移向量% x(16:35,41:nelx)=0.001; c = c + x(ely,elx)penal*ue*ke*ue;%目標(biāo)函數(shù) dc(ely,elx) = -penal*x(ely,elx)(penal-1)*ue*ke*ue;%敏度值dc end end4.1.2 拓?fù)鋬?yōu)化結(jié)果及分析 實(shí)驗(yàn)結(jié)果分析：從該實(shí)驗(yàn)結(jié)果來看，在我們給定的體積保留率的情況下，每經(jīng)過一次拓?fù)浣Y(jié)構(gòu)優(yōu)化，該優(yōu)化程序就將c形夾的拓?fù)浣Y(jié)構(gòu)中強(qiáng)度要求不高處材料的密度減小，直到所有無用的材料都將被

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

18、算法的收斂速度較慢，并且還會出現(xiàn)不同程度的震蕩，并且體積壓縮率越大，該算法的振動也震蕩。4.2 基于optistruct的c形夾拓?fù)鋬?yōu)化4.2.1有限元模型的建立在optistruct中建立有限元模型并劃分網(wǎng)格。相關(guān)參數(shù)為：;添加邊界調(diào)節(jié)和載荷，得到有限元模型如下圖5.1所示：圖5.1 有限元模型4.2.2 基于optistruct的拓?fù)鋬?yōu)化結(jié)果及分析當(dāng)體積分?jǐn)?shù)volfrac=0.35時，拓?fù)浣Y(jié)果圖如圖5.2所示：圖5.2 c形夾拓?fù)鋬?yōu)化結(jié)果圖（volfrac=0.35）當(dāng)體積分?jǐn)?shù)volfrac=0.5時，拓?fù)浣Y(jié)果圖如圖5.3所示：圖5.3 c形夾拓?fù)鋬?yōu)化結(jié)果圖（volfrac=0.5）拓?fù)鋬?yōu)

19、化結(jié)構(gòu)分析：在以變形最小為目標(biāo)函數(shù)，即結(jié)構(gòu)剛度最大，一定的材料保留比率volfrac為約束調(diào)節(jié)，當(dāng)volfrac取不同的值時，優(yōu)化得到的拓?fù)浣Y(jié)構(gòu)基本保持不變，說明該拓?fù)浣Y(jié)構(gòu)存在一理想狀態(tài)，通過不斷嘗試，得到如下發(fā)現(xiàn)：在本例中當(dāng)volfrac=0.35時所得的拓?fù)浣Y(jié)構(gòu)最為規(guī)律，并且該拓?fù)浣Y(jié)構(gòu)趨向與桁架結(jié)構(gòu)。（具體拓?fù)浣Y(jié)構(gòu)演變過程見附件視頻動畫）第6章 結(jié)論及總結(jié)本文中，從兩條獨(dú)立的途徑來分別對該問題進(jìn)行研究。hyperworks中，利用hypermesh自帶的的模塊對該問題進(jìn)行了建模，網(wǎng)格劃分，得到有限元模型，然后將有限元模型導(dǎo)入拓?fù)鋬?yōu)化模塊optistruct，添加邊界條件及約束，定義目標(biāo)函數(shù)

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

21、行拓?fù)浣Y(jié)構(gòu)優(yōu)化是所得到的結(jié)果的拓?fù)浣Y(jié)構(gòu)是一致的。這驗(yàn)證了我們的matlab數(shù)學(xué)模型是對的。參考文獻(xiàn)1o.sigmund. a 99 line topology optimization code written in matlab.educational article, 20012羅震.基于變密度法的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計技術(shù)研究.博士學(xué)位論文.20053孫靖民.機(jī)械優(yōu)化設(shè)計第三版.機(jī)械工業(yè)出版社.20044白新理等.結(jié)構(gòu)優(yōu)化設(shè)計.黃河水利出版社.2008.45郭仁生.機(jī)械工程設(shè)計分析和matlab應(yīng)用.機(jī)械工業(yè)出版社.2006附件：c形夾拓?fù)鋬?yōu)化matlab程序function top(n

22、elx,nely,volfrac,penal,rmin);%變量初始化nelx=60;nely=50;volfrac=0.35;penal=3;rmin=2;x(1:nely,1:nelx) = volfrac; loop = 0; change = 1.;% 開始迭代while change 0.01 loop = loop + 1; xold = x; u=fe(nelx,nely,x,penal); %調(diào)用子函數(shù)得到整體位移向量 % 將求解域分成三塊矩形進(jìn)行目標(biāo)函數(shù)定義和敏度分析 ke = lk; c = 0.; for ely = 1:nely for elx = 1:(nelx-20

23、) n1 = (nely+1)*(elx-1)+ely; n2 = (nely+1)* elx +ely; ue = u(2*n1-1;2*n1; 2*n2-1;2*n2; 2*n2+1;2*n2+2; 2*n1+1;2*n1+2,1); x(16:35,41:nelx)=0.001; c = c + x(ely,elx)penal*ue*ke*ue;%目標(biāo)函數(shù) dc(ely,elx) = -penal*x(ely,elx)(penal-1)*ue*ke*ue;%敏度 end end for ely = 1:15 for elx = 41:nelx n1 = (nely+1)*(elx-1)+

24、ely; n2 = (nely+1)* elx +ely; ue = u(2*n1-1;2*n1; 2*n2-1;2*n2; 2*n2+1;2*n2+2; 2*n1+1;2*n1+2,1); x(16:35,41:nelx)=0.001; c = c + x(ely,elx)penal*ue*ke*ue;%目標(biāo)函數(shù) dc(ely,elx) = -penal*x(ely,elx)(penal-1)*ue*ke*ue; end end for ely = 36:nely for elx = 41:nelx n1 = (nely+1)*(elx-1)+ely; n2 = (nely+1)* elx

25、+ely; ue = u(2*n1-1;2*n1; 2*n2-1;2*n2; 2*n2+1;2*n2+2; 2*n1+1;2*n1+2,1); x(16:35,41:nelx)=0.001; c = c + x(ely,elx)penal*ue*ke*ue;%目標(biāo)函數(shù) dc(ely,elx) = -penal*x(ely,elx)(penal-1)*ue*ke*ue; end end% 敏度過濾，避免拓?fù)浣Y(jié)構(gòu)邊界鋸齒化 dc = check(nelx,nely,rmin,x,dc); % 通過oc法更新單元密度 x = oc(nelx,nely,x,volfrac,dc); for j =(n

26、elx-19):(nelx-1) %指定非優(yōu)化區(qū)域 x(15,j)=0.9; end for j =(nelx-19):(nelx-1) x(36,j)=0.9; end change = max(max(abs(x-xold); disp( it.: sprintf(%4i,loop) obj.: sprintf(%10.4f,c) . vol.: sprintf(%6.3f,sum(sum(x)/(nelx*nely) . ch.: sprintf(%6.3f,change )%顯示結(jié)果 colormap(gray); imagesc(-x); axis equal; axis tight

27、; axis off;pause(1e-6);%繪制密度圖end % 最優(yōu)性判據(jù)更新單元密度值 %function xnew=oc(nelx,nely,x,volfrac,dc) l1 = 0; l2 = 100000; move = 0.2;while (l2-l1 1e-4) lmid = 0.5*(l2+l1); xnew = max(0.001,max(x-move,min(1.,min(x+move,x.*sqrt(-dc./lmid); if sum(sum(xnew) - volfrac*nelx*nely 0; l1 = lmid; else l2 = lmid; endend

28、% 單元密度過濾 %function dcn=check(nelx,nely,rmin,x,dc)dcn=zeros(nely,nelx);for i = 1:(nelx-20) for j = 1:nely sum=0.0; for k = max(i-floor(rmin),1):min(i+floor(rmin),nelx) for l = max(j-floor(rmin),1):min(j+floor(rmin),nely) fac = rmin-sqrt(i-k)2+(j-l)2); sum = sum+max(0,fac); dcn(j,i) = dcn(j,i) + max(0

29、,fac)*x(l,k)*dc(l,k); end end dcn(j,i) = dcn(j,i)/(x(j,i)*sum); endendfor i = (nelx-19):nelx for j = 1:15 sum=0.0; for k = max(i-floor(rmin),1):min(i+floor(rmin),nelx) for l = max(j-floor(rmin),1):min(j+floor(rmin),nely) fac = rmin-sqrt(i-k)2+(j-l)2); sum = sum+max(0,fac); dcn(j,i) = dcn(j,i) + max(

30、0,fac)*x(l,k)*dc(l,k); end end dcn(j,i) = dcn(j,i)/(x(j,i)*sum); endendfor i = (nelx-19):nelx for j = (nely-14):nely sum=0.0; for k = max(i-floor(rmin),1):min(i+floor(rmin),nelx) for l = max(j-floor(rmin),1):min(j+floor(rmin),nely) fac = rmin-sqrt(i-k)2+(j-l)2); sum = sum+max(0,fac); dcn(j,i) = dcn(

31、j,i) + max(0,fac)*x(l,k)*dc(l,k); end end dcn(j,i) = dcn(j,i)/(x(j,i)*sum); endend% 有限元分析 %function u=fe(nelx,nely,x,penal)ke = lk; k = sparse(2*(nelx+1)*(nely+1), 2*(nelx+1)*(nely+1);f = sparse(2*(nely+1)*(nelx+1),1); u = zeros(2*(nely+1)*(nelx+1),1);for elx = 1:(nelx-20) for ely = 1:nely n1 = (nel

32、y+1)*(elx-1)+ely; n2 = (nely+1)* elx +ely; edof = 2*n1-1; 2*n1; 2*n2-1; 2*n2; 2*n2+1; 2*n2+2; 2*n1+1; 2*n1+2; k(edof,edof) = k(edof,edof) + x(ely,elx)penal*ke; endendfor elx = (nelx-19):nelx for ely = 1:15 n1 = (nely+1)*(elx-1)+ely; n2 = (nely+1)* elx +ely; edof = 2*n1-1; 2*n1; 2*n2-1; 2*n2; 2*n2+1; 2*n2+2; 2*n1+1; 2*n1+2; k(edof,edof) = k(edof,edof) + x(e