华中科技大学研究生课程考试答题本 基于MATLAB和Optistruct的C形夹拓扑优化学 院（部）： 机械科学与工程学院 课 程 名 称： 工程优化设计 学 生 姓 名： 班 级： 机硕1107 班 学 号： 2012年01月10日目录第1章 选题背景介绍及问题描述 31.1选题背景及意义 31.1.1工程背景及基本原理 31.1.2 本文研究意义 31.2研究现状 31.2.1 理论研究现状 31.2.2 应用研究现状 41.3 该研的意义 4第2章 SIMP变密度法理论基础 52.1 SIMP密度刚度插值法理论基础 52.2 拓扑优化的数学模型 52.3 优化准则的基本理论 6第3章 优化设计的数学模型及解决方案 73.1问题描述 73.2 优化问题的数学模型 73.3 模型分析求解方法选择 8第四章 拓扑优化步骤及结果 94.1 基于matlab的变密度拓扑优化 94.1.1 问题求解的关键技术及代码 94.1.2 拓扑优化结果及分析 104.2 基于Optistruct的C形夹拓扑优化 114.2.1有限元模型的建立 114.2.2 基于Optistruct的拓扑优化结果及分析 12第6章 结论及总结 14参考文献 15附件：悬臂梁拓扑优化matlab程序 15第1章 选题背景介绍及问题描述1.1选题背景及意义1.1.1工程背景及基本原理通常把结构优化按设计变量的类型划分成三个层次：结构尺寸优化、形状优化和拓扑优化。
尺寸优化和形状优化已得到充分的发展，但它们存在着不能变更拓扑结构的缺陷在这样的背景下，人们开始研究拓扑优化拓扑优化的基本思想是将寻求结构的最优拓扑问题转化为在给定的设计区域内寻求最优材料的分布问题寻求一个最佳的拓扑结构形式有两种基本的原理：一种是退化原理，另一种是进化原理退化原理的基本思想是在优化前将结构所有可能杆单元或所有材料都加上，然后构造适当的优化模型，通过一定的优化方法逐步删减那些不必要的结构元素，直至最终得到一个最优化的拓扑结构形式进化原理的基本思想是把适者生存的生物进化论思想引入结构拓扑优化，它通过模拟适者生存、物竞天择、优胜劣汰等自然机理来获得最优的拓扑结构 1.1.2 本文研究意义目前，结构优化大部分集中在尺寸设计变量 (如板厚、杆的剖面积及管梁的直径 )拓扑结构优化较尺寸优化复杂 ,但对于有些问题拓扑结构优化比尺寸优化有效 ,C形夹是其中的例子之一本文讨论C形夹的拓扑优化问题 ,围绕这一问题,怎样使结构具有最大刚度的设计占有相当重要的地位;怎样优化结构的形状使材料的分布，更加合理从而达到使结构具有最小柔度的目的是本文要研究的问题1.2研究现状1.2.1 理论研究现状结构拓扑优化是近20年来从结构优化研究中派生出来的新分支，它在计算结构力学中已经被认为是最富挑战性的一类研究工作。
目前有关结构拓扑优化的工程应用研究还很不成熟，在国外处在发展的初期，尤其在国内尚属于起步阶段1904年Michell在桁架理论中首次提出了拓扑优化的概念自1964年Dorn等人提出基结构法，将数值方法引入拓扑优化领域，拓扑优化研究开始活跃20世纪80年代初，程耿东和N.Olhoff在弹性板的最优厚度分布研究中首次将最优拓扑问题转化为尺寸优化问题，他们开创性的工作引起了众多学者的研究兴趣1988年Bendsoe和Kikuchi发表的基于均匀化理论的结构拓扑优化设计，开创了连续体结构拓扑优化设计研究的新局面1993年XieYM和StevenGP提出了渐进结构优化法1999年Bendsoe和Sigmund证实了变密度法物理意义的存在性2002年罗鹰等提出三角网格进化法，该方法在优化过程中实现了退化和进化的统一，提高了优化效率1.2.2 应用研究现状在前人提出的重要理论基础上，后人也将其跟其他现代设计的方法相结合，衍生出了其他一些拓扑结构优化方法：如与可靠性相结合的情况下，MAUTE等应用变密度法并结合可靠性分析对一微机电系统进行了基于可靠性的拓扑优化设计，PAPADRAKAKIS等将遗传算法应用于具有可靠性约束的桁架结构拓扑优化设计中，国内学者马洪波也对基于遗传算法的结构可靠性优化问题进行了讨论。
华南理工大学机械工程学院欧阳高飞等对基于水平集方法的结构可靠性拓扑优化进行了研究1.3 本文研究的意义通过这次的作业加深对工程优化算法的学习和使用，提高对拓扑优化的方法和过程的了解和学习另外对相关软件软件的应用能够达到一个新的高度这些不仅能使我们现在的知识体系得到充实和优化，而且也是我们今后人生的财富第2章 SIMP变密度法理论基础2.1 SIMP密度刚度插值法理论基础SIMP模型主要通过引入惩罚因子，在材料的弹性模量和单元相对密度之间建立起一种显示的非线性对应关系它的作用是当设计变量的值在(0,1)之间时，对中间密度值进行惩罚，使中间密度值逐渐向0/1两端聚集，这样可以使连续变量的拓扑优化模型能很好地逼近原来0-1离散变量的优化模型这时中间密度单元对应一个很小的弹性模量，对结构刚度矩阵的影响将变得很小，可以忽略不计SIMP材料模型的数学表达形式： （2.1） （2.2）其中：为两数学模型中对中间密度材料的惩罚因子惩罚因子的作用是当设计变量的值在(0,1)之间时，通过逐渐增加的值对设计变量的中间值进行惩罚，随着值的增大，设计逐渐接近0/1设计。
为有效压缩中间密度材料，要求表示插值以后的弹性模量，和分别为固体和空洞部分材料的弹性模量，，/1000表示单元的设计变量表示插值以后的刚度矩阵，表示第个单元固体材料的刚度矩阵2.2 拓扑优化的数学模型以结构的柔度最小化(或刚度最大化、应变能最小化)作为优化的目标函数，以结构整体的体积约束作为优化的约束条件刚度优化的数学模型表示为： （2.3）SIMP对应的柔度函数和敏度形式： （2.4） （2.5）其中：以上各式中，表示第个单元的刚度矩阵表示结构的位移向量；表示设计变量，为避免总刚度矩阵奇异，取=0.001为单元数目，表示结构的柔度，表示柔度关于设计变量的敏度2.3 优化准则的基本理论刚度拓扑优化问题，是典型的具有不等式约束的非线性规划问题不等式约束多元函数极值的必要条件是Kuhn-Tucker条件，它是采用优化准则法求解非线性优化问题的重要理论引入对设计变量上下限约束的拉格朗日乘子、以及对体积约束的拉格朗日乘子，构造Lagrange函数为如下形式：(2.6)对于Lagrange函数，当时，设计变量的上下限侧面约束均不起作用，设计变量是主动变量。
主动变量在迭代过程中作为设计变量允许发生改变，；当时，仅设计变量下限约束起作用，设计变量为被动变量，；当时，仅设计变量上限约束起作用，设计变量也为被动变量，被动变量在迭代过程中不能变化，只能由侧面约束的边界值来确定[2]第3章 优化设计的数学模型及解决方案3.1问题描述如图所示，C形夹在自由端口受到三角形分布力F的作用，要求在保持对原材料体积一定缩减比的情况下对原实体悬臂梁做结构拓扑优化设计，优化目标是使结构刚度最大优化的结果应该使原设计区域产生孔洞，使结构拓扑发生变化原实体C形夹为如下图3.1所示的C形，尺寸如下图所示，材料为45钢，密度ρ为，弹性模量，泊松比图3.1 C形夹的尺寸和受力（单位：mm）3.2 优化问题的数学模型该问题中，要求同时满足刚度最大，质量最轻，这两个变量若同时改变，则问题复杂度太大，并且可能导致问题不可求解所以我们采用在确定的质量下，来讨论刚度最大的问题由于对特定的材料，其质量和体积有一定的关系，并且我们采用去除法的思想来建立模型的，故我们可以采用给优化后的体积与优化前的体积比赋确定的值，来达到在给定质量条件下满足刚度最大的问题其数学模型如下： (1)注：其中为结构变形能，U为结构变形总位移矩阵，K为结构总刚度矩阵，N为划分单元总数， 为单元位移向量，为单元刚度，（由于划分单元的时候，我们采用等分矩形单元，所以每个单元的刚度可用一个常量来处理）是拓扑结构优化过程中变化着的体积，为未经过优化前悬臂梁的体积。
F为结构所受的三角形载荷为悬臂梁的相对密度3.3 模型分析求解方法选择对该问题是用两种方法求解方法一：基于matlab的变密度拓扑优化法；方法二：是用成熟的有限元拓扑优化软件Optistruct进行优化基于matlab的变密度拓扑优化该问题的优化方法有很多种，常用的有如下方法：Optimality Criteria(OC) methods,（优化准则方法）Sequential Linear Programming (SLP) methods（序列线性规划法）Method of Moving Asymptotes (MMA bySvanberg 1987)等为了简化问题的复杂度，此处我们采用standard OC-method.方法来实现在处理过程中，关于设计变量相对密度x每一步的更新，我们采用在1995年提出的如下算法来实现第四章 拓扑优化步骤及结果4.1 基于matlab的变密度拓扑优化4.1.1 问题求解的关键技术及代码一般而言，由于OC法所使用的单元是矩形，所以OC法很适合求解求解域为矩形的优化问题，而本文选题为一C形结构，若使用OC法，则需将C形结构划分成为三个矩形的集合，因而在整体刚度矩阵的组装方面应该考虑如何进行，由于本程序是参考经典99行OC法拓扑优化程序改变，对于程序中避免边界锯齿现象所使用的check函数，该如何进行修改，以及如何进行边界条件的施加。
优化问题的初值条件：Nelx=60 x方向单元的数目为60nely=50 y方向单元的数目为50volfrac=0.35 保留原材料的体积分数为0.35penal=3.0 抑制权值为3.0（该取值是资料建议的典型值）rmin=1.2 过滤大小为1.2（该取值是资料建议的典型值）关键代码：function [xnew]=OC(nelx,nely,x,volfrac,dc) %%定义OC优化准则函数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)))));%%根据和敏度dc更新x的值%%% if sum(sum(xnew)) - volfrac*nelx*nely > 0;%%根据迭代过程中体积比是否达到预设的体积比，判断迭代是否继续进行%% l1 = lmid; else l2 = lmid; endend for ely = 1:nely for elx = 1:(nelx-20) %%取C形夹的左半矩形区域 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;%目标函数 dc(ely,elx) = -penal*x(ely,elx)^(penal-1)*Ue*KE*Ue;%%敏度值dc end end4.1.2 拓扑优化结果及分析 实验结果分析：从该实验结果来看，在我们给定的体积保留率的情况下，每经过一次拓扑结构优化，该优化程序就将C形夹的拓扑结构中强度要求不高处材料的密度减小，直到所有无用的材料都将被去除为止。
我们的拓扑结构优化模型是建立在结构变形能最小、体积去除率自己给定的基础上进行的，故我们可以根据实际情况，自行确定体积去除率在拓扑优化的过程中，我们可以观察到，我无论体积压缩率如何变化，C形夹模型最终都向桁架结构进化这说明，在结构件中，在自身材料多少相同的条件下，桁架具有很高的刚度和强度，其实这也就是为什么拓扑结构优化首先在桁架结构领域提出故工程上，我们常见工程人员采用桁架结构来作为一些工程的支撑结构，如塔吊等在验证不同的体积压缩率时，在不同的给定体积压缩率下，算法的有效性也不同，但在验证过程中，发现算法一直会收敛体积压缩率小的时候，该算法能很快终止；体积压缩率较大的时候，该算法的收敛速度较慢，并且还会出现不同程度的震荡，并且体积压缩率越大，该算法的振动也震荡4.2 基于Optistruct的C形夹拓扑优化4.2.1有限元模型的建立在Optistruct中建立有限元模型并划分网格相关参数为：;添加边界调节和载荷，得到有限元模型如下图5.1所示：图5.1 有限元模型4.2.2 基于Optistruct的拓扑优化结果及分析当体积分数volfrac=0.35时，拓扑结果图如图5.2所示：图5.2 C形夹拓扑优化结果图（volfrac=0.35）当体积分数volfrac=0.5时，拓扑结果图如图5.3所示：图5.3 C形夹拓扑优化结果图（volfrac=0.5）拓扑优化结构分析：在以变形最小为目标函数，即结构刚度最大，一定的材料保留比率volfrac为约束调节，当volfrac取不同的值时，优化得到的拓扑结构基本保持不变，说明该拓扑结构存在一理想状态，通过不断尝试，得到如下发现：在本例中当volfrac=0.35时所得的拓扑结构最为规律，并且该拓扑结构趋向与桁架结构。
具体拓扑结构演变过程见附件视频动画）第6章 结论及总结本文中，从两条独立的途径来分别对该问题进行研究Hyperworks中，利用hypermesh自带的的模块对该问题进行了建模，网格划分，得到有限元模型，然后将有限元模型导入拓扑优化模块Optistruct，添加边界条件及约束，定义目标函数和约束条件，对C形夹进行拓扑结构优化所得到的拓扑结构趋向于桁架，说明桁架结构有很高的刚度Matlab中，首先对该问题构建数学模型，采用SIMP(OC)算法对数学模型进行优化求解，由于本文所研究的问题的求解域并不是规则的矩形，而是由三部分矩形叠加，所以在使用SIMP法时，对总体刚度矩阵组装，添加边界和约束条件成为了本文的关键技术，涉及到较多的有限元理论在matlab中将每次迭代所得到的处理结果进行了图像动态显示，以此来清晰的观察拓扑结构优化的动态过程，给人以直观的印象在用Optistruct进行拓扑结构优化的时候，我们发现，当材料去除率为60%时，其所得到的拓扑结构与我们用Matlab进行拓扑结构优化是所得到的结果的拓扑结构是一致的这验证了我们的matlab数学模型是对的参考文献[1]O.sigmund. A 99 line topology optimization code written in Matlab.Educational article, 2001[2]罗震.基于变密度法的连续体结构拓扑优化设计技术研究.博士学位论文.2005[3]孙靖民.机械优化设计第三版.机械工业出版社.2004[4]白新理等.结构优化设计.黄河水利出版社.2008.4[5]郭仁生.机械工程设计分析和MATLAB应用.机械工业出版社.2006附件：C形夹拓扑优化matlab程序function top(nelx,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); %调用子函数得到整体位移向量 % 将求解域分成三块矩形进行目标函数定义和敏度分析 [KE] = lk; c = 0.; for ely = 1:nely for elx = 1:(nelx-20) 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;%目标函数 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)+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;%目标函数 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 +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;%目标函数 dc(ely,elx) = -penal*x(ely,elx)^(penal-1)*Ue*KE*Ue; end end% 敏度过滤，避免拓扑结构边界锯齿化 [dc] = check(nelx,nely,rmin,x,dc); % 通过OC法更新单元密度 [x] = OC(nelx,nely,x,volfrac,dc); for j =(nelx-19):(nelx-1) %指定非优化区域 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 )])%显示结果 colormap(gray); imagesc(-x); axis equal; axis tight; axis off;pause(1e-6);%绘制密度图end %%%%%%%%%% 最优性判据更新单元密度值 %%%%%%%%%%%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%%%%%%%%%% 单元密度过滤 %%%%%%%%%%%%%%%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,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(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(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 = (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(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(ely,elx)^penal*KE; endendfor elx = (nelx-19):nelx for ely = 36:50 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(ely,elx)^penal*KE; endend% 定义载荷和边界条件i=0;for elx = (nelx-19):nelx ely=15; n2 = (nely+1)* elx +ely;F(2*n2,1) = i;i=i+1;end i=0;for elx = (nelx-19):nelx ely=36; n1 = (nely+1)*(elx-1)+ely; F(2*n1,1) = -i;i=i+1;endfixeddofs = [1,2,2*(nely+1)-1,2*(nely+1)];alldofs = [1:2*(nely+1)*(nelx+1)];freedofs = setdiff(alldofs,fixeddofs);U(freedofs,:) = K(freedofs,freedofs) \ F(freedofs,:); % 求解G U(fixeddofs,:)= 0;%%%%%%%%%% 单元刚度矩阵 %%%%%%%%%%%%%%%%%%%%%%%function [KE]=lkE = 200000; nu = 0.3;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-nu^2)*[ k(1) k(2) k(3) 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)];。