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

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

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

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

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

5、有相当重要的 地位; 怎样优化结构的形状使材料的分布，更加合理从而达到使结构具有最大刚 度的目的是本文要研究的问题。研究现状2.2.1理论研究现状： 结构拓扑优化是近 20年来从结构优化研究中派生出来 的新分支， 它在计算结构力学中已经被认为是最富挑战性的一类研究工作。 目前 有关结构拓扑优化的工程应用研究还很不成熟， 在国外处在发展的初期， 尤其在 国内尚属于起步阶段。 1904年Michell 在桁架理论中首次提出了拓扑优化的概念。 自 1964年 Dorn等人提出基结构法， 将数值方法引入拓扑优化领域， 拓扑优化研究 开始活跃。20世纪80年代初，程耿东和在弹性板的最优厚度分布研究中首次

6、将最 优拓扑问题转化为尺寸优化问题，他们开创性的工作引起了众多学者的研究兴 趣。1988年 Bendsoe和Kikuchi 发表的基于均匀化理论的结构拓扑优化设计， 开创 了连续体结构拓扑优化设计研究的新局面。 1993年XieYM和 StevenGP提出了渐进 结构优化法。 1999年Bendsoe和Sigmund证实了变密度法物理意义的存在性。 2002 年罗鹰等提出三角网格进化法，该方法在优化过程中实现了退化和进化的统一， 提高了优化效率。2.2.2应用研究现状： 在前人提出的重要理论基础上， 后人也将其跟其他现代设 计的方法相结合， 衍生出了其他一些拓扑结构优化方法： 如与可靠性相结合

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

8、结构拓扑发生变化。）优化目标分析：要使梁满足以上条件情况必须同时兼顾重量和 刚度问题。 重量太轻不能满足刚度问题。 刚度太大又必须足够的材 料以满足其要求。 将两个因素进行耦合分析， 我们可以得出的结论 是必须通过结构层次来进行优化。 连续体结构拓扑优化较成熟的方 法主要有均匀化方法、 变密度方法和渐进结构优化法， 通过综合比 较我们最终选用变密度法。2.3.3实现方法： 本文在对题目进行分析后，决定从两个独立的方向来分别 研究，以此来互相验证，保证结果的可靠性。实验中要使用 ansys 和 matlab 两 种软件分别独立进行优化分析，然后将两种结果来进行比较。2.3.4实践目标： 通过这次

9、的作业我们期望通过搜索资料，团体讨论，分组 作业的方式， 以期望达到对拓扑优化的方法和过程有个具体的了解和学习。 另外 对分析软件的应用能够达到一个新的高度。 这些不仅能使我们现在的知识体系得 到充实和优化，而且也是我们今后人生的财富研究内容2.4.1问题描述如图所示，悬臂梁在自由端受到竖直向下集中力 P的作用，要求在保持悬臂 梁长度尺寸不变的情况下对原实体悬臂梁做结构拓扑优化设计， 优化目标是使结 构刚度最大、重量最轻。 （优化的结果应该使原设计区域产生孔洞，使结构拓扑 发生变化。）原实体悬臂梁横截面为矩形，矩形宽度为 0.2m，悬臂梁的材料为 45 钢，密 度 为10-6kg/mm3，弹性

10、模量 E=2105MPa，许用切应力 =60MPa。竖直 向下集中力为 P 1.5 104N 。图 悬臂梁受集中力作用（单位： m）2.4.2问题转化本研究中，要求在保持悬臂梁在结构尺寸不变的情况下， 使悬臂梁的刚度最 大、质量最轻。一般悬臂梁都为三维立体结构， 载荷是在竖直方向上加载的， 所以在垂直纸 面方向上， 并且在悬臂梁的对称平面上加载， 根据材料力学知识可知， 所加载在 和只对加载区周围较小范围内产生影响， 其余部分与加载形式无关。 其受力与如 图所示平面的手里相同，故可以将问题抽象为平面问题来处理。当占有区域 S的一个结构上作用有强度为 P （ S ）的载荷时, 使结构具有最大 刚

11、度的问题和使结构所受外力功 W 具有最小值的问题是等价的。 结构所受的外力 功W 与结构变形能 C之间的关系为 W = 2C，所以我们可以将该问题处理为求在一 定的约束条件下，该悬臂梁的的最小变形能最小的问题。2.4.3数学模型： 该问题中，要求同时满足刚度最大，质量最轻，这两个变量若同时改变， 则问题复杂度太大， 并亲切可能导致问题不可求解。 所以我们采用在确定的质量 下，来讨论刚度最大的问题。由于对特定的材料，其质量和体积有一定的关系，并且我们采用去除法的思想来建立模型的， 故我们可以采用给优化后的体积与优 化前的体积比赋确定的值， 来达到在给定质量条件下满足刚度最大的问题。 其数 学模型

12、如下：注：其中 C（x）为结构变形能， U为结构变形总位移矩阵， K 为结构 总刚度矩阵， N 为划分单元总数， ue 为单元位移向量， ko 为单元刚度， （由于划分单元的时候，我们采用等分举行单元，所以个单元的刚度可用 一个常量来处理） V（x）是拓扑结构优化过程中变化着的体积， V0 为未经 过优化前悬臂梁的体积。 F为结构所受的总载荷。 x 为悬臂梁的相对密度。2.4.3模型分析求解： 该问题的优化方法有很多种，常用的有如下方法： Optimality Criteria（OC） methods,（优化准则方法）Sequential Linear Programming （SLP） me

13、thods（序列线性规划法）Method of Moving Asymptotes （MMA bySvanberg 1987） 等 为了简化问题的复杂度，此处我们采用 standard OC-method. 方法来实 现。在处理过程中，关于设计变量相对密度 x每一步的更新，我们采用在 1995 年提出的如下算法来实现：注：其中 m 是移动限制量，即相对密度变化的最大步长，更新 条件中， 加上最小边界条件判断， 是防止相对密度变为 0，使问题出 现奇异；加上相对密度最大边界条件判断，是为了防止出现相对密 度大于 1，是问题是去物理意义。 是数值阻尼系数。Be 可以由以下条件来确定：其中 为拉格朗

14、日乘数，由双向切片算法来确定。目标函数对单元相对密度的变化率，即单元敏感度如下确定：为了确保该问题的拓扑结构优化解的存在，在处理技术上，我们得对该 问题的求解过程加些限制。此处我们采用过滤技术来进行处理。我们采用网格独立性滤波器来对单元敏感度的计算来进行改进，其改进 算法如下：其中 Hf 为卷积因子（即权重因子），其由以下表达式来确定： 采用该算法对单元敏感度的更新进行改进， 是为了在更新单元 e的敏感度 时，将距其距离不超过 rmin 的临域内的单元来进行均匀化，防止各乡邻单元 的敏感度变化过大。注：其中 dist（e,f） 定义为单元体 e的中心到单元体 f 中心之间的 距离， rmin为

15、过滤大小。在优化过程中， 我们采用改进的敏感度表达 式（ 5）来代替表达式（ 4）来对（ 3）是进行计算。三、实验设计及结果分析：技术路线：在本实验中，由于我们提出了两种独立的方法，故我们将小组分为两部分 ;小组一：该组成员采用 Ansys 软件来对该悬臂梁模型进行建模分析，并进行拓 扑结构优化。小组二： 该组成员采用上面的数学模型，用 Matlab 原来来实现整个数学模型 的算法。并对处理结果进行动态显示，以便让我们清晰地看到整个拓扑结 构优化的过程。试验结果及分析：3.2.1 Ansys 分析的拓扑结构优化结果：3.2.1.2 建立几何模型：3.2.1. 几何参数设置如下：3.2.1. 在

16、 Ansys/workbench 中建立如下几何模型：3.2.1.2 网格化： 将几何模型直接调入有限元环境。 在 workbench 里建立的模型可 以无缝导入到有限元分析环境里。采用自动网格划分工具，得到的 模型如下：其中，网格单元和节点数设置如下：StatisticsNodes 节点数目3120Elements 单元数 目5613.2.1.3 施加载荷和设置边界条件选择优化类型为 “ shape optimization ”形状优化在上截面施加 15000N 的力3.2.1.3.2.1.4 求解：设置材料属性：StructuralYoungs Modulus 杨氏模量+011 PaPoi

17、ssons RatioDensity 密度7850. kg/m3Thermal Expansion1/ CTensile Yield Strength+008 PaCompressive Yield Strength+008 PaTensile Ultimate Strength+008 PaCompressive Ultimate Strength0. PaThermalThermal ConductivityW/mCSpecific Heat434. J/kg CElectromagneticsRelative Permeability10000ResistivityOhmm3.2.1.

18、定义期望目标值的减少量，即期望减少的质量或体积分数：ScopeGeometryAll BodiesDefinitionTargetReduction50. %3.2.1. 拓扑结果如下：3.2.1. 其他两种结果：定义质量减少率为 80%ScopeGeometryAll BodiesDefinitionTargetReduction80. %拓扑结果：定义质量减少率为 60%Reduction60%的拓扑结果如下3.2.1.4 试验结果分析及结论：从拓扑结构我们可以看到，在不同的材料去除率下，总的趋势都 是边界面梁，其中包含着等强度设计的思想，这种拓扑结果较为符合 现实情况，在工程上具有很强的

19、实用性。在不同的材料去除率下，我们发现结果具有不同的拓扑结构。3.2.2 拓扑结构优化： 注：以下结果是在材料去除率为 40%的情况下进行的优化。拓扑结构优化前，该结构的各项参数如下：3.2.2.1 优化分析：在上面 Ansys 图中，我们可以发现，改图朴结构十分不规则。 工程实际中，这样的结构在加工上很难实现。基于此，我们决定对 上面的结果进行优化，得到较为规则的拓扑结构，以便工程上的加 工。3.2.2.1 有限元分析：优化前材料去除率为 34%BoxLength XmLength Y1. mLength ZmPropertiesVolumem3MasskgOriginal Mass2512

20、 kgCAD Parametersds_d520ds_s150ds_D1320ds ss250ds dd990总位移变形图 应变图应力图优化后的分析结果，Object Name Total Deformation Equivalent Elastic Strain Equivalent StressMinimum0. mm/m10911 PaMaximummm/m+005 Pa3.2.2.1 优化设计：根据拓扑结果， 我们设置了 5 个参数作为优化设计变量， 其各参数 取值如下：我们采用多目标优化，将质量和结构最大变形作为目标函数，使其取值最小。在 将多目标转化为单目标时，我们对质量和最大变形

21、分别赋予了不同的权重，具体设置如 下表：优化结果：通过对比以上三种结果，我们发现结果 A 较为理想。Rating definitions:The most positive rating of means that the design is excellent in terms of satisfying the goal defined for the parameter.The most negative rating of means that the design is poor in termsof satisfying the goal.A rating of is neutra

22、l with respect to the goal.响应面 :3.2.3Matlab 语言实现的拓扑结构优化结果：3.2.2.1 重要参数设计：Nelx=32 x 方向单元的数目为 32 nely=20 y 方向单元的数目为 20 注：该悬臂梁的长宽比为 16：10f/volfrac= penal=保留原材料的体积分数为抑制权值为（该取值是资料建议的典型值）rmin=过滤大小为（该取值是资料建议的典型值）3.2.2.2优化实验结果：3.2.2.3 试验结果分析 从该实验结果来看，在我们给定的体积保留率的情况下， 每经过一次拓扑结构优化，该优化程序就将悬臂梁的拓扑结构 中强度要求不高处材料的密

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

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

25、人以直观的印象。在用 Ansys 进行拓扑结构优化的时候，我们发现，当材料去除率为 60%时， 其所得到的拓扑结构与我们用 Matlab 进行拓扑结构优化是所得到的结果的拓扑 结构是一致的。这验证了我们的数学模型是对的。但是我们的模型也存在着一定的问题，即在材料去除率很小的时候，预加 载荷未变的情况下， 该模型的解均存在。 出现这一情况的原因是， 我们在计算悬 臂梁的位移矩阵， 一直认为其应力未超过破坏强度， 整个过程均在弹性变形范围 内。故在实际应用的时候， 需要我们结合一定的工程经验来确定个重要参数的取四、参考文献及附录：主要参考文献：1.Sigmund, O. 1994: Design

26、of material structures using topology optimization. . Thesis, Department of Solid Mechanics, Technical University of Denmark2. Sigmund, O. 1997: On the design of compliant mechanismsusing topology optimization. Mech. Struct. Mach. 25, 4955263.ANSYS/WORKBEN设C计H、仿真与优化 李兵 何正嘉 陈雪峰等 清华大学出版社 2008 年 8 月第一版

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

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

29、= 1:nelx %内层循环， x 方向循环次数为 nelxn1 = (nely+1)*(elx-1)+ely; %在全局坐标系中，左上方的节点数目n2 = (nely+1)* elx +ely;%在全局坐标系中，右上方的节点数目 %n1 、n2 作用是从全局位移矢量提取局部或单元位移矢量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);% 返回单元位移c = c + x(ely,elx)penal*Ue*KE*Ue;dc(ely,elx) = -penal*x(ely,elx)(penal-1)*Ue*KE*Ue

30、;endend% FILTERING OF SENSITIVITIES灵敏度检验dc = check(nelx,nely,rmin,x,dc); % 调用网格滤波子函数% DESIGN UPDATE BY THE OPTIMALITY CRITERIA METH标O准D 优化x = OC(nelx,nely,x,volfrac,dc); %调用优化子函数% PRINT RESULTS拓 扑过程显示change = max(max(abs(x-xold);disp( 循环次数 It.: sprintf(%4i,loop) 优化目标值 Obj.: sprintf(% ,c) . 拓扑质量分数 Vo

31、l.: sprintf(%,sum(sum(x)/(nelx*nely) . 终止条件 ch.: sprintf(%,change ) % 优化过程指标显示% PLOT DENSITIES 密度分布colormap(gray); imagesc(-x); axis equal; axis tight; axis off; pause(1e-6); % 图像显示设置end% MESH-INDEPENDENCY FILTER % function dcn=check(nelx,nely,rmin,x,dc)dcn=zeros(nely,nelx);for i = 1:nelxfor j = 1:n

32、ely sum=;for k = max(i-round(rmin),1):min(i+round(rmin),nelx)for l = max(j-round(rmin),1):min(j+round(rmin), nely)fac = rmin-sqrt(i-k)2+(j-l)2);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);endend% FE-ANALYSIS % function U=FE(nelx,nely,x,

33、penal)KE=lk;K = sparse(2*(nelx+1)*(nely+1), 2*(nelx+1)*(nely+1);F = sparse(2*(nely+1)*(nelx+1),1); U =sparse(2*(nely+1)*(nelx+1),1);for ely = 1:nely % 全局刚度矩阵的形成for elx = 1:nelxn1 = (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,edo

34、f) = K(edof,edof) +x(ely,elx)penal*KE;endend% DEFINE LOADSAND SUPPORTS(HALF MBB-载BE荷A的M施) 加和边界条件的定 义%F(2,1) = -1; % 在单元力在左上方%fixeddofs = union(1:2:2*(nely+1),2*(nelx+1)*(nely+1);F(2*(nelx)*(nely+1)+2,1) = -1; fixeddofs = 1:2*(nely+1);alldofs = 1:2*(nely+1)*(nelx+1);freedofs = setdiff(alldofs,fixeddo

35、fs);% SOLVINU(freedofs,:) = K(freedofs,freedofs) F(freedofs,:); %添加约束条件，限制自由度U(fixeddofs,:)= 0; %freedofs 没有约束的自由度， fixeddofs 受到约束 的自由度% ELEMENT STIFFNESS MATRIX % function KE=lkE = 2.;nu = ;k=1/2-nu/6 1/8+nu/8 -1/4-nu/12 -1/8+3*nu/8 .-1/4+nu/12 -1/8-nu/8 nu/6 1/8-3*nu/8;KE = E/(1-nu2)*k(1) k(2) k(3

36、) k(4) k(5) k(6) k(7) k(8)k(2) k(1) k(8) k(7) k(6) k(5) k(4) k(3)k(3) k(8) k(1) k(6) k(7) k(4) k(5) k(2)k(4) 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) k(4)k(6) k(5) k(4) k(3) k(2) k(1) k(8) k(7)k(7) k(4) k(5) k(2) k(3) k(8) k(1) k(6)k(8) k(3) k(2) k(5) k(4) k(7) k(6) k(1);

37、% OPTIMALITY CRITERIA UPDATE %基%于%标%准%优%化%方 法的优化程序function xnew=OC(nelx,nely,x,volfrac,dc) % 返回最新的变量值，不 断更新变量l1 = 0; l2 = 100000; move = ;while (l2-l1 1e-4)lmid = *(l2+l1);xnew = max,max(x-move,min(1.,min(x+move,x.*sqrt(-dc./lmid);if sum(sum(xnew) - volfrac*nelx*nely 0; %bi-sectioning algorithml1 =

38、lmid; else l2 = lmid; end end% A 99 LINE TOPOLOGY OPTIMIZATION CODE BY OLESIGMUND, OCTOBER 1999 %function top(nelx,nely,volfrac,penal,rmin);% INITIALIZE x(1:nely,1:nelx) = volfrac;loop = 0;change = 1.;% START ITERATION while change loop = loop + 1;xold = x;% FE-ANALYSIS 有限元 U=FE(nelx,nely,x,penal);%

39、 OBJECTIVE FUNCTION AND SENSITIVITY ANALYSISKE = lk; c = 0.; for ely = 1:nely for elx = 1:nelx 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);c = c + x(ely,elx)penal*Ue*KE*Ue;dc(ely,elx) = -penal*x(ely,elx)(penal-1)*Ue*KE*Ue;endend

40、% FILTERING OF SENSITIVITIES灵敏度检验dc = check(nelx,nely,rmin,x,dc);% DESIGN UPDATE BY THE OPTIMALITY CRITERIA METH标O准D 优化x = OC(nelx,nely,x,volfrac,dc);% PR_NT RESULT也芋K就越引 change H max(max(abs*xo-d)八disp(r -*- spinff(-%4-00p) - obj:_sprinm%-o) :-stw幷s 午w幷_mid :(_2 土kxxnewH ma vo-fac*ne_x*ne_y V 0八二 H -mid八e-se_2 上mid-endend% MES HND EPENDENCYzlLTER % function Ldcnllcheck(ne_2+0-)2x su