拓扑优化经典99行

合集下载

拓扑优化_精品文档

-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),在该系统中小规模的计算使得很难利用刚体结构来实现铰链、轴承以及滑块处的机动性。
如果我们只考虑线性弹性材料（只发生微小变形）的分析问题，则决定输出位移的的有限元方法公式为：

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

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行程序解读

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

拓扑优化方法

结构尺寸优化设计
结构构型优化设计
结构形状优化设计
结构优化设计的数学描述
具有有限维的结构，其结构优化设计的数学模型的一般形式为
结构优化的目标函数
min f (X , Z, )
结构优化的约束条件
X 设计变量
Z
位移变量
频率变量
静力平衡条件
K(X )Z F(X )
固有频率条件
K(X )Y M (X )Y
wi
i
此时，平均频率的灵敏度计算公式为ຫໍສະໝຸດ m d2 ml
wi
l wi i
2
i1 i
d
其中
i
d
iT
K d
i
M d
i
i 1
此时，单元的刚度和质量的灵敏度计算公式为
Ke
d
BeT
De H d
Be
d
Ke d
e
d
NeT Ne
d
上式中 Ne 为单元形状函数,i为振型,为单元密度。
均匀化理论
连续体结构拓扑优化方法由于其优化模型描述方法的困难以及数值优化算法的巨大计算量而发展缓慢,其蓬勃发展的起点以 1988年kikuchi和bendsoe等人提出的均匀化算法(The Homogenization Method)为标志。
正是由于kikuchi和bendsoe的介绍后,拓扑优化方法在学术界得到了广泛地普及,并应用到材料设计、机构设计、ＭＥＭＳ器件设计、柔性微机构的设计和别的更复杂的结构设计中。
对于由n个杆件组成的桁架结构，其满应力条件为
Fi Ai
i
i 1, 2, , n
由此可构造如下的迭代公式
A(k1)
i

声结构耦合系统双材料模型的拓扑优化设计

声结构耦合系统双材料模型的拓扑优化设计商林源;赵国忠;陈刚【摘要】A bi-material topology optimization approach was investigated based on the microstructure-based design domain method with the combination of an adjoint method and a relaxed form of optimality criteria.An adjoint sensitivity method whereby sound pressure level derivative with respect to topology variables was proposed and deduced.A relaxed form of optimality criteria was deduced and used to solve the optimization problem of the coupled systems.Numerical examples show the high efficiency and the high accuracy of the adjoint sensitivity analysis,and the quick convergence and the high stability of the relaxed form of optimality criteria.The results also show that the topology optimization method of bi-material reduces internal noise and validates the optimization method.%基于微结构设计域法，并结合伴随法与放松形式准则法，研究了针对声结构耦合系统的双材料拓扑优化方法。

拓扑优化

¾ ¾
¾ ¾
TOSCA应用的常见问题
¾
¾
¾
TOSCA是否支持组合件结构拓扑优化？ TOSCA通过提取结构分析结果进行优化求解，从原则上说，任何支持的结构分析类型结构形式下的拓扑优化功能均能完成。对于软件系统来说，组合件和单个零部件结构表现为刚度阵规模大小不同，故而TOSCA支持组合件结构拓扑优化。 TOSCA软件的学习方式？ TOSCA.help是最好的学习教程，帮助内容包括安装方式和支持的求解器，命令行解释，算法适用范围，学习算例等。 TOSCA优化结果能否导入到CAD软件中？ TOSCA.Smooth模块对结果进行平滑处理并可生成igs，stl等几何格式，可导入到CAD软件模型中进行处理。
TOSCA应用的常见问题
¾
TOSCA是否支持强度优化？应力具有局部性特点，即在强度较弱的结构部位进行加强即可满足工程要求，通常工程优化流程为在得到刚度最大化结构的基础上进一步进行优化设计以满足应力约束要求。同时刚度最大化结果可理解为“全局应力”最小化结构，强度拓扑优化结果通常与刚度拓扑优化结果较大的相似性。 TOSCA 拓扑优化模块不具备强度优化功能，而在形状优化中设置应力约束条件。
( SIMP − model)
连续体结构拓扑优化建模方式
体积比约束：建模方式对应的工程含义为“结构体积一定的情况下，目标函数如柔顺度、固有频率值最小或最大”。建模方法的好处在于：由于体积敏度值保持恒正（或恒负），有利于优化求解，通常情况下优化迭代收敛性好。体积最小化目标: 建模方式对应的工程含义为“在满足结构响应量如节点位移、固有频率值不大于（或不小于）许用值的情况下，结构体积最小化”。建模方法的好处在于：具有处理多种不同性质约束的能力。

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

拓扑优化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)）。

拓扑优化99行算法解读

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

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

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

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

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

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

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

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

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

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

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

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

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

第三步：更新连接矩阵。

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

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

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

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

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

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

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

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

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

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

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

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

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

基于SIMP方法的最小柔度拓扑优化问题的优化准则算法研究

基于SIMP方法的最小柔度拓扑优化问题的优化准则算法研究李颖雄【摘要】拓扑优化的SIMP（Solid Isotropic Material Penalization）模型是对拓扑优化均匀化方法的改进与简化，假设材料密度在单元内是常数并以其为设计变量，而材料特性用单元密度的指数函数来模拟，从而减少了设计变量的个数，提高了计算效率。

本文主要是通过一些算例对基于SIMP方法的最小柔度拓扑优化问题的优化准则算法进行研究。

【期刊名称】《科技风》【年(卷),期】2011(000)015【总页数】1页(P93-93)【关键词】SIMP;拓扑优化;99行Matlab程序;准则优化法;最小柔度【作者】李颖雄【作者单位】华南理工大学,广东广州510000【正文语种】中文【中图分类】O342固体各向同性材料惩罚模型（SIMP:Solid Isotropic Material with Penalization）是一种常用的密度-刚度插值模型，该模型假设材料密度在单元内是常数并以其为设计变量，材料特性用单元密度的指数函数来模拟，对拓扑优化均匀化方法进行改进与简化，以达到方便计算，提高效率的目的。

SIMP模型认为材料的弹性模量与构成材料的密度呈指数关系，不需要均匀化方法所需的细观结构计算宏观物理、力学参量，其设计变量数目一般等价于结构的有限单元数目。

本文的主要内容为基于SIMP材料插值理论，在推导最小柔度拓扑优化问题的优化准则算法中，使用M atlab程序进行研究分析。

式中，ρ表示设计变量密度的列向量；C（ρ）表示目标函数；U，F，K分别表示全局位移向量，全局荷载向量，总刚度矩阵；P表示惩罚因子。

以悬臂梁结构为例，假定材料为各向同性：杨氏模量E=200N/mm，泊松比μ=0.3，厚度=1mm。

梁右下部受垂直向下作用力F=1N，模型近似为二维平面应力问题。

结构离散为60×45个四节点四边形单元，30%体积约束。

惩罚因子P取不同数值时的拓扑优化结果，网格过滤参数rmin取1.0由计算结果可知：当P≤2.5时存在大量多孔材料，计算结果没有可制造性。