下载文档原格式
拓扑优化的99行程序学习报告4月19日2011《结构优化设计》课程学习报告任课教师：李书一、前言：在最近的结构优化设计课程上学习了O.Sigmund的《A 99 line topology optimization code written in Matlab》一文，对拓扑优化的理论原理与实际的计算机程序实现都有了一定的理解，文章主要是通过拓扑优化的原理来实现对简单结构的静力学问题的优化求解，而编写的代码仅有99行，包括36行的主程序，12行的OC优化准则代码，16行的网格过滤代码和35行的有限元分析代码。
自1988 年丹麦学者Bendsoe与美国学者Kikuchi提出基于均匀化方法的结构拓扑优化设计基本理论以来,均匀化方法应用到具有周期性结构的材料分析中,近几年该方法已经成为分析夹杂、纤维增强复合材料、混凝土材料等效模量,以及材料的细观结构拓扑优化常用的手段之一。
其基本思想是在组成拓扑结构的材料中引入微结构，优化过程中以微结构的几何尺寸作为设计变量，以微结构的消长实现其增删，并产生介于由中间尺寸微结构组成的复合材料，从而实现了结构拓扑优化模型与尺寸优化模型的统一。
文章就是通过均匀化的基础，结合拓扑结构优化的工程实际，以计算机模拟的方法将拓扑优化的一般过程呈现出来，有助于初涉拓扑优化的读者对拓扑优化有个基础的认识。
二、拓扑优化问题描述为了简化问题的描述，文中假设设计域是简单的矩形形式，且在进行有限元离散的时候采用正方形单元对其进行离散。
这样不仅便于进行单元离散和单元编号，也利于对结构进行几何外形的描述。
一般说来，基于指数逼近法的拓扑优化最小化的问题可作如下描述：文中采用的对结构材料属性的描述是所谓的“指数逼近法”或者称为SIMP 逼近法，即（Solid Isotropic Material with Penalization带惩罚因子的各项同性材料模型法），该方法是拓扑优化中常用的变密度材料插值模型中最具代表性的一种。
如何利用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的单元，才执行优化计算；否则，就不执行优化计算。
拓扑优化是指形状优化，有时也称为外型优化。
拓扑优化的目标是寻找承受单载荷或多载荷的物体的最佳材料分配方案。
这种方案在拓扑优化中表现为“最大刚度”设计。
与传统的优化设计不同的是，拓扑优化不需要给出参数和优化变量的定义。
目标函数、状态变量和设计变量（参见“优化设计”一章）都是预定义好的。
用户只需要给出结构的参数（材料特性、模型、载荷等）和要省去的材料百分比。
给每个有限元的单元赋予内部伪密度来实现。
这些伪密度用PLNSOL ，TOPO 命令来绘出。
拓扑优化的目标——目标函数——是在满足结构的约束（V ）情况下减少结构的变形能。
减小结构的变形能相当于提高结构的刚度。
这个技术通过使用设计变量。
结构拓扑优化的基本思想是将寻求结构的最优拓扑问题转化为在给定的设计区域内寻求最优材料分布的问题。
通过拓扑优化分析，设计人员可以全面了解产品的结构和功能特征，可以有针对性地对总体结构和具体结构进行设计。
特别在产品设计初期，仅凭经验和想象进行零部件的设计是不够的。
只有在适当的约束条件下，充分利用拓扑优化技术进行分析，并结合丰富的设计经验，才能设计出满足最佳技术条件和工艺条件的产品。
连续体结构拓扑优化的最大优点是能在不知道结构拓扑形状的前提下，根据已知边界条件和载荷条件确定出较合理的结构形式，它不涉及具体结构尺寸设计，但可以提出最佳设计方案。
拓扑优化技术可以为设计人员提供全新的设计和最优的材料分布方案。
拓扑优化基于概念设计的思想，作为结果的设计空间需要被反馈给设计人员并做出适当的修改。
最优的设计往往比概念设计的方案结构更轻，而性能更佳。
经过设计人员修改过的设计方案可以再经过形状和尺寸优化得到更好的方案。
5.1.2优化拓扑的数学模型优化拓扑的数学解释可以转换为寻求最优解的过程，对于他的描述是：给定系统描述和目标函数，选取一组设计变量及其范围，求设计变量的值，使得目标函数最小（或者最大）。
一种典型的数学表达式为：()()()12,,0,,0min ,g x x v g x x v f x v ⎧=⎪⎪≤⎨⎪⎪⎩式中，x -系统的状态变量；12g g 、-一等式和不等式的结束方程；(),f x v -目标函数；v -设计变量。
拓扑优化的99行程序学习报告4月19日2011《结构优化设计》课程学习报告任课教师：李书一、 前言：在最近的结构优化设计课程上学习了O.Sigmund 的《A 99 line topology optimization code written in Matlab 》一文，对拓扑优化的理论原理与实际的计算机程序实现都有了一定的理解，文章主要是通过拓扑优化的原理来实现对简单结构的静力学问题的优化求解，而编写的代码仅有99行，包括36行的主程序，12行的OC 优化准则代码，16行的网格过滤代码和35行的有限元分析代码。
自1988 年丹麦学者Bendsoe 与美国学者Kikuchi 提出基于均匀化方法的结构拓扑优化设计基本理论以来,均匀化方法应用到具有周期性结构的材料分析中,近几年该方法已经成为分析夹杂、纤维增强复合材料、混凝土材料等效模量,以及材料的细观结构拓扑优化常用的手段之一。
其基本思想是在组成拓扑结构的材料中引入微结构，优化过程中以微结构的几何尺寸作为设计变量，以微结构的消长实现其增删，并产生介于由中间尺寸微结构组成的复合材料，从而实现了结构拓扑优化模型与尺寸优化模型的统一。
文章就是通过均匀化的基础，结合拓扑结构优化的工程实际，以计算机模拟的方法将拓扑优化的一般过程呈现出来，有助于初涉拓扑优化的读者对拓扑优化有个基础的认识。
二、 拓扑优化问题描述为了简化问题的描述，文中假设设计域是简单的矩形形式，且在进行有限元离散的时候采用正方形单元对其进行离散。
这样不仅便于进行单元离散和单元编号，也利于对结构进行几何外形的描述。
一般说来，基于指数逼近法的拓扑优化最小化的问题可作如下描述：01min :()()():0min 1NT p Te e e Xe c X U KU x u k u V X subjecttof V KU FX x =⎫==⎪⎪⎪=⎬⎪⎪=⎪<≤≤⎭∑文中采用的对结构材料属性的描述是所谓的“指数逼近法”或者称为SIMP逼近法，即（Solid Isotropic Material with Penalization 带惩罚因子的各项同性材料模型法），该方法是拓扑优化中常用的变密度材料插值模型中最具代表性的一种。
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-文件中，修改蓝色部分的格式，保存。
%%%% 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 ITERA TIONwhile 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 CRITERIAUPDATE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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-INDEPENDENCYFILTER %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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;%%%%%%%%%% ELEMENT STIFFNESSMATRIX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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)];%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%。
TopOpt针对99⾏改进的88⾏拓扑优化程序完全注释 博客搬家到⾃⼰搭建的 啦q(≧▽≦q)，⼤家快来逛逛鸭！The program can be promoted by line:top88(120,40,0.5,3.0,3.5,1)代码注释88⾏程序为了提⾼效率，逻辑、流程没有99⾏程序那么清楚，建议初学先读99⾏。
%%%%%%%%%%%% AN 88 LINE TOPOLOGY OPTIMIZATION CODE Nov 2010 %%%%%%%%%%%% %%%%%%%%%%%% COMMENTED - OUT BY HAOTIAN_W AUGUST 2020 %%%%%%%%%%%%function top88(nelx,nely,volfrac,penal,rmin,ft)% ===================================================================================% 88⾏程序在99⾏程序上的主要改进:% 1) 将for循环语句向量化处理，发挥MATLAB矩阵运算优势;% 2) 为不断增加数据的数组预分配内存，避免MATLAB花费额外时间寻找更⼤的连续内存块;% 3) 尽可能将部分程序从循环体⾥抽出，避免重复计算;% 4) 设计变量不再代表单元伪密度，新引⼊真实密度变量xphys;% 5) 将原先的所有⼦程序都集成在主程序⾥，避免频繁调⽤;% 总体上，程序的效率有显著提升（近百倍）、内存占⽤降低，但是对初学者来说可读性不如99⾏% ===================================================================================% nelx : ⽔平⽅向上的离散单元数;% nely : 竖直⽅向上的离散单元数;%% volfrac : 容积率，材料体积与设计域体积之⽐，对应的⼯程问题就是"将结构减重到百分之多少";%% penal : 惩罚因⼦，SIMP⽅法是在0-1离散模型中引⼊连续变量x、系数p及中间密度单元，从⽽将离% 散型优化问题转换成连续型优化问题，并且令0≤x≤1，p为惩罚因⼦，通过设定p>1对中间密% 度单元进⾏有限度的惩罚，尽量减少中间密度单元数⽬，使单元密度尽可能趋于0或1;%% 合理选择惩罚因⼦的取值，可以消除多孔材料，从⽽得到理想的拓扑优化结果:% 当penal<=2时存在⼤量多孔材料，计算结果没有可制造性;% 当penal>=3.5时最终拓扑结果没有⼤的改变;% 当penal>=4时结构总体柔度的变化⾮常缓慢，迭代步数增加，计算时间延长;%% rmin : 敏度过滤半径，防⽌出现棋盘格现象;% ft : 与99⾏程序不同的是，本程序提供了两种滤波器，% ft=1时进⾏灵敏度滤波，得到的结果与99⾏程序⼀样；ft=2时进⾏密度滤波。
