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

1、華中科技大學(xué)研究生課程考試答題本基于matlab和optistruct的c形夾拓?fù)鋬?yōu)化學(xué) 院（部）： 機(jī)械科學(xué)與工程學(xué)院 課 程 名 稱： 工程優(yōu)化設(shè)計(jì) 學(xué) 生 姓 名： 班 級(jí)： 機(jī)碩1107 班 學(xué) 號(hào)： 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è)計(jì)的數(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è)計(jì)變量的類型劃分成三個(gè)層次：結(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è)計(jì)區(qū)域內(nèi)尋求最優(yōu)材料的分布問題。尋求一個(gè)最佳的拓?fù)浣Y(jié)構(gòu)形式有兩種基本的原理：一種是退化原理，另一種是進(jìn)化原理。退化原理的基本思想是在優(yōu)化前將結(jié)構(gòu)所有可能桿單元或所有材料都加上，然后構(gòu)造適當(dāng)?shù)膬?yōu)化模型，通過(guò)一定的優(yōu)化方法逐步刪減那些不必要的結(jié)構(gòu)元素，直至最終得到一個(gè)最優(yōu)化的拓?fù)浣Y(jié)構(gòu)形式。進(jìn)化原理的基本思想是把適者生存的生物進(jìn)化論思想引入結(jié)構(gòu)拓?fù)鋬?yōu)化，它通過(guò)模擬適者生存、物競(jìng)天擇、優(yōu)勝劣汰等自然機(jī)理來(lái)獲得最優(yōu)的拓?fù)浣Y(jié)構(gòu)。 1.1.2 本文

4、研究意義目前，結(jié)構(gòu)優(yōu)化大部分集中在尺寸設(shè)計(jì)變量 (如板厚、桿的剖面積及管梁的直徑 )。拓?fù)浣Y(jié)構(gòu)優(yōu)化較尺寸優(yōu)化復(fù)雜 ,但對(duì)于有些問題拓?fù)浣Y(jié)構(gòu)優(yōu)化比尺寸優(yōu)化有效 ,c形夾是其中的例子之一。本文討論c形夾的拓?fù)鋬?yōu)化問題 ,圍繞這一問題,怎樣使結(jié)構(gòu)具有最大剛度的設(shè)計(jì)占有相當(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年來(lái)從結(jié)構(gòu)優(yōu)化研究中派生出來(lái)的新分支，它在計(jì)算結(jié)構(gòu)力學(xué)中已經(jīng)被認(rèn)為是最富挑戰(zhàn)性的一類研究工作。目前有關(guān)結(jié)構(gòu)拓?fù)鋬?yōu)化的工程應(yīng)用研究還很不成熟，在國(guó)外處在發(fā)展的初期，尤其在國(guó)內(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è)計(jì)，開創(chuàng)了連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)研究的新局面。1993年xieym和stevengp提出了漸進(jìn)結(jié)構(gòu)優(yōu)化法。1999年bendsoe和sigmund證實(shí)了變密度法物理意義的存在性。2002年羅鷹等提出三角網(wǎng)格

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

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

8、型的數(shù)學(xué)表達(dá)形式： （2.1） （2.2）其中：為兩數(shù)學(xué)模型中對(duì)中間密度材料的懲罰因子。懲罰因子的作用是當(dāng)設(shè)計(jì)變量的值在(0,1)之間時(shí)，通過(guò)逐漸增加的值對(duì)設(shè)計(jì)變量的中間值進(jìn)行懲罰，隨著值的增大，設(shè)計(jì)逐漸接近0/1設(shè)計(jì)。為有效壓縮中間密度材料，要求。表示插值以后的彈性模量，和分別為固體和空洞部分材料的彈性模量，/1000。表示單元的設(shè)計(jì)變量。表示插值以后的剛度矩陣，表示第個(gè)單元固體材料的剛度矩陣。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對(duì)應(yīng)的柔度函數(shù)和敏度

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

10、動(dòng)變量在迭代過(guò)程中作為設(shè)計(jì)變量允許發(fā)生改變，；當(dāng)時(shí)，僅設(shè)計(jì)變量下限約束起作用，設(shè)計(jì)變量為被動(dòng)變量，；當(dāng)時(shí)，僅設(shè)計(jì)變量上限約束起作用，設(shè)計(jì)變量也為被動(dòng)變量，。被動(dòng)變量在迭代過(guò)程中不能變化，只能由側(cè)面約束的邊界值來(lái)確定2。第3章 優(yōu)化設(shè)計(jì)的數(shù)學(xué)模型及解決方案3.1問題描述如圖所示，c形夾在自由端口受到三角形分布力f的作用，要求在保持對(duì)原材料體積一定縮減比的情況下對(duì)原實(shí)體懸臂梁做結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)，優(yōu)化目標(biāo)是使結(jié)構(gòu)剛度最大。（優(yōu)化的結(jié)果應(yīng)該使原設(shè)計(jì)區(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é)模型該問題中，要求同時(shí)滿足剛度最大，質(zhì)量最輕，這兩個(gè)變量若同時(shí)改變，則問題復(fù)雜度太大，并且可能導(dǎo)致問題不可求解。所以我們采用在確定的質(zhì)量下，來(lái)討論剛度最大的問題。由于對(duì)特定的材料，其質(zhì)量和體積有一定的關(guān)系，并且我們采用去除法的思想來(lái)建立模型的，故我們可以采用給優(yōu)化后的體積與優(yōu)化前的體積比賦確定的值，來(lái)達(dá)到在給定質(zhì)量條件下滿足剛度最大的問題。其數(shù)學(xué)模型如下： (1)注：其中為結(jié)構(gòu)變形能，u為結(jié)構(gòu)變形總位移矩陣，k為結(jié)構(gòu)總剛度矩陣，n為劃分單元總數(shù)， 為單元位移向量，為單元?jiǎng)偠?，（由于劃分單元的時(shí)候，我們采用等分矩形單元，所以每個(gè)單元的剛度可用一個(gè)常量來(lái)處理）

12、是拓?fù)浣Y(jié)構(gòu)優(yōu)化過(guò)程中變化著的體積，為未經(jīng)過(guò)優(yōu)化前懸臂梁的體積。f為結(jié)構(gòu)所受的三角形載荷。為懸臂梁的相對(duì)密度。3.3 模型分析求解方法選擇對(duì)該問題是用兩種方法求解。方法一：基于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)等為了簡(jiǎn)化問題的復(fù)雜度，此處我們采用standard oc-method.方法來(lái)實(shí)現(xiàn)。在處理過(guò)程中，關(guān)于設(shè)計(jì)變量相對(duì)密度x每一步的更新，我們采用在1995年提出的如下算法來(lái)實(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)劃分成為三個(gè)矩形的集合，因而在整體剛度矩陣的組裝方面應(yīng)該考慮如何進(jìn)行，由于本程序是參考經(jīng)典99行oc法拓?fù)鋬?yōu)化程序改變，對(duì)于程序中避免邊界鋸齒現(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 過(guò)濾大小為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ù)迭代過(guò)程中體積比是否達(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é)果來(lái)看，在我們給定的體積保留率的情況下，每經(jīng)過(guò)一次拓?fù)浣Y(jié)構(gòu)優(yōu)化，該優(yōu)化程序就將c形夾的拓?fù)浣Y(jié)構(gòu)中強(qiáng)度要求不高處材料的密度減小，直到所有無(wú)用的材料都將被

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

18、算法的收斂速度較慢，并且還會(huì)出現(xiàn)不同程度的震蕩，并且體積壓縮率越大，該算法的振動(dòng)也震蕩。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時(shí)，拓?fù)浣Y(jié)果圖如圖5.2所示：圖5.2 c形夾拓?fù)鋬?yōu)化結(jié)果圖（volfrac=0.35）當(dāng)體積分?jǐn)?shù)volfrac=0.5時(shí)，拓?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取不同的值時(shí)，優(yōu)化得到的拓?fù)浣Y(jié)構(gòu)基本保持不變，說(shuō)明該拓?fù)浣Y(jié)構(gòu)存在一理想狀態(tài)，通過(guò)不斷嘗試，得到如下發(fā)現(xiàn)：在本例中當(dāng)volfrac=0.35時(shí)所得的拓?fù)浣Y(jié)構(gòu)最為規(guī)律，并且該拓?fù)浣Y(jié)構(gòu)趨向與桁架結(jié)構(gòu)。（具體拓?fù)浣Y(jié)構(gòu)演變過(guò)程見附件視頻動(dòng)畫）第6章 結(jié)論及總結(jié)本文中，從兩條獨(dú)立的途徑來(lái)分別對(duì)該問題進(jìn)行研究。hyperworks中，利用hypermesh自帶的的模塊對(duì)該問題進(jìn)行了建模，網(wǎng)格劃分，得到有限元模型，然后將有限元模型導(dǎo)入拓?fù)鋬?yōu)化模塊optistruct，添加邊界條件及約束，定義目標(biāo)函數(shù)

20、和約束條件，對(duì)c形夾進(jìn)行拓?fù)浣Y(jié)構(gòu)優(yōu)化。所得到的拓?fù)浣Y(jié)構(gòu)趨向于桁架，說(shuō)明桁架結(jié)構(gòu)有很高的剛度。matlab中，首先對(duì)該問題構(gòu)建數(shù)學(xué)模型，采用simp(oc)算法對(duì)數(shù)學(xué)模型進(jìn)行優(yōu)化求解，由于本文所研究的問題的求解域并不是規(guī)則的矩形，而是由三部分矩形疊加，所以在使用simp法時(shí)，對(duì)總體剛度矩陣組裝，添加邊界和約束條件成為了本文的關(guān)鍵技術(shù)，涉及到較多的有限元理論。在matlab中將每次迭代所得到的處理結(jié)果進(jìn)行了圖像動(dòng)態(tài)顯示，以此來(lái)清晰的觀察拓?fù)浣Y(jié)構(gòu)優(yōu)化的動(dòng)態(tài)過(guò)程，給人以直觀的印象。在用optistruct進(jìn)行拓?fù)浣Y(jié)構(gòu)優(yōu)化的時(shí)候，我們發(fā)現(xiàn)，當(dāng)材料去除率為60%時(shí)，其所得到的拓?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é)模型是對(duì)的。參考文獻(xiàn)1o.sigmund. a 99 line topology optimization code written in matlab.educational article, 20012羅震.基于變密度法的連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)技術(shù)研究.博士學(xué)位論文.20053孫靖民.機(jī)械優(yōu)化設(shè)計(jì)第三版.機(jī)械工業(yè)出版社.20044白新理等.結(jié)構(gòu)優(yōu)化設(shè)計(jì).黃河水利出版社.2008.45郭仁生.機(jī)械工程設(shè)計(jì)分析和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% 敏度過(guò)濾，避免拓?fù)浣Y(jié)構(gòu)邊界鋸齒化 dc = check(nelx,nely,rmin,x,dc); % 通過(guò)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、% 單元密度過(guò)濾 %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