基于MATLAB和Optistruct的C形夹拓扑优化

1、华中科技大学研究生课程考试答题本谨于MATLAB和Optistruct的C形夹拓扑优化学院（部）：机械科学与工程学院课程名称： 工程优化设计学生姓名：范利洪班 级：机硕11 07班学M201170602号：2012年01月1 0日目录第1章选题背景介绍及问题描述。错误!未定义书签。1 .1选题背景及意义。错误!未定义书签。1.1.1 工程背景及基本原理。错误!未定义书签。1.1.2 本文研究意义。错误!未定义书签。1.2 研究现状错误!未定义书签。1.2. 1理论研究现状错误!未定义书签。1.3. 2应用研究现状。错误!未定义书签。1.3 该研的意义4第2章SIMP变密度法理论基础错误!未定义

2、书签。2. 1SIM P密度刚度插值法理论基础错误!未定义书签。2.1 拓扑优化的数学模型。错误!未定义书签。2.2 优化准则的基本理论。错误!未定义书签。 第3章优化设计的数学模型及解决方案。错误!未定义书签。3. 1问题描述错误!未定义书签。3.1 优化问题的数学模型错误!未定义书签。3.2 模型分析求解方法选择错误!未定义书签。第四章拓扑优化步骤及结果错误!未定义书签。4. 1基于matlab的变密度拓扑优化。错误!未定义书签。4.1 1问题求解的关键技术及代码。错误!未定义书签。4. 1.2拓扑优化结果及分析错误!未定义书签。4.2基于Optistruct的C形夹拓扑优化。错误!未定义

3、书签。1 .2.1有限元模型的建立。错误!未定义书签。4 . 2 .2基于Opti s truct的拓扑优化结果及分析。错误!未定义书签。第6章 结论及总结错误!未定义书签。参考文献错误!未定义书签。附件：悬臂梁拓扑优化mat 1ab程序错误!未定义书签。第1章 选题背景介绍及问题描述1.1 选题背景及意义1.1.1 工程背景及基本原理通常把结构优化按设计变量的类型划分成三个层次:结构尺寸优化、形状优 化和拓扑优化。尺寸优化和形状优化已得到充分的发展，但它们存在着不能变更 拓扑结构的缺陷。在这样的背景下，人们开始研究拓扑优化。拓扑优化的基本思 想是将寻求结构的最优拓扑问题转化为在给定的设计区域

4、内寻求最优材料的分 布问题。寻求一个最佳的拓扑结构形式有两种基本的原理:一种是退化原理，另一 种是进化原理。退化原理的基本思想是在优化前将结构所有可能杆单元或所有材 料都加上，然后构造适当的优化模型,通过一定的优化方法逐步删减那些不必要 的结构元素，直至最终得到一个最优化的拓扑结构形式。进化原理的基本思想是 把适者生存的生物进化论思想引入结构拓扑优化，它通过模拟适者生存、物竞天 择、优胜劣汰等自然机理来获得最优的拓扑结构。1.1.2 本文研究意义目前，结构优化大部分集中在尺寸设计变量（如板厚、杆的剖面积及管梁的 直径）。拓扑结构优化较尺寸优化复杂，但对于有些问题拓扑结构优化比尺寸 优化有效，c

5、形夹是其中的例子之一。本文讨论c形夹的拓扑优化问题，围绕这 一问题，怎样使结构具有最大刚度的设计占有相当重要的地位;怎样优化结构的形 状使材料的分布,更加合理从而达到使结构具有最小柔度的目的是本文要研究的 问题。1-2研究现状1.2 . 1理论研究现状结构拓扑优化是近20年来从结构优化研究中派生出来的新分支，它在计算 结构力学中已经被认为是最富挑战性的一类研究工作。目前有关结构拓扑优化的 工程应用研究还很不成熟，在国外处在发展的初期，尤其在国内尚属于起步阶段。 1 904年M i chell在桁架理论中首次提出了拓扑优化的概念。自1964年Dorn 等人提出基结构法，将数值方法引入拓扑优化领域

6、，拓扑优化研究开始活跃。20世 纪8 0年代初，程耿东和N.O 1 ho f f在弹性板的最优厚度分布研究中首次将最 优拓扑问题转化为尺寸优化问题，他们开创性的工作引起了众多学者的研究兴 趣。1 9 88年Be n d so e和Ki ku c h i发表的基于均匀化理论的结构拓扑优化设 计，开创了连续体结构拓扑优化设计研究的新局面。I 9 93年Xi e YM和S t evenGP提出了渐进结构优化法。1999年Bendsoe和S i gmund证实了变密度法 物理意义的存在性。2002年罗鹰等提出三角网格进化法，该方法在优化过程中实 现了退化和进化的统一，提高了优化效率。1.3 . 2应用

7、研究现状在前人提出的重要理论基础上，后人也将其跟其他现代设计的方法相结合， 衍生出了其他一些拓扑结构优化方法：如与可靠性相结合的情况下，MAUTE 等应用变密度法并结合可靠性分析对一微机电系统进行了基于可靠性的拓扑优 化设计,PAPADRA KAKIS等将遗传算法应用于具有可靠性约束的桁架结构拓 扑优化设计中，国内学者马洪波也对基于遗传算法的结构可靠性优化问题进行了 讨论。华南理工大学机械工程学院欧阳高飞等对基于水平集方法的结构可靠性拓 扑优化进行了研究。1.4 本文研究的意义通过这次的作业加深对工程优化算法的学习和使用，提高对拓扑优化的方法 和过程的了解和学习。另外对相关软件软件的应用能够达

8、到一个新的高度。这些 不仅能使我们现在的知识体系得到充实和优化，而且也是我们今后人生的财富。第2章 SIMP变密度法理论基础2.1 SIMP密度刚度插值法理论基础S I M P模型主要通过引入惩罚因子，在材料的弹性模量和单元相对密度之 间建立起一种显示的非线性对应关系。它的作用是当设计变量的值在(0, 1 )之 间时.，对中间密度值进行惩罚,使中间密度值逐渐向0/1两端聚集，这样可以使 连续变量的拓扑优化模型能很好地逼近原来01离散变量的优化模型。这时中 间密度单元对应一个很小的弹性模量，对结构刚度矩阵的影响将变得很小，可以 忽略不计。SIMP材料模型的数学表达形式：¥但)二炉加+必

9、(£° E巧JJ1) .1)的=之(即+叼3储(2.2 )其中:p为两数学模型中对中间密度材料的惩罚因子。惩罚因子的作用是当 设计变量的值在(0, 1)之间时，通过逐渐增加的值对设计变量的中间值进行惩 罚，随着值的增大，设计逐渐接近0/1设计。为有效压缩中间密度材料，要 求 2 2。EP表示插值以后的弹性模量,石。和Em,n分别为固体和空洞部分材料的 弹性模量，AE = E°-Em,n,y'-5。/ ooo。,(/ = 1,2,表示单元/的设计 变量。K表示插值以后的刚度矩阵,K/表示笫，个单元固体材料的刚度矩阵。2) 2拓扑优化的数学模型以结构的柔度最小

10、化(或刚度最大化、应变能最小化)作为优化的目标函数， 以结构整体的体积约束作为优化的约束条件。刚度优化的数学模型表示为：'Minimize :C(X) = UrKUX = $ 也,./ P 龄广小Subject to J 0 < x而门-xj - 1，KU=F.(2.3)S I MP对应的柔度函数和敏度形式：c(x) = X(£min+ 考 AE)U%j=i(2.4)Cx) = £ px/TUAEj=i(2. 5)其中:以上各式中，KJ表示第j个单元的刚度矩阵。U表示结构的位移向 量；X表示设计变量，为避免总刚度矩阵奇异，取Xmin=0-00 1。为单元数目，

11、 。表示结构的柔度,C'表示柔度关于设计变量的敏度。2.3优化准则的基本理论刚度拓扑优化问题，是典型的具有不等式约束的非线性规划问题。不等式约 束多元函数极值的必要条件是Kuhn-Tucker条件，它是采用优化准则法求解 非线性优化问题的重要理论。引入对设计变量上下限约束的拉格朗日乘子力、万 以及对体积约束的拉格朗日乘子A ,构造Lagrange函数L为如下形式： 乙=。(%)+ 八43厂 0) +，(KU-用)+ £用(-七)+ 支片区72) j=ij=j=i(2.6)对于La grange函数，当x)时,设计变量的上下限侧面约束均不起 作用，设计变量是主动变量。主动变量在

12、迭代过程中作为设计变量允许发生改 变,右=片=0;当/加之与时，仅设计变量下限约束起作用，设计变量为被动变 量,*20,片=0；当心a、",时，仅设计变量上限约束起作用，设计变量也为被动 变量，2>0,2；>0o被动变量在迭代过程中不能变化，只能由侧面约束的边界 值来确定。第3章优化设计的数学模型及解决方案31问题描述如图所示，C形夹在自由端口受到三角形分布力F的作用，要求在保持对原 材料体积一定缩减比的情况下对原实体悬臂梁做结构拓扑优化设计，优化目标是 使结构刚度最大。（优化的结果应该使原设计区域产生孔洞，使结构拓扑发生变 化。）原实体C形夹为如下图3. 1所示的C形，

13、尺寸如下图所示，材料为4 5钢, 密度p为7.8x/加3，弹性模量E = 2e5MPa，泊松比 =0.3。F = 20N。3.2优化问题的数学模型该问题中，要求同时满足刚度最大，质量最轻，这两个变量若同时改变，则问 题复杂度太大，并且可能导致问题不可求解。所以我们采用在确定的质量下，来讨 论刚度最大的问题。由于对特定的材料，其质量和体积有一定的关系，并且我们 采用去除法的思想来建立模型的，故我们可以采用给优化后的体积与优化前的体 积比赋确定的值，来达到在给定质量条件下满足刚度最大的问题。其数学模型如 下：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)为结构变形能,U为结构变形总位移矩阵，K为结构总刚度矩 阵,N为划分单元总数,% 为单元位移向量,。为单元刚度，(由于划分单元的时 候，我们采用等分矩形单元，所以每个单元的刚度可用一个常量来处理)匕是拓 扑结构优化过程中变化着的体积,0为未经过优化前悬臂梁的体积。F为结构所 受的三角形载荷。，为悬臂梁的相对密度。3.3模型分析求解方法选择对该问题是用两种方法求解。方法一:基于mat lab的变密度拓扑优化法; 方法二:是用成熟

15、的有限元拓扑优化软件Optis t ruct进行优化。基于matlab的变密度拓扑优化该问题的优化方法有很多种，常用的有如下方法：Optima 1 it y Crite r ia (OC) methods,(优化准则方法)Sequentia 1 Linear P rogramm i ng (SLP ) meth o ds (序列线性规划法) Meth o d of Moving Asymp t o tes (MM A by S van ber g 1987)等 为了简化问题的复杂度,此处我们采用stand a rd 0 C-meth。d .方法来 实现。在处理过程中，关于设计变量相对密度x每一

16、步的更新，我们采用灰公四 在1 99 5年提出的如下算法来实现。第四章拓扑优化步骤及结果4.1基于ma t la b的变密度拓扑优化4.1. 1问题求解的关键技术及代码一般而言，由于OC法所使用的单元是矩形，所以0C法很适合求解求解域 为矩形的优化问题，而本文选题为一 C形结构，若使用OC法，则需将C形结 构划分成为三个矩形的集合，因而在整体刚度矩阵的组装方面应该考虑如何进 行，由于本程序是参考经典9 9行OC法拓扑优化程序改变，对于程序中避免边界 锯齿现象所使用的check函数，该如何进行修改，以及如何进行边界条件的施加。优化问题的初值条件：Nelx=6 0 x方向单元的数目为60ne 1

17、y=50y方向单元的数目为50vo 1 f rac = 0.35保留原材料的体积分数为0 . 35p e nal=3. 0抑制权值为3 .0 (该取值是资料建议的典型值)rmin=l. 2过滤大小为1. 2 (该取值是资料建议的典型值)关键代码：funct i o n xnew=OC(n e lx, nely,x, volfrac,de)%定义 OC 优化准则函数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)；%根据Xj和敏度de更新X的值%if s u m(sum(xnew) - volfrac*n e 1 x * n e ly > 0 ;%根据迭代过程 中体积比是否达到预设的体积比，判断迭代是否继续进行%11 = Imid;else1 2 = Im i d;endendfor e 1 y = 1 :nelyfor elx = I :( n e 1 x-20)%取C形夹的左半矩形区域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;% 目标函数 dc(el y , e 1 x ) = - p enal* x (eljse 1 x)A( penal-1 )*Ue* KE*Ue;% 敏度值deendend4.1.2拓扑优化结果及分析实验结果分析：从该实验结果

20、来看，在我们给定的体积保留率的情况下,每经过一次拓扑结构 优化，该优化程序就将C形夹的拓扑结构中强度要求不高处材料的密度减小，直 到所有无用的材料都将被去除为止。我们的拓扑结构优化模型是建立在结构变形 能最小、体积去除率自己给定的基础上进行的，故我们可以根据实际情况，自行确 定体积去除率。在拓扑优化的过程中，我们可以观察到，我无论体积压缩率如何变化,c形 夹模型最终都向桁架结构进化。这说明，在结构件中，在自身材料多少相同的条 件下,桁架具有很高的刚度和强度，其实这也就是为什么拓扑结构优化首先在桁架 结构领域提出。故工程上，我们常见工程人员采用桁架结构来作为一些工程的支 撑结构，如塔吊等。在验证

21、不同的体积压缩率时，在不同的给定体积压缩率下，算法的有效性也 不同，但在验证过程中，发现算法一直会收敛。体积压缩率小的时候，该算法能 很快终止;体积压缩率较大的时候，该算法的收敛速度较慢,并且还会出现不同程 度的震荡，并且体积压缩率越大，该算法的振动也震荡。4.2基于Optis truct的C形夹拓扑优化4.2.1 有限元模型的建立在Optis true t中建立有限元模型并划分网格。相关参数为：E=2xl05：4=03；密度片7800口/ m3添加边界调节和载荷，得到有限元模型如下图5. 1所示：图5.1有限元模型4. 2.2基于O p t i s t met的拓扑优化结果及分析当体积分数v

22、ol frac =0. 35时，拓扑结果图如图5. 2所示:Rssult E /tmp/CMAX/cxingjial-des.KJdDesign Iteration 30Fram31图5.2 C形夹拓扑优化结果图(vol f r a c=0.3 5 ) 当体积分数v olf r ac= 0 . 5时，拓扑结果图如图5 .3所示：Result E./tmp/CMAX/ cxinqj ia1 _dc5. h3dDesign ： Iteration 30LFrame 31图5.3 C形夹拓扑优化结果图(volfrac= 0 .5)拓扑优化结构分析：在以变形最小为目标函数，即结构刚度最大，一定的材 料

23、保留比率volf r a c为约束调节，当v o Ifrac取不同的值时，优化得到的拓扑结 构基本保持不变，说明该拓扑结构存在一理想状态，通过不断尝试，得到如下发现: 在本例中当vol frac=O. 35时所得的拓扑结构最为规律，并且该拓扑结构趋向 与桁架结构。(具体拓扑结构演变过程见附件视频动画)第6章结论及总结本文中，从两条独立的途径来分别对该问题进行研究。Hyper works中，利用hype r mesh自带的的模块对该问题进行了建模，网 格划分，得到有限元模型，然后将有限元模型导入拓扑优化模块Opt i st r u c t, 添加边界条件及约束，定义目标函数和约束条件，对C形夹进

24、行拓扑结构优化。 所得到的拓扑结构趋向于桁架，说明桁架结构有很高的刚度。Matlab中，首先对该问题构建数学模型，采用SI MP(OC)算法对数学模 型进行优化求解，由于本文所研究的问题的求解域并不是规则的矩形，而是由三 部分矩形叠加，所以在使用SIMP法时，对总体刚度矩阵组装，添加边界和约束 条件成为了本文的关键技术，涉及到较多的有限元理论。在mat 1 a b中将每次迭 代所得到的处理结果进行了图像动态显示，以此来清晰的观察拓扑结构优化的动 态过程，给人以直观的印象。在用Optistru c t进行拓扑结构优化的时候，我们发现，当材料去除率为 60%时，其所得到的拓扑结构与我们用Matla

25、b进行拓扑结构优化是所得到的结 果的拓扑结构是一致的。这验证了我们的matlab数学模型是对的。参考文献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罗震.基于变密度法的连续体结构拓扑优化设计技术研究.博士学位论文.20 0 53 孙靖民.机械优化设计第三版.机械工业出版社.200 4白新理等.结构优化设计.黄河水利出版社.2 0 0 8.4郭仁生.机械工程设计分析和MATLAB应用.机械工业出版

26、社.200 6附件：C形夹拓扑优化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) ;%调用子函数得到整体位移

27、向量%将求解域分成三块矩形进行目标函数定义和敏度分析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区)、口1*110'*1<£*110；%目标函

28、数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;% 目标函数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 ; % 目标函数 dc(elyxlx) = -penal * x (e 1 y,el x ) A(pe n al-l)*U e "*KE*U e;end end%敏度过滤，避免拓扑结构边界锯齿传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 ) %指定非优化区域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%最优性判据更新单元密度值%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