基于MATLAB和ANSYS的悬臂梁拓扑优化

1、现代设计理论和方法大作业-基于MATLA和ANSYS勺悬臂梁拓扑优化指导老师：志冈H 项目组长：薛亚波 项目成员：机自 66 学院：机械工程学院基于MATLA和ANSYS勺悬臂梁拓扑优化、计划和任务安排表1.1学习时间安排表:任务时间查阅整理资料5月03日一5月08日学习相关知识和论文5月09日一5月15日设计实验方案5月16日一5月17日学习优化理论5月18日一5月20日撰写论文5月21日一5月22日集体交流学习心得5月23日一5月24日1.2任务分配表学号任务06011136齐寰宇查阅资料、答辩PPT制作06011137王晓庆查阅资料、答辩PPT制作06011139冰查阅资料、研究问题分析

2、06011141康辰龙查阅资料、ansys拓扑模 块研究06011143梁伟云查阅资料、项目背景研究06011145卢岗查阅资料ansys拓扑模块 研究06011153清伟查阅资料ansys拓扑模块 研究06011154薛亚波查阅资料、后期工作总结06011155严宏伟查阅资料ansys拓扑模块 研究06011156志亮查阅资料、项目背景研究06011157龙查阅资料、答辩PPT制作06012020乔百杰查阅资料、Matlab结构优 化06073012邓磊查阅资料、项目背景研究06182047立峰查阅资料、研究问题分析06055164杜博查阅资料、研究问题分析06011146牛玉辉查阅资料、研

3、究问题分析06011147朴群星查阅资料、项目背景研究、项目背景介绍及问题描述2.1项目背景及意义：2.1.1工程背景及基本原理：通常把结构优化按设计变量的类型划分成三 个层次：结构尺寸优化、形状优化和拓扑优化。尺寸优化和形状优化已得到充分 的发展，但它们存在着不能变更拓扑结构的缺陷。在这样的背景下，人们开始研究拓扑优化。拓扑优化的基本思想是将寻求结构的最优拓扑问题转化为在给定的 设计区域寻求最优材料的分布问题。寻求一个最佳的拓扑结构形式有两种基本的 原理：一种是退化原理，另一种是进化原理。退化原理的基本思想是在优化前将 结构所有可能杆单元或所有材料都加上，然后构造适当的优化模型，通过一定的优

4、化方法逐步删减那些不必要的结构元素，直至最终得到一个最优化的拓扑结构 形式。进化原理的基本思想是把适者生存的生物进化论思想引入结构拓扑优化， 它通过模拟适者生存、物竞天择、优胜劣汰等自然机理来获得最优的拓扑结构。2.1.2本文研究意义：目前，结构优化大部分集中在尺寸设计变量 （如板 厚、杆的剖面积及管梁的直径）。拓扑结构优化较尺寸优化复杂，但对于有些问 题拓扑结构优化比尺寸优化有效，悬臂梁是其中的例子之一。本文讨论悬臂梁的 拓扑优化问题，围绕这一问题，怎样使结构具有最大刚度的设计占有相当重要的 地位；怎样优化结构的形状使材料的分布，更加合理从而达到使结构具有最大刚 度的目的是本文要研究的问题。

5、2.2研究现状2.2.1理论研究现状：结构拓扑优化是近20年来从结构优化研究中派生出来 的新分支，它在计算结构力学中已经被认为是最富挑战性的一类研究工作。目前有关结构拓扑优化的工程应用研究还很不成熟，在国外处在发展的初期，尤其在 国尚属于起步阶段。1904年Michell在桁架理论中首次提出了拓扑优化的概念。 自1964年 Dorn等人提出基结构法，将数值方法引入拓扑优化领域，拓扑优化研究 开始活跃。20世纪80年代初，程耿东和N.OIhoff在弹性板的最优厚度分布研究中 首次将最优拓扑问题转化为尺寸优化问题，他们开创性的工作引起了众多学者的 研究兴趣。1988年Bendsoe和Kikuchi

6、发表的基于均匀化理论的结构拓扑优化设 计，开创了连续体结构拓扑优化设计研究的新局面。1993年XieYM和Steve nGP提出了渐进结构优化法。1999年 Ben dsoe和Sigm unc证实了变密度法物理意义的存在 性。2002年罗鹰等提出三角网格进化法，该方法在优化过程中实现了退化和进化 的统一，提高了优化效率。2.2.2应用研究现状：在前人提出的重要理论基础上，后人也将其跟其他现代设 计的方法相结合，衍生出了其他一些拓扑结构优化方法：如与可靠性相结合的情况下，MAUT等应用变密度法并结合可靠性分析对一微机电系统进行了基于可靠 性的拓扑优化设计，PAPADRAKAKIS将遗传算法应用于

7、具有可靠性约束的桁架结 构拓扑优化设计中，国学者马洪波也对基于遗传算法的结构可靠性优化问题进行 了讨论。华南理工大学机械工程学院欧阳高飞等对基于水平集方法的结构可靠性 拓扑优化进行了研究。2.3研究目标：2.3.1设计目标：设计一悬臂梁在自由端受到竖直向下集中力P的作用，要求保持悬臂梁长度尺寸不变。2.3.2优化目标：在满足以上条件情况下对原实体悬臂梁做结构拓扑优化设 计，使其结构刚度最大、重量最轻。（优化的结果应该使原设计区域产生孔洞， 使结构拓扑发生变化。）优化目标分析：要使梁满足以上条件情况必须同时兼顾重量和 刚度问题。重量太轻不能满足刚度问题。刚度太大又必须足够的材 料以满足其要求。将

8、两个因素进行耦合分析，我们可以得出的结论 是必须通过结构层次来进行优化。连续体结构拓扑优化较成熟的方 法主要有均匀化方法、变密度方法和渐进结构优化法，通过综合比 较我们最终选用变密度法。2.3.3实现方法：本文在对题目进行分析后，决定从两个独立的方向来分别 研究，以此来互相验证，保证结果的可靠性。实验中要使用ansys和matlab两种软件分别独立进行优化分析，然后将两种结果来进行比较。2.3.4实践目标：通过这次的作业我们期望通过搜索资料，团体讨论，分组作业的方式，以期望达到对拓扑优化的方法和过程有个具体的了解和学习。另外对分析软件的应用能够达到一个新的高度。这些不仅能使我们现在的知识体系得

9、 到充实和优化，而且也是我们今后人生的财富2.4研究容2.4.1问题描述如图所示，悬臂梁在自由端受到竖直向下集中力 P的作用，要求在保持悬臂 梁长度尺寸不变的情况下对原实体悬臂梁做结构拓扑优化设计，优化目标是使结构刚度最大、重量最轻。（优化的结果应该使原设计区域产生孔洞，使结构拓扑 发生变化。）原实体悬臂梁横截面为矩形，矩形宽度为 0.2m,悬臂梁的材料为45钢，密 度p为7.8 x 10-6kg/mm3,弹性模量E=2X 105MPa许用切应力t =60MPa竖 直向下集中力为P 1.5 104n。图 悬臂梁受集中力作用(单位： m242问题转化本研究中，要求在保持悬臂梁在结构尺寸不变的情况

10、下， 使悬臂梁的刚度最 大、质量最轻。一般悬臂梁都为三维立体结构，载荷是在竖直方向上加载的，所以在垂直纸 面方向上，并且在悬臂梁的对称平面上加载，根据材料力学知识可知，所加载在 和只对加载区周围较小围产生影响， 其余部分与加载形式无关。其受力与如图所 示平面的手里相同，故可以将问题抽象为平面问题来处理。当占有区域S的一个结构上作用有强度为P ( S )的载荷时,使结构具有最大 刚度的问题和使结构所受外力功 W具有最小值的问题是等价的。结构所受的外力 功W与结构变形能C之间的关系为W = 2C,所以我们可以将该问题处理为求在一 定的约束条件下，该悬臂梁的的最小变形能最小的问题。243数学模型：该

11、问题中，要求同时满足刚度最大，质量最轻，这两个变量若同时改变， 则问题复杂度太大，并亲切可能导致问题不可求解。所以我们采用在确定的质量 下，来讨论刚度最大的问题。由于对特定的材料，其质量和体积有一定的关系， 并且我们采用去除法的思想来建立模型的， 故我们可以采用给优化后的体积与优 化前的体积比赋确定的值，来达到在给定质量条件下满足刚度最大的问题。 其数 学模型如下：Nmin : c(x) = UtKU 二uj k0 ueC=1subject to: f:KU= F:0 Xm 迫 xl(1)注：其中C (x)为结构变形能，U为结构变形总位移矩阵，K为结构 总刚度矩阵，N为划分单元总数，Ue为单元

12、位移向量，ko为单元刚度， (由于划分单元的时候，我们采用等分举行单元，所以个单元的刚度可用 一个常量来处理)V(x)是拓扑结构优化过程中变化着的体积，Vo为未经过优化前悬臂梁的体积。F为结构所受的总载荷。x为悬臂梁的相对密度。243模型分析求解：该问题的优化方法有很多种，常用的有如下方法：Optimality Criteria(OC) methods, (优化准则方法)Seque ntial Linear Programmi ng (SLP) methods(序列线性规划法)Method of Movi ng Asymptotes (MMA bySva nberg 1987) 等 为了简化问

13、题的复杂度，此处我们采用 standard OC-method.方法来实 现。在处理过程中，关于设计变量相对密度x每一步的更新，我们采用在1995 年提出的如下算法来实现：碍榨=:max(厲金亡m):if %石2 一 m),I if m) xeB niin( 1, re + m),:+ m)i if mm(l5 xe + m T注：其中m是移动限制量，即相对密度变化的最大步长，更新条件中，加上最小边界条件判断，是防止相对密度变为0,使问题出 现奇异；加上相对密度最大边界条件判断，是为了防止出现相对密度大于1,是问题是去物理意义。是数值阻尼系数。Be可以由以下条件来确定:其中 为拉格朗日乘数，由

14、双向切片算法来确定目标函数对单元相对密度的变化率，即单元敏感度如下确定：deb丁二-PMP uf k0 % (4)为了确保该问题的拓扑结构优化解的存在，在处理技术上，我们得对该 问题的求解过程加些限制。此处我们采用过滤技术来进行处理。我们采用网格独立性滤波器来对单元敏感度的计算来进行改进，其改进 算法如下：de1Nf=i工丑宀却/=1丁其中Hf为卷积因子(即权重因子)，其由以下表达式来确定：/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设置材料属性: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 定义期望目标值的减少量，即期望减少的质量或体积分数:ScopeGeometryAll BodiesDefin iti onTargetReduct ion50. %32143 拓扑结果如下:32144 其他两种结果:定义质量减少率为 80%ScopeGeometryAll BodiesDefin iti onTargetReducti on80. %拓扑结果:60%勺

18、拓扑结果如下定义质量减少率为 60%ScopeGeometryAll BodiesDefin iti onTarget Reducti on60. %321.4 试验结果分析及结论：从拓扑结构我们可以看到，在不同的材料去除率下，总的趋势都 是边界面梁，其中包含着等强度设计的思想，这种拓扑结果较为符合 现实情况，在工程上具有很强的实用性。在不同的材料去除率下，我们发现结果具有不同的拓扑结构。3.2.2拓扑结构优化：注：以下结果是在材料去除率为40%勺情况下进行的优化。拓扑结构优化前，该结构的各项参数如下：Model Shape Optimization Solution ResultsObjec

19、t NameShape FinderStateSolvedScopeGeometryAll BodiesDefinitionTarget Reduction60. %ResultsOriginal Mass2512. kgOptimized Mass1188.9 kgMarginal Mass35.046 kg3.2.2.1 优化分析：在上面Ansys图中，我们可以发现，改图朴结构十分不规则 工程实际中，这样的结构在加工上很难实现。基于此，我们决定对 上面的结果进行优化，得到较为规则的拓扑结构，以便工程上的加 工。3.2.2.1有限元分析：优化前材料去除率为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总位移变形图应变图应力图优化后的分析结果,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 优化设计：根据拓扑结果，我们设置了 5个参数作为优化设计变量，其各参数 取值如下：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.我们采用多目标优化，将质量和结构最大变形作为目标函数，使其取值最小。在 将多目标转化为单目标时，我们对质量和最大变形分别赋予了不同的权重，具体设置如 下表：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优化结果：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通过对比以上三种结果，我们发现结果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.响应面：3.222 优化实验结果:3.221 重要参数设计：Nelx=32 x 方向单元的数目为 32nely=20 y 方向单元的数目为 20注：该悬臂梁的长宽比为16: 10f/volfrac=0.4 pen al=3.0保留原材料的体积分数为0.4抑制权值为3.0 （该取值是资料建议的典型值）=1.2 过滤大小为1.2 （该取值是资料建议的典型值）3.2.2.3试验结果分析从该实验结果来看，在我们给定的体积保留率的情况下， 每经过一次拓扑结构优化，该优化程序就将悬臂梁的拓扑结构 中强度要求不高处材料的密度减小，直到所有无用

26、的材料都将 被去除为止。我们的拓扑结构优化模型是建立在结构变形能最 小、体积去除率自己给定的基础上进行的，故我们可以根据实 际情况，自行确定体积去除率。在拓扑优化的过程中，我们可以观察到，我无论体积压缩 率如何变化，悬臂梁模型最终都向桁架结构进化。这说明，在 结构件中，在自身材料多少相同的条件下，桁架具有很高的刚 度和强度，其实这也就是为什么拓扑结构优化首先在桁架结构 领域提出。故工程上，我们常见工程人员采用桁架结构来作为 一些工程的支撑结构，如塔吊等。在实验中，在验证不同的体积压缩率时，我们发现，在不 同的给定体积压缩率下，算法的有效性也不同，但在验证过程 中，发现算法一直会收敛。体积压缩率

27、小的时候，该算法能很 快终止；体积压缩率较大的时候，该算法的收敛速度较慢，并 且还会出现不同程度的震荡，并且体积压缩率越大，该算法的 振动也震荡。3.3实验结论及总结：实验中，为了充分发挥我们小组人员的优势，我们从两条独立的途径来 分别对该问题进行研究。An sys 中，我们利用其自带的的模块对该问题进行了建模，拓扑结构优化。 为了工程上的加工方便，我们在拓扑结构优化的基础上，对其形状进行了规则化， 与此同时，进行再次优化。Matlab 中，我们将建立的数学模型用数值模拟的形式进行了处理，并将所 得到的处理结果进行了图像动态显示，以此来清晰的观察拓扑结构优化的动态过 程，给人以直观的印象。在用

28、Ansys进行拓扑结构优化的时候，我们发现，当材料去除率为60%寸,其所得到的拓扑结构与我们用 Matlab进行拓扑结构优化是所得到的结果的拓扑 结构是一致的。这验证了我们的数学模型是对的。但是我们的模型也存在着一定的问题，即在材料去除率很小的时候，预加 载荷未变的情况下，该模型的解均存在。出现这一情况的原因是，我们在计算悬 臂梁的位移矩阵，一直认为其应力未超过破坏强度，整个过程均在弹性变形围。 故在实际应用的时候，需要我们结合一定的工程经验来确定个重要参数的取值。四、参考文献及附录：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设计、仿真与优化兵何正嘉雪峰等清华大学2008年8月第一版4. .topopt.dtu.dk .5

30、. 基于结构最大刚度的形状优化方法欧阳高飞建民机械工程学报2008年10月4.2相关附件：附件1:悬臂梁拓扑结构优化Matlab程序：% A 99 LINE TOPOLOGY OPTIMIZATION CODE BY OLESIGMUND, OCTOBER 1999 %fun ctio n top( nelx, nely,volfrac,pe nal,rmi n);%主程序，即调用程序。其中的参数需要用户提供。依次为x方向单元数、y方向单元数、期望优化的体 积或质量百分比、抑制权值、过滤因子% INITIALIZE 初始化x(1: nely,1: nelx) = volfrac; %在整个求解

31、域均分材料loop = 0;cha nge = 1.; % 新旧改变之差% START ITERATIONwhile cha nge 0.01 %第一层，设置图像显示次数，即优化的次数loop = loop + 1;xold = x;% FE-ANALYSIS 有限元U=FE( nelx, nely,x,pe nal); %调用有限元子程序，返回位移矢量% OBJECTIVE FUNCTION AND SENSITIVITY ANALYSISKE = lk; %调用单元刚度阵c = 0.;for ely = 1:nely %16-24行一次循环。外层循环，y方向循环nely次数for elx

32、= 1:nelx %层循环，x方向循环次数为nelxn1 = (n ely+1)*(elx-1)+ely; %在全局坐标系中，左上方的节点数目n2 = (n ely+1)* elx +ely;%在全局坐标系中，右上方的节点数目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); %调用网格滤波子函数% DESIGN UPDATE BY THE OPTIMALITY CRITERIA METH 优化x = OC( nelx, nely,x,volfrac,dc); %调用优化子函数% PRINT RESULTS拓扑过程显示cha nge = max(max(abs(x-xold);disp( 循环次数 It.: sprintf(%4i,loop)优化目标值 Obj.:spri ntf(%10.4f ,c).拓扑质

34、量分数 Vol.: sprintf(%6.3f,sum(sum(x)/(nelx*nely).终止条件ch.: sprintf(%6.3f,change ) %优化过程指标显示% PLOT DENSITIES 密度分布colormap(gray); imagesc(-x); axis equal; axis tight; axis off;pause(1e-6); %图像显示设置 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 %基于标准优化方法的优化程序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