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

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

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

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