下载文档原格式
拓扑优化（topology optimization）1. 基本概念拓扑优化是结构优化的一种。
结构优化可分为尺寸优化、形状优化、形貌优化和拓扑优化。
其中尺寸优化以结构设结构优化类型的差异计参数为优化对象，比如板厚、梁的截面宽、长和厚等；形状优化以结构件外形或者孔洞形状为优化对象，比如凸台过渡倒角的形状等；形貌优化是在已有薄板上寻找新的凸台分布，提高局部刚度；拓扑优化以材料分布为优化对象，通过拓扑优化，可以在均匀分布材料的设计空间中找到最佳的分布方案。
拓扑优化相对于尺寸优化和形状优化，具有更多的设计自由度，能够获得更大的设计空间，是结构优化最具发展前景的一个方面。
图示例子展示了尺寸优化、形状优化和拓扑优化在设计减重孔时的不同表现。
2. 基本原理拓扑优化的研究领域主要分为连续体拓扑优化和离散结构拓扑优化。
不论哪个领域，都要依赖于有限元方法。
连续体拓扑优化是把优化空间的材料离散成有限个单元（壳单元或者体单元），离散结构拓扑优化是在设计空间内建立一个由有限个梁单元组成的基结构，然后根据算法确定设计空间内单元的去留，保留下来的单元即构成最终的拓扑方案，从而实现拓扑优化。
3. 优化方法目前连续体拓扑优化方法主要有均匀化方法[1]、变密度法[2]、渐进结构优化法[3]（ESO）以及水平集方法[4]等。
离散结构拓扑优化主要是在基结构方法基础上采用不同的优化策略（算法）进行求解，比如程耿东的松弛方法[5]，基于遗传算法的拓扑优化[6]等。
4. 商用软件目前，连续体拓扑优化的研究已经较为成熟，其中变密度法已经被应用到商用优化软件中，其中最著名的是美国Altair公司Hyperworks系列软件中的Optistruc t和德国Fe-design公司的Tosca等。
前者能够采用Hypermesh作为前处理器，在各大行业内都得到较多的应用；后者最开始只集中于优化设计，而没有自己的有限元前处理器，操作较为麻烦，近年来和Ansa联盟，开发了基于Ansa的前处理器，但在国内应用的较少。
第二章拓扑优化什么是拓扑优化？拓扑优化是指形状优化，有时也称为外型优化。
拓扑优化的目标是寻找承受单载荷或多载荷的物体的最佳材料分配方案。
这种方案在拓扑优化中表现为“最大刚度”设计。
与传统的优化设计不同的是，拓扑优化不需要给出参数和优化变量的定义。
目标函数、状态变量和设计变量（参见“优化设计”一章）都是预定义好的。
用户只需要给出结构的参数（材料特性、模型、载荷等）和要省去的材料百分比。
拓扑优化的目标——目标函数——是在满足结构的约束（V）情况下减少结构的变形能。
减小结构的变形能相当于提高结构的刚度。
这个技术通过使用设计变量（ i）给每个有限元的单元赋予内部伪密度来实现。
这些伪密度用PLNSOL，TOPO命令来绘出。
例如，给定V=60表示在给定载荷并满足最大刚度准则要求的情况下省去60%的材料。
图2-1表示满足约束和载荷要求的拓扑优化结果。
图2-1a表示载荷和边界条件，图2-2b表示以密度云图形式绘制的拓扑结果。
图2-1 体积减少60%的拓扑优化示例如何做拓扑优化拓扑优化包括如下主要步骤：1．定义拓扑优化问题。
2．选择单元类型。
3．指定要优化和不优化的区域。
4．定义和控制载荷工况。
5．定义和控制优化过程。
6．查看结果。
拓扑优化的细节在下面给出。
关于批处理方式和图形菜单方式不同的做法也同样提及。
定义拓扑优化问题定义拓扑优化问题同定义其他线性，弹性结构问题做法一样。
用户需要定义材料特性（杨氏模量和泊松比），选择合适的单元类型生成有限元模型，施加载荷和边界条件做单载荷步或多载荷步分析。
参见“ANSYS Analysis Procedures Guides”第一、二章。
选择单元类型拓扑优化功能可以使用二维平面单元，三维块单元和壳单元。
要使用这个功能，模型中只能有下列单元类型：二维实体单元：SOLID2和SOLID82三维实体单元：SOLID92和SOLID95壳单元：SHELL93二维单元用于平面应力问题。
如何利用AN SYS进行拓扑优化就目前而言，利用有限元进行优化主要分成两个阶段：（1）进行拓扑优化，明确零件最佳的外形、刚度、体积，或者合理的固有频率，主要目的是确定优化的方向；（2）进行尺寸优化，主要目的是确定优化后的的零件具体尺寸值，通常是在完成拓扑优化之后，再执行尺寸优化。
在ANSYS中，利用拓扑优化，可以完成以下两个目的：（1）在特定载荷和约束的条件下，确定零件的最佳外形，或者最小的体积（或者质量）；（2）利用拓扑优化，使零件达到需要的固有频率，避免在使用过程中产生共振等不利影响。
1.ANSYS进行拓扑优化的进行拓扑优化的过程在ANSYS中，执行优化，通常分为以下6个步骤：1.1定义需要求解的结构问题对于结构进行优化分析，定义结构的物理特性必不可少，例如，需要定义结构的杨氏模量、泊松比（其值在0.1～0.4之间）、密度等相关的结构特性方面定义需要求解的结构问题，选择合理的优化单元类型，设定优化和非优化的区域定义载荷步或者需要提取的频率对优化过程进行定义和控制，计算并查看结果的信息，以供结构计算能够正常执行下去。
1.2选择合理的优化单元类型，在ANSYS中，不是所有的单元类型都可以执行优化的，必须满足如下的规定：（1）2D平面单元：PLANE82单元和P L ANE183单元；（2）3D实体单元：SOLID92单元和S O LID95单元；（3）壳单元：SHELL93单元。
上述单元的特性在帮助文件中有详细的说明，同时对于2D单元，应使用平面应力或者轴对称的单元选项。
1.3指定优化和非优化的区域在ANSYS中规定，单元类型编号为1的单元，才执行优化计算；否则，就不执行优化计算。
拓扑优化方法拓扑优化方法是一种有效的优化方法，目前被广泛应用于求解复杂优化问题。
本文通过介绍拓扑优化方法的基本原理、典型案例、优势与应用等方面，来深入探讨拓扑优化的相关知识。
一、什么是拓扑优化方法拓扑优化方法（Topology Optimization，简称TO）是一种解决复杂最优化问题的有效优化方法，它是利用拓扑的可变性，用于求解复杂拓扑结构组合优化问题的一种新兴方法。
拓扑优化方法既可以用来求解有限元分析（Finite Element Analysis，简称FEA）中有序结构问题，也可以用来求解无序结构问题。
二、拓扑优化方法的基本原理拓扑优化方法的基本原理是：在设定的最优化目标函数及运算范围内，利用优化技术，使得复杂结构拓扑结构达到最优，从而达到最优化设计目标。
拓扑优化方法的优势主要体现在重量最小化、强度最大化、结构疲劳极限优化等多种反向设计问题上。
此外，由于拓扑优化方法考虑到结构加工、安装、维护等方面，其结构设计更加实用性好。
三、拓扑优化方法的典型案例1、航空外壳优化：目前，航空外壳的拓扑优化设计可以使得外壳的重量减轻50%以上，同时提升外壳的强度和耐久性。
2、机械联轴器优化：拓扑优化方法可以有效的提高机械联轴器长期使用的耐久性，减少其体积和重量，满足高性能要求。
3、结构优化：通过拓扑优化方法，可以有效地减少刚性框架结构的重量，优化结构设计，改善结构性能，大大降低制造成本。
四、拓扑优化方法的优势1、灵活性强：拓扑优化方法允许在设计过程中改变结构形态，可以有效利用具有局部不稳定性的装配元件；2、更容易操作：拓扑优化方法比传统的有序结构模型更容易实现，不需要做过多的运算；3、成本低：拓扑优化方法可以有效降低产品的工艺制造成本，在改进出色性能的同时，可以节省大量人力物力；4、可重复性高：拓扑优化方法可以实现由抽象到具体的可重复的设计，可以实现大量的应用系统。
五、拓扑优化方法的应用拓扑优化方法目前被广泛应用在机械、航空航天、汽车等机械工程领域，具体应用包括但不限于：机械手和夹具的设计优化，汽车机架优化，电器结构优化，机械外壳优化，振动优化，和结构强度优化等等。
3188-1-1.htmlSigmund教授所编写的top优化经典99行程序，可以说是我们拓扑优化研究的基础；每一个新手入门都会要读懂这个程序，才能去扩展，去创新；99行程序也有好多个版本，用于求解各种问题，如刚度设计、柔顺机构、热耦合问题，但基本思路大同小异；本文拟对其中的一个版本进行解读，愿能对新手有点小小的帮助。
不详之处，还请论坛内高手多指点读懂了该程序，只能说是略懂拓扑优化理论了，我手里就有一些水平集源程序是成千上万行，虽然在99行的基础上成熟了很多，但依然还有很多的发展空间。
源程序如下：%%%% A 99 LINE TOPOLOGY OPTIMIZATION CODE BY OLE SIGMUND, JANUARY 2000 %%%%%%% CODE MODIFIED FOR INCREASED SPEED, September 2002, BY OLE SIGMUND %%%function top(nelx,nely,volfrac,penal,rmin);nelx=80;nely=20;volfrac=0.4;penal=3;rmin=2;% INITIALIZEx(1:nely,1:nelx) = volfrac;loop = 0;change = 1.;% START ITERATIONwhile change > 0.01loop = loop + 1;xold = x;% FE-ANAL YSIS[U]=FE(nelx,nely,x,penal);% OBJECTIVE FUNCTION AND SENSITIVITY ANAL YSIS[KE] = lk;c = 0.;for ely = 1:nelyfor elx = 1:nelxn1 = (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% FILTERING OF SENSITIVITIES[dc] = check(nelx,nely,rmin,x,dc);% DESIGN UPDA TE BY THE OPTIMALITY CRITERIA METHOD[x] = OC(nelx,nely,x,volfrac,dc);% PRINT RESULTSchange = 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 )])% PLOT DENSITIEScolormap(gray); imagesc(-x); axis equal; axis tight; axis off;pause(1e-6);end%%%%%%%%%% OPTIMALITY CRITERIA UPDATE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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;elsel2 = lmid;endend%%%%%%%%%% MESH-INDEPENDENCY FILTER %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [dcn]=check(nelx,nely,rmin,x,dc)dcn=zeros(nely,nelx);for i = 1:nelxfor j = 1:nelysum=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);endenddcn(j,i) = dcn(j,i)/(x(j,i)*sum);endend%%%%%%%%%%FE-ANAL YSIS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%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:nelxfor ely = 1:nelyn1 = (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% DEFINE LOADS AND SUPPORTS (HALF MBB-BEAM)F(2*(nelx/2+1)*(nely+1),1) = 1;fixeddofs = [2*(nely/2+1),2*nelx*(nely+1)+2*(nely/2+1)];alldofs = [1:2*(nely+1)*(nelx+1)];freedofs = setdiff(alldofs,fixeddofs);% SOLVINGU(freedofs, ：）= K(freedofs,freedofs) \ F(freedofs,：）;U(fixeddofs,：）= 0; % 这两行复制后应换成英文字符，我这里为了防止生成QQ表% 情修改了一下格式%%%%%%%%%% ELEMENT STIFFNESS MATRIX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [KE]=lkE = 1.;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)];% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%程序执行方法：（just for matlab new users）打开matlab，点开new M-file，将上述源程序复制粘贴到M-文件中，修改蓝色部分的格式，保存。
3188-1-1.htmlSigmund 教授所编写的top 优化经典99 行程序，可以说是我们拓扑优化研究的基础；每一个新手入门都会要读懂这个程序，才能去扩展，去创新；99 行程序也有好多个版本，用于求解各种问题，如刚度设计、柔顺机构、热耦合问题，但基本思路大同小异；本文拟对其中的一个版本进行解读，愿能对新手有点小小的帮助。
不详之处，还请论坛内高手多指点读懂了该程序，只能说是略懂拓扑优化理论了，我手里就有一些水平集源程序是成千上万行，虽然在99 行的基础上成熟了很多，但依然还有很多的发展空间。
源程序如下：%%%% A 99 LINE TOPOLOGY OPTIMIZATION CODE BY OLE SIGMUND, JANUARY 2000 %%%%%%% CODE MODIFIED FOR INCREASED SPEED, September 2002, BY OLE SIGMUND %%%function top(nelx,nely,volfrac,penal,rmin);nelx=80;nely=20;volfrac=0.4;penal=3;rmin=2;% INITIALIZEx(1:nely,1:nelx) = volfrac;loop = 0;change = 1.;% START ITERATIONwhile change > 0.01loop = loop + 1;xold = x;% FE-ANAL YSIS[U]=FE(nelx,nely,x,penal);% OBJECTIVE FUNCTION AND SENSITIVITY ANAL YSIS[KE] = lk;c = 0.;for ely = 1:nelyfor elx = 1:nelxn1 = (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)Ape nal*Ue'*KE*Ue;dc(ely,elx) = -pe nal*x(ely,elx)A(pe nal-1)*Ue'*KE*Ue;endend% FILTERING OF SENSITIVITIES[dc] = check(nelx,nely,rmin,x,dc);% DESIGN UPDA TE BY THE OPTIMALITY CRITERIA METHOD[x] = OC(nelx,nely,x,volfrac,dc);% PRINT RESULTSchange = 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 )])% PLOT DENSITIEScolormap(gray); imagesc(-x); axis equal; axis tight; axis off;pause(1e-6);end%%%%%%%%%% OPTIMALITY CRITERIA UPDATE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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;elsel2 = lmid;endend%%%%%%%%%% MESH-INDEPENDENCY FILTER %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [dcn]=check(nelx,nely,rmin,x,dc)dcn=zeros(nely,nelx);for i = 1:nelxfor j = 1:nelysum=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)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);endend%%%%%%%%%%FE-ANAL YSIS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%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:nelxfor ely = 1:nelyn1 = (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)A pe nal*KE;endend% DEFINE LOADS AND SUPPORTS (HALF MBB-BEAM)F(2*(nelx/2+1)*(nely+1),1) = 1;fixeddofs = [2*(nely/2+1),2*nelx*(nely+1)+2*(nely/2+1)];alldofs = [1:2*(nely+1)*(nelx+1)];freedofs = setdiff(alldofs,fixeddofs);% SOLVINGU(freedofs, ： )= K(freedofs,freedofs) \ F(freedofs, ： );U(fixeddofs, ： ) = 0; % 这两行复制后应换成英文字符，我这里为了防止生成QQ 表% 情修改了一下格式%%%%%%%%%% ELEMENT STIFFNESS MATRIX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function [KE]=lkE = 1.;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-nuA2)*[ 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)]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%程序执行方法： (just for matlab new users )打开matlab，点开new M-file，将上述源程序复制粘贴到M-文件中，修改蓝色部分的格式，保存。
