拓扑优化经典99行程序解读

合集下载

北航拓扑优化程序学习报告

拓扑优化的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带惩罚因子的各项同性材料模型法），该方法是拓扑优化中常用的变密度材料插值模型中最具代表性的一种。

拓扑优化_精品文档

-1整数变量问题变为0～1间的连续变量优化模型，获得方程（在设计变
量上松弛整数约束）的最直接方式是考虑以下问题：
min u,
uout
N
s.t.: min 1 min e Ke u f e1
N
vee V
e1
0 e 1, e 1,2,, N
其中 e 可取0-1之间的值
（6）
然而这种方程会导致较大区域内 e 是在0-1之间的值，所以必须添加额外的约束来避免这种“灰色”区域。要求是优化结果基本上都在 e 1 或
而对于结构拓扑优化来说，其所关心的是离散结构中杆件之间的最优连接关系或连续体中开孔的数量及位置等。拓扑优化力图通过寻求结构的最优拓扑布局(结构内有无孔洞，孔洞的数量、位置、结构内杆件的相互联接方式)，使得结构能够在满足一切有关平衡、应力、位移等约束条件的情形下，将外荷载传递到支座，同时使得结构的某种性能指标达到最优。拓扑优化的主要困难在于满足一定功能要求的结构拓扑具有无穷多种形式，并且这些拓扑形式难以定量的描述即参数化。
结构渐进优化法(简称ESO法)
通过将无效的或低效的材料一步步去掉，获得优化拓扑，方法通用性好，可解决尺寸优化，还可同时实现形状与拓扑优化(主要包括应力，位移／刚度和临界应力等约束问题的优化)。
2.问题的设定
柔顺机构的拓扑优化
首先假设线性弹性材料有微小的变形
柔顺结构的一个重要运用在于机电系统（MicroElectroMechanical Systems(MEMS),在该系统中小规模的计算使得很难利用刚体结构来实现铰链、轴承以及滑块处的机动性。
如果我们只考虑线性弹性材料（只发生微小变形）的分析问题，则决定输出位移的的有限元方法公式为：

(完整版)ANSYS拓扑优化原理讲解以及实例操作

拓扑优化是指形状优化，有时也称为外型优化。

拓扑优化的目标是寻找承受单载荷或多载荷的物体的最佳材料分配方案。

这种方案在拓扑优化中表现为“最大刚度”设计。

与传统的优化设计不同的是，拓扑优化不需要给出参数和优化变量的定义。

目标函数、状态变量和设计变量（参见“优化设计”一章）都是预定义好的。

用户只需要给出结构的参数（材料特性、模型、载荷等）和要省去的材料百分比。

给每个有限元的单元赋予内部伪密度来实现。

这些伪密度用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 -设计变量。

拓扑优化学习报告_北理工_王路

i 1 j 1 i 1 j 1
s.t
KU F
V fV0 xi , j i , j
i 1 j 1 m n
0 xmin xi , j xmax 1
其中：
X ——单元相对密度矢量
C ——结构的柔度
F ——载荷矢量
U ——位移矢量
北京理工大学车辆工程王路
FILTER %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% function [dcn]=check(nelx,nely,rmin,x,dc) dcn=zeros(nely,nelx); for i = 1:nelx 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
E( xi ) Emin ( xi ) p ( E0 Emin ), xi 0,1
其中：
E(xi ) ——插值以后的弹性模量 E0 ——实体部分材料的弹性模量 Emin ——孔洞部分材料的弹性模量
（1）
xi ——单元相对密度，取值为
1 时表示有材料，为 0 时表示无材料即孔洞
p ——惩罚因子
（二）拓扑优化问题的描述
（1）材料插值模型

北航拓扑优化程序学习报告

拓扑优化的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 带惩罚因子的各项同性材料模型法），该方法是拓扑优化中常用的变密度材料插值模型中最具代表性的一种。

拓扑优化经典99行

% Please sent your comments to the author: [email]sigmund@fam.dtu.dk[/email]%
% %
%%%% A 99 LINE TOPOLOGY OPTIMIZATION CODE BY OLE SIGMUND, JANUARY 2000 %%%
%%%% CODE MODIFIED FOR INCREASED SPEED, September 2002, BY OLE SIGMUND %%%
% DESIGN UPDATE BY THE OPTIMALITY CRITERIA METHOD
[x] = OC(nelx,nely,x,volfrac,dc);
% PRINT RESULTS
change = max(max(abs(x-xold)));
disp([' It.: ' sprintf('%4i',loop) ' Obj.: ' sprintf('%10.4f',c) ...
K(edof,edof) = K(edof,edof) + x(ely,elx)^penal*KE;
end
end
% DEFINE LOADS AND SUPPORTS (HALF MBB-BEAM)
F(2,1) = -1;
fixeddofs = union([1:2:2*(nely+1)],[2*(nelx+1)*(nely+1)]);
% The author reserves all rights but does not guaranty that the code is %

拓扑优化经典99行程序解读

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-文件中，修改蓝色部分的格式，保存。

Sigmund_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 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行拓扑优化程序完全注释

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时进⾏密度滤波。

拓扑优化算法及其实现

>Imagesc(-x)
敏度分析
目的：消除棋盘格效应及网格依赖性
方法：1、高阶单元法（计算量大） 2、周长约束法（周长约束的上限值需要依靠经验来确定，因为局部尺寸和周长边
界间没有直接的关系。如果周长约束边界定得太紧，则可能导致没有计算结果，如果定的太松又达不到预期的效果。因此约束边界很难确定，这种情况在三维问题下特别明显）
内单元的应变能密度是恒定的常数。由此建立更新设计变量迭代格式：
n1enB,
Bp( e)p1ue Tk0ue ve
n 1 m e n a B x (m in,
e nm ) if
if e nB m ax (m in, e nm ) m ax (m in, e nm )e nB m in (1 , e nm )
RAMP（ Rational Approximation of Material Properties ）（材料属性的理性近似模型）
Level Set法 (水平集法)
ICM(独立映射法)
ESO(进化法)
……
优化求解方法
OC法(优化准则法) MMA法(移动渐进线法) SLP(序列线性规划法） SQP(序列二次规划法) …………
优化求解 OC法优化求解
基本原理变量约束优化的KT条件：定义在闭区间上的一元函数的优化（仅有变量上下限约束）
基本原理
OC法数学模型对应于目标函数的拉格朗日函数为：
OC法数学模型
OC法数学模型
u
优化设计准则
OC法数学模型
优化设计准则
上式即为设计变量的迭代准则。由该式可以看出当柔度取得极值时,在整个设计区域
解决方法：线性和非线性 Pratt,(1991) 证明了非线性技术优于线性技术。但是非线性技术不能应用于拓扑优化，因为非线性方法常会造成拓扑优化问题的不光滑。

拓扑优化经典99行程序解读

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-文件中，修改蓝色部分的格式，保存。

单边缆车

拓扑优化设计:单边缆车1.OC 法简介：在机械优化设计中常常用到拓扑优化设计中的OC 法,拓扑优化是在给定的设计域，约束和载荷条件下，以结构的某种性能指标达到最优为目标确定材料的拓扑结构，确定结构构件的连接方式，结构内有无空洞、空洞数量及位置等拓扑形式，是结构轻量化设计的关键技术。

Bendsoe 在 1989 年提出的（Solid Isotropic Material with Penalization ）SIMP 法，以单元密度为设计变量，通过引入惩罚因子对中间密度值进行惩罚，而使单元材料接近 0 或 1，材料为 0 的单元予以删除，材料为 1 的单元则予以保留，最后得到材料的拓扑结构。

拓扑优化的求解方法很多比如（变密度法、水平集法、独立映射法、进化法），但是常用的是优化准则法（Optimality Criteria method: OC 法）。

常用的优化求解方法有OC 法，即优化准则法、MMA 法，即移动渐近线法、SLP 法，即序列线性规划法、SQP ，即序列二次规划法。

OC 法是一种启发式优化算法，并根据数学规划理论中的 Kuhn Tucker 条件导出满足条件的优化解。

传统的 OC 法得到的拓扑结构中往往存在一定的中间密度单元和棋盘格现象，如图 1 所示。

图 1 拓扑优化中的数值不稳定现象O.Sigmund 对拓扑优化中的数值不稳定问题进行了综述，其中过滤法是消除数值不稳定现象的一种有效方法，过滤因子的选取对拓扑结果的影响较大，过大或过小都无法得到理想的拓扑结构。

为了得到比较光滑的拓扑结构，文献利用节点变量和投影函数进行边界的光滑处理，但结果中仍存在一定比例的中间密度单元。

在得到较理想的拓扑后，提高 SIMP 法的惩罚因子对消除中间密度单元影响不大。

为了提高惩罚效率，得到更为理想的拓扑结构，本文对 OC 法优化求解得到的密度进行惩罚。

本文通过最小柔度设计，以 SIMP 法为基础，引入惩罚因子对 OC 法中的单元密度进行惩罚，以节点密度为设计变量，建立过滤因子和单元大小的关系，对99 行程序进行适当的改变进行结构优化设计。

拓扑优化99行算法解读

拓扑优化99行算法解读拓扑优化是一种解决网络优化问题的算法。

该算法通过对给定网络中各节点之间的关系进行建模，来确定节点之间的最优连接方式，以使得整个网络的性能达到最优化。

拓扑优化算法的主要目标是通过重新组织网络节点之间的连接关系，来最小化网络中的延迟、能耗或其他性能指标。

这种算法可以应用于各种不同类型的网络，包括计算机网络、通信网络以及物联网等。

99行拓扑优化算法是一种经典的基于贪心策略的算法。

该算法首先根据网络中各节点之间的距离来初始化一个连接矩阵。

然后，通过遍历网络中的各个节点，依次选择最佳的连接方式来更新连接矩阵。

具体的算法步骤如下：第一步：初始化连接矩阵。

对于给定的网络，可以通过计算各个节点之间的距离来初始化连接矩阵。

距离可以根据节点之间的物理距离或其他性能指标来计算。

第二步：对网络中的各个节点进行遍历。

对于每个节点，计算其与其他节点之间的连接费用，并选择最佳的连接方式。

连接费用可以根据节点之间的距离、带宽或其他性能指标来计算。

第三步：更新连接矩阵。

在选择了最佳连接方式后，更新连接矩阵中对应节点之间的连接信息。

第四步：重复以上步骤，直到所有节点都被遍历完毕。

在进行拓扑优化时，算法会优先选择具有低延迟、高带宽或其他性能指标的连接方式。

这样可以使得整个网络的性能达到最佳化。

拓扑优化算法的优点在于它具有较好的可伸缩性和灵活性。

算法可以应用于各种不同规模的网络，并且可以根据不同的性能指标进行优化。

此外，算法的实现相对简单，计算复杂度较低，因此适用于实际应用。

然而，拓扑优化算法也存在一些限制。

首先，算法可能会导致网络中的某些节点之间的连接断开，从而影响网络的连通性。

其次，算法没有考虑节点之间的地理位置和网络负载等因素，这可能会导致某些节点之间的连接过于拥挤。

总之，拓扑优化算法是一种解决网络优化问题的经典算法。

它通过重新组织网络节点之间的连接关系，来最小化网络中的延迟、能耗或其他性能指标。

虽然算法存在一些限制，但其简单的实现和较低的计算复杂度使其成为一种实际应用广泛的算法。

拓扑优化学习报告-北理工-王路

基于99行程序的拓扑优化学习报告（一）背景和前言随着汽车工业的飞速发展以及日益突出的能源问题，汽车工业面临的挑战以及竞争环境也越来越激烈，对汽车产品提出了降低其制造成本及燃油经济性的新要求。

在提高汽车安全性、减少汽车排放和解决能源消耗的背景下，提出了汽车轻量化技术。

实现汽车轻量化的途径包括三个方面：结构优化技术、新型材料和先进性制造工艺。

其中，我们所讨论的是结构优化技术，其中结构优化设计分为三个层次：尺寸优化（Size Optimization）、形状优化（Shape Optimization）和拓扑优化（Topology Optimization）。

本文我们基于99行matlab程序初步学习拓扑优化技术中的理论和优化方法。

拓扑优化技术指的是在给定的设计空间内寻求最佳的材料布局，同时在满足平衡方程、物理关系、几何关系和边界约束条件下使得结构达到某种性能最优的应用技术。

拓扑优化的理论研究最早可以追溯到Michel提出的桁架理论，连续体结构的拓扑优化由于描述和数值计算得困难，发展一直相对缓慢，直到Bendsoe和Kikuchi在1988年提出的均匀化方法之后才得到迅速的发展，其基本思想是在组成拓扑结构中引入微结构，通过微结构的几何参数作为设计变量，通过微结构的增加和删减实现结构的拓扑形状的改变，实现拓扑优化和尺寸优化的统一。

在微结构的基础上，我们介绍变密度法的应用，变密度法是在均匀化方法的基础上产生的，把材料引入微结构代之以密度在0~1之间变化的假想材料，把密度作为设计变量，从而实现材料的删减，因其模型简单、计算变量相对较少成为目前广泛采用的方法。

根据不同的插值模式，变密度法又有不同的插值模型：SIMP法（Solid Isotropic Material with Penalization）、Hashin-Shtrikman法，以及RAMP法（Rational Approximation of Material Properties）。

拓扑优化经典99行程序解读

拓扑优化经典99行程序解读3188-1-1.htmlSigmund教授所编写的top优化经典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-ANALYSIS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%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-文件中，修改蓝色部分的格式，保存。

拓扑优化99行代码翻译

拓扑优化中的99行matlab代码——o.sigmund摘要这篇文章描述了用matlab语言来简洁的实现在静态负载下符合最小化原理的拓扑优化。

总共只需要输入99行代码，包括优化程序和有限元分析子程序。

这99行代码中，其中36行为主程序，12行为基于最优控制器的优化程序，16行为敏度过滤分析，其余35行代码作为有限元分析。

实际上，除去注释行以及输出行、有限元分析行，仅有49行matlab代码输入用于解决一个适定的拓扑优化问题。

再加上3行补充程序代码，这个程序就可以解决多种负载工况问题。

这个代码主要是以教育指导为用，完整的matlab代码在附录中给出，同时也可以在网页http://www.topopt.dtu.dk上下载。

关键词拓扑优化教育最优准则万维网matlab代码1 简介文中展示的matlab代码主要是为工程教育所用。

学生和在拓扑领域的新手可以在网页http://www.topopt.dtu.dk上下载。

这个代码可以用于结构最优化课程学习，学生们可以在多重负载工况、独立网格选择策略、无源场进行扩展应用。

另一种可能就是用来激发学生们的直觉来进行最优化设计。

研究生可以推测探究在给定边界条件和容量的情况下的拓扑优化并、比较得出最优策略。

在文献中，你可以找到很多处理拓扑优化问题的方法。

在一篇Bendsøe and Kikuchi (1988)的原创论文中，基于对现存结论的学习，所谓微观结构或均化作用的方法被使用。

均化作用方法在很多文章中都被采用，但它也存在一些缺点，比如对微观结构最优化方法果断的评估与决策很麻烦的，而且结果很难获得如果没有对微观结构进行确定的长度衡量。

然而，在这个意义上来说均化作用方法对拓扑优化也是很重要的，它可以在结构的理论分析上提供一定的界限。

拓扑优化的另一种方法叫做“幂律法则”或者SIMP法（(Solid Isotropic Material with Penalization) (Bendsøe1989; Zhou and Rozvany 1991; Mlejnek 1992)）。