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

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

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 %

9-拓扑优化方法PPT课件

➢ 对这若干个子设计区域进行结构分析和灵敏度分析，建立设计变量与结构位移、应力、频率等关系，从而形成目标函数和约束条件；
➢ 按某种优化策略和准则从这若干个子设计区域中删除某些单元，用保留下来的单元描述结构的最优拓扑。
16
四、拓扑优化方法分类
从其物理模型的描述方法上一般分为 ➢ 基结构法(The Ground Structural Method) ➢ 均匀化方法(The Homogenization Method) ➢ 渐进结构优化方法 (The Evolutionary Structural
x fi u pu g xu i 0 i=1,2, ,n
由此可构造如下的迭代公式
x(k1) i
c(k)xi(k)
i=1,2,
,n
其中c(k)＝－1－ p
f u
ugxui
为小于1的因子
xi
7
x fi u pu g xu i 0 i=1,2, ,n
x fi u pu g xu i
i=1,2, ,n
对于由n个杆件组成的桁架结构，其满应力条件为
Fi Ai
i
i1,2,
,n
由此可构造如下的迭代公式
(k)
A(k1) i
i
A(k) i
i1,2,
,n
i
6
2. 基于K－T条件的准则法对于结构优化设计问题：
m in f(X ) X R n
s .t.g u ( X ) 0u 1 ,2 ,,p
极值点X*应满足的Kuhn－Tucker条件
结构优化与材料优化
第一节概述第二节结构优化设计的准则法第三节结构的拓扑优化方法第四节功能材料优化设计第五节柔性机构优化设计第六节结构多学科设计优化

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

拓扑优化PPT课件

N
KlineKe e1
（3）
其中[Ke]为单元刚度矩阵，下标lin表示它与设计变量是线性关系。
如果我们希望将问题设定成一个标准的嵌套式方程，其中要求平衡条
件能排除使刚度矩阵奇异的情况，可以用一个很小但非零的值 min来代替 0 ，其刚度矩阵可以写成：
N
K af f m in 1mineK e （4） e 1 16
这不仅意味着需要处理大量的设计变量，而且也影响到有限元分析的计算成本。些为了得到高精度的设计，运用模拟退火法、遗传算法、或是确定性方法计算成本都是很高的，而且这些方法只适用于相对较小的规模，或是些特定的设计问题，如最小柔顺性问题。
在式（2）的连续性问题假设中可以看出，寻求结构拓扑的基本思想是通过寻找
一个在定义域的子集上定义的指示函数来得到的，很明显这一问题很难解决，我
们可以通过限制子集的等级或是扩展设计集来获得一个适当的模式。对于柔度，均
匀的多尺度层状微结构组成了一个扩展的设计空间，同时也意味着整数约束松弛
为连续约束。
18
2.2 解决灰色尺度：差值模式
由于整数模将0
这包括大量的拓扑方面的著作特别是在所谓的均质化方法和多样性方法的方量的拓扑方面的著作特别是在所谓的均质化方法和多样性方法的方55三种优化的直观区别三种优化的直观区别66尺寸优化的设计变量是板的厚度二力杆的截面积以及梁截面的高度等尺寸优化的设计变量是板的厚度二力杆的截面积以及梁截面的高度等结构的尺寸参数尺寸优化的目的是要在满足结构的力学控制方程周长约束结构的尺寸参数尺寸优化的目的是要在满足结构的力学控制方程周长约束以及诸多性态约束条件的前提下寻求一组最优的结构尺寸参数使得关于结以及诸多性态约束条件的前提下寻求一组最优的结构尺寸参数使得关于结构性能的某种指标函数达到最优

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

声结构耦合系统双材料模型的拓扑优化设计商林源;赵国忠;陈刚【摘要】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.%基于微结构设计域法，并结合伴随法与放松形式准则法，研究了针对声结构耦合系统的双材料拓扑优化方法。

拓扑优化经典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行拓扑优化算法是一种经典的基于贪心策略的算法。

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

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

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

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

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

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

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

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

第三步：更新连接矩阵。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

（完整版）ANSYS拓扑优化原理讲解以及实例操作拓扑优化是指形状优化，有时也称为外型优化。

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

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

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

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

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

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

这些伪密度⽤PLNSOL ，TOPO 命令来绘出。

拓扑优化的⽬标——⽬标函数——是在满⾜结构的约束（V ）情况下减少结构的变形能。

减⼩结构的变形能相当于提⾼结构的刚度。

这个技术通过使⽤设计变量。

结构拓扑优化的基本思想是将寻求结构的最优拓扑问题转化为在给定的设计区域内寻求最优材料分布的问题。

通过拓扑优化分析，设计⼈员可以全⾯了解产品的结构和功能特征，可以有针对性地对总体结构和具体结构进⾏设计。

特别在产品设计初期，仅凭经验和想象进⾏零部件的设计是不够的。

只有在适当的约束条件下，充分利⽤拓扑优化技术进⾏分析，并结合丰富的设计经验，才能设计出满⾜最佳技术条件和⼯艺条件的产品。

连续体结构拓扑优化的最⼤优点是能在不知道结构拓扑形状的前提下，根据已知边界条件和载荷条件确定出较合理的结构形式，它不涉及具体结构尺⼨设计，但可以提出最佳设计⽅案。

拓扑优化技术可以为设计⼈员提供全新的设计和最优的材料分布⽅案。

拓扑优化基于概念设计的思想，作为结果的设计空间需要被反馈给设计⼈员并做出适当的修改。

最优的设计往往⽐概念设计的⽅案结构更轻，⽽性能更佳。

经过设计⼈员修改过的设计⽅案可以再经过形状和尺⼨优化得到更好的⽅案。

5.1.2优化拓扑的数学模型优化拓扑的数学解释可以转换为寻求最优解的过程，对于他的描述是：给定系统描述和⽬标函数，选取⼀组设计变量及其范围，求设计变量的值，使得⽬标函数最⼩（或者最⼤）。

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

基于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）。