下载文档原格式
拓扑优化的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带惩罚因子的各项同性材料模型法），该方法是拓扑优化中常用的变密度材料插值模型中最具代表性的一种。
matlab拓扑优化99行代码以下是一份 MATLAB 拓扑优化的代码，该代码共99行，带有注释，共计1000字左右。
%% MATLAB 拓扑优化代码% 该代码实现了一种基于拓扑优化的算法，用于最小化结构体积并满足力学约束% 本代码适用于 MATLAB （版本 R2020a 或以上）。
使用前请确保已安装优化工具箱（Optimization Toolbox）。
% 作者：[作者姓名]% 时间：[编写时间]%% 1. 定义初始构型% 定义结构参数L = 1.0; % 长度H = 0.2; % 高度W = 0.01; % 宽度n_x = 30; % 网格数目（x-方向）n_y = 10; % 网格数目（y-方向）% 定义初始密度变量 rho（所有单元的密度都相等）rho_min = 1e-3; % 最小密度（设定为一个小值，以避免矩阵不可逆）rho = ones(n_x * n_y, 1) * 0.5; % 初始密度变量 rhorho(1:n_y) = 1.0; % 左侧固定边界rho((n_x-1)*n_y+1:n_x*n_y) = 1.0; % 右侧固定边界% 将 rho 转化为二维数组rho = reshape(rho, n_y, n_x)';rho = rho(:);% 定义单元尺寸变量 x 和 ydx = L / n_x;dy = H / n_y;x = repmat([0:dx:L]', n_y+1, 1);y = repmat([0:dy:H], n_x+1, 1)';% 绘制初始构型figurepcolor(x, y, reshape(rho, n_y, n_x))axis equalshading interpcolormap(gray)title('Initial Topology')%% 2. 定义边界和 FEA 参数% 定义边界和 FEA 网格xmin = 0; ymin = 0;xmax = L; ymax = H;nelx = n_x-1; nely = n_y-1;[xx,yy] = meshgrid(linspace(xmin,xmax,nx),linspace(ymin,ymax,ny));[II,JJ] = meshgrid(1:nelx,1:nely);% 定义单元编号矩阵（每个单元有8个节点）ind = reshape(1:(n_x+1)*(n_y+1), n_x+1, n_y+1);ind = ind(1:end-1, 1:end-1);eleInd = kron(ind, ones(8,1)) + kron(ones(size(ind)), [0 1 (n_x+1)*[1 1]+[0 1] (n_x+1)*(n_y+1)+[0 1] -n_x-1 -n_x]);%% 3. 定义拓扑优化参数% 定义拓扑优化参数volfrac = 0.45; % 允许的最大材料体积分数penal = 3; % Sigmoid 惩罚因子% 定义 Sigmoid 函数H = @(x, a) 1.0 ./ (1.0 + exp(-a .* x));% 定义材料与空气的弹性模量Emin = 1e-9; % 材料弹性模量E = H(rho, penal) * (E0-Emin) + Emin; % 节点弹性模量 % 定义 Sigmoid 求导矩阵H1rho = spdiags(H1(rho, penal), 0, dof, dof);% 定义全局刚度矩阵K = sparse(dof, dof);for el = 1:size(eleInd, 1)% 获取当前单元的节点编号nodes = eleInd(el,:);% 获取当前单元的节点坐标X = xy(nodes,:);% 计算单元刚度矩阵[Ke,~] = elestiff(X,E(nodes),nu);% 将单元刚度矩阵加入全局刚度矩阵K(nodes,nodes) = K(nodes,nodes) + Ke;end% 定义移动限制矩阵（仅移动 rho < 1 的节点）fixeddofs = union(zetax,zetay(:));freedofs = setdiff(1:dof,fixeddofs);% 定义循环参数iter = 1;change = 1.0;tol = 1e-6;maxIter = 100;rho0 = zeros(dof, 1);rho0(freedofs) = 1; % 初始密度为 1ksym = symrcm(K); % 使用对称逆消结构进行前置处理[L,U] = ilu(K(ksym,ksym)); % 使用不完全 LU 分解进行矩阵逆求解前置处理 while change > tol && iter <= maxIter% 保存当前的 rhorho_last = rho;% 使用 Sigmoud 函数更新单元弹性模量E = H(rho, penal) * (E0-Emin) + Emin;% 更新全局刚度矩阵K = sparse(dof, dof);for el = 1:size(eleInd, 1)% 获取当前单元的节点编号nodes = eleInd(el,:);% 获取当前单元的节点坐标X = xy(nodes,:);% 计算单元刚度矩阵[Ke,~] = elestiff(X,E(nodes),nu);% 将单元刚度矩阵加入全局刚度矩阵K(nodes,nodes) = K(nodes,nodes) + Ke;end% 移动节点[v,d] = eigs(K(freedofs,freedofs), 1, 'sm',struct('tol',1e-8,'maxit',500,'P',L(freedofs,freedofs),'D',U(freedofs,freedofs )));rho0(freedofs) = rho(freedofs) + v;rho = min(max(rho0, 0), 1);% 更新拓扑结构rho = update(rho, volfrac, x, y, n_x, n_y, H, H1, H1rho, Emin, penal);% 绘制拓扑结构figure(3)pcolor(reshape(rho, n_y, n_x)')shading interpaxis equaltitle(['Iteration: ', num2str(iter)])% 计算变化量change = norm(rho_last - rho);% 更新迭代计数器iter = iter + 1;end% 输出结果fprintf('Volume fraction: %f\n', sum(rho) / numel(rho)) fprintf('Compliance: %f\n', F' * full(K) * F)。
Simwe 仿真论坛---(邀请注册) » 我的资料用户组: 3级会员在线时间: 总计在线 159.33 小时, 本月在线 55.33 小时 (升级剩余时间 1 小时)积分: 8技术积分: 8 , 仿真币: 937信用信用评评价 (查看详细详细的信用的信用的信用记录记录) 买家信用评价: 0卖家信用评价: 0417332551(UID: 456465) 性别:男自我介绍:200 字节以内不支持自定义 Discuz! 代码QQ:417332551MSN:发帖数级别: 实习记者阅读权限: 15帖子: 523 篇(占全部帖子的 0.02%)平均每日发帖: 1.07 篇 精华: 0 篇 页面访问量: 3121注册日期: 2010-3-27 上次访问: 昨天 15:36 最后发表: 3 小时前 注册 IP: 121.8.210.55 - 上次访问IP: 58.217.103.11 -发短消息 加为好友 使用道具 搜索帖子 个人空间417332551 | 我的帖子 空间 短消息 | 个人中心 退出综合首页个人空间电子期刊技术论坛搜索插件帮助导航Google 搜索搜索WWW搜索百度搜索搜索WWW搜索 网站已站已经经在公安机在公安机关备关备关备案案并处于监督之中督之中，，请不要不要发发布违反法律的反法律的内内容。
违反者可能反者可能会会被追究刑事被追究刑事责责1、密码问题：有一些网友注册后，无法收到密码，原因可能是：a.没有填写正确的电子信箱或本来就是假电子信箱。
b.采用境外的 的信箱，密码信发出时可能被拒收。
c.登陆时密码没有区分大小写，例如密码中A 和a 是不同的。
d.修改密码只需点击右上角码”和“确认密码”两个选项改为一致的新密码即可。
2、登陆问题：a.有一些网友注册后无法登陆本站.很有可能是因为IE 的隐私政策设置问题.论坛需要cookies 支持，解决方法如下:把可.b.没有填入正确的安全码或者安全码已经过期。
图像压缩编码 image compression encoding推挽放大器 push-pull amplifier退n步arq go-back-n arq吞吐量 throughput脱网工作 talkaround拓扑学 topology拓扑优化 topological optimizationw外部抗扰性 external immunity外差跟踪 heterodyne tracking外差接收 heterodyne reception完成到示忙用户呼叫completion of calls to busy subscribers (ccbs)完成接续 complete connection完全故障 complete failure完全握手证实 full handshake authentication网管中心 network management center (nmc)网间互通 interworking between network网络保密 network security网络层 network layer网络传输时延 network transfer delay网络管理 network management (nm)网络管理要求 network management requirements网络互连 network interworking网络控制中心 network control center (ncc)网络口 port (of a network)网络体系结构 network architecture网络通路 network path网络拓扑 network topology网络行政管理 network administration网络拥挤 network congestion网络优化 network optimization网络指定判据 network directed criteria网络资源 network resource微巴 microbar微波 microwave微波暗室（无反射室） microwave unreflected chamber微波单片集成电路microwave monolithic integrated circuit (mmic)微波低躁声放大器 microwave low noise amplifier微波电路 microwave circuit微波功率放大器 microwave power amplifier微波混合集成电路 microwave hybrid integrated circuit微波集成电路 microwave integrated circuit微波通信 microwave communication微波中继站 microwave relay station微波终端站 microwave terminal station微处理机开发系统 microprocessor development system微电子技术 microelectronic technology微电子组装 microelectronic packaging微处理机 microprocessor微区 microcell微微区 picocell微特比算法 viterbi algorithm维修 maintenance维修层次 indenture level of maintenance维修性 maintainability伪比特 dummy bits伪三进制 pseudo-ternary signal伪随机码 pseudorandom code伪随机数 rand伪随机序列 pseudorandom sequence伪突发 dummy burst (db)尾部比特 tail bits卫星航空移动业务 aeronautical mobile-satellite service卫星间链路 inter-satellite link卫星陆地移动业务 land mobile-satellite service卫星水上移动业务 maritime mobile satellite service卫星网络 satellite network卫星移动通信系统 satellite mobile communication system卫星移动业务 mobile-satellite service卫星转发器 satellite transponder未来公用陆地移动通信系统future public land mobile telecommunication system (fplmts) 未完成呼叫 unsuccessful call位置撤消登记 location deregistration位置登记 location registration位置登记器 location register (lr)位置跟踪 location tracking位置更新规程 location updating procedure位置区 location registration zone位置区识别 location area identification (lai)位置取消规程 location cancellation procedure位置信息 location information位置信息恢复规程 location information retrieval procedure位置信息请求规程 location information requested procedure 温度补偿瓷介电容器temperature compensated ceramic capacitor温度循环试验 temperature cycling test文件图像（业务） telematics (services)文件用户电报 teletex纹波 ripple稳定度 stability稳压电源 regulated voltage supply稳压器 voltage regulator稳压器的基准电压 reference voltage of regulator。
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 提出基于均匀化方法的结构拓扑优化设计基本理论以来,均匀化方法应用到具有周期性结构的材料分析中,近几年该方法已经成为分析夹杂、纤维增强复合材料、混凝土材料等效模量,以及材料的细观结构拓扑优化常用的手段之一。
其基本思想是在组成拓扑结构的材料中引入微结构，优化过程中以微结构的几何尺寸作为设计变量，以微结构的消长实现其增删，并产生介于由中间尺寸微结构组成的复合材料，从而实现了结构拓扑优化模型与尺寸优化模型的统一。
文章就是通过均匀化的基础，结合拓扑结构优化的工程实际，以计算机模拟的方法将拓扑优化的一般过程呈现出来，有助于初涉拓扑优化的读者对拓扑优化有个基础的认识。
二、 拓扑优化问题描述为了简化问题的描述，文中假设设计域是简单的矩形形式，且在进行有限元离散的时候采用正方形单元对其进行离散。
这样不仅便于进行单元离散和单元编号，也利于对结构进行几何外形的描述。
一般说来，基于指数逼近法的拓扑优化最小化的问题可作如下描述：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)];%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%。
拓扑优化中的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)）。
这里，假设物质性能使恒定不变的同时每个元素是设计区域离散化，变量是元素的相对密度。
物质属性在相对物质密度增加到固体材料的物质属性的很多倍时被模板化。
这种方法曾一度引起争议因为它认为没有任何物理材料的属性特征能被幂律法则所描述。
然而，Bendsøe and Sigmund (1999)最近发表的文章证实只要在单一条件下能量能够满足，幂律法则在物理上就是可行的。
为了确认这个结论的存在性，幂律法则必须与周界约束、梯度约束或者过滤技术相结合。
这个拓扑优化的幂律法则方法已经应用于多重约束，多属性，多材料的问题中了。
然而，上面提及到的解决方法是基于数学规划法和连续设计变量法，一系列的文章都有涉及到解决拓扑优化的整数问题。
Beckers (1999)通过双重途径成功的解决了大规模服从最小化问题，但是其他方法大多基于遗传算法或者其他为了几个元素需要成千上万的函数求值半随机方法，当然这很可能是不切实际的。
除了上述提到的方法外，他们都能解决目标明确的问题，一些减少服从或者目标函数的启发式或直觉的方法也已经被提出来。
这些理论都被统称为进化设计理论。
除了很容易的理解和运用外，进化性分析方法主要的动机似乎是基于数学或连续变量法“涉及一些微积分应用和数学分析”，并且他们包含“一些复杂问题的数学理论”，反之进化性的方法能很好的应用强大的计算处理技术和自然进化过程中的直觉理念。
两件事可以反对他。
第一，曾经由于更多的比服从最小化理论复杂的目标被考虑进去。
进化论方法自身已经变得非常复杂。
第二，正如文中所说，以数学理论为基础的方法解决服从最小化问题很容易实现并且计算处理同样很有效率。
不仅如此，基于数学理论的方法很容易扩展解决像非共轭和多物理量的无服从问题以及多约束问题。
而用扩展进化方法来处理这些问题似乎更加不可行。
完整的matlab代码在附录中给出，文章剩余部分包括对优化问题的定义和讨论（第2部分），matlab实现的详细解释（第3部分），接着是扩展问题的讨论（第4部分）和最后总结（第5部分）。
2 拓扑优化问题有许多简化方法都是用来简化matlab的代码。
第一，假设设计区域是矩形且被平方有限元离散化。
这种情况下，元素的数目和结点就很容易表示（一列一列从左上角开始）并且结构的纵横比通过水平（nelx）和垂直（nely）方向元素比率来确定。
一个拓扑优化问题基于指数定律法，目标是实现最小化，可以如下表示（1）式中U 和F 分别表示整体变形和力的向量，K 表示整体刚度矩阵，e u 和e k 分别表示元素的位移矢量和刚度矩阵，x 是设计变量的向量，min x 是相对密度的最小向量（非零以避免奇点）N （=nelx * nely ）是设计区域离散化的元素数目，p 是惩罚因子（通常p=3），(x)V 和0V 是材料体积和设计区域体积，f (volfrac)为规定的容积率。
该拓扑优化问题（1）可以用几种不同的方法来解决，如Optimality Criteria （OC 算法），Sequential Linear Programming (SLP)，Method of Moving Asymptotes (MMA by Svanberg 1987)或者其他的方法。
为简单起见，我们这里用标准的OC 算法。
根据Bendsøe (1995)对设计变量应用启发式的调整法，公式表示为其中m （移动量）是一个正的移动界限，η (= 1/2)是数值阻尼系数，e B 在最优化条件中表示为其中λ是拉格朗日乘子，在bi-sectioning algorithm中获得。
目标函数的敏感度被表示为更多的详细信息将通过最优化准则理论的推导和实现，在文献(e.g. Bendsøe 1995).中提到。
为了确保对拓扑优化问题（1）解决方法的存在性，将会介绍一些作为结果的限制条件。
这里我们使用过滤技术，需要强调的是这个过滤技术并不能证明方法的可行性，但是通过作者大量应用证明过滤技术在实际中能产生独立性网格。
这个网格独立性滤波器是来修改元素敏度，如下:H表示为这个卷积算子（重量系数）f控制函数dist(e, f)定义为元素中心e到元素中心f的距离。
在过滤面积以外卷积H为0.卷积算子随着元素f的距离成线性减少。
除了原始的敏度（4），修算子f改的元素敏度（5）被用于OC算法（3）。
3 matlab算法实现Matlab代码(见附录),建立了作为一个标准拓扑优化的代码。
他的主程序为其中nelx和nely分别是横轴和纵轴方向上元素的个数，volfrac是容积率，penal是惩罚因子，rmin是过滤半径。
其他的变量如边界条件在matlab代码中定义且可根据实际需要进行调整。
在拓扑优化每一次循环迭代中，程序都会生成一幅当前密度分布图来。
图1显示了附录代码在给定以下输入条件下的结果密度分布默认的边界条件对应于“黑梁”的一半(图1)。
垂直方向的负载承受在左上角，水平方向结构左边和右边为对称边界。
Matlab代码中重要的细节将在下面各小节讨论。
3.1 主程序（1-36行）主程序的开始是材料设计区域的均匀分配（第4行）。
一些基本初始化后，主要的循环开始调用有限元分析子程序（12行）并返回位移矢量U。
因为固体材料的元素刚度矩阵对任何元素来说都是一样的，元素刚度矩阵子程序仅被调用一次（14行）。
接下来，一个循环遍历所有元素(第16-24行)来确定目标函数和敏感性(4)。
变量n1和n2表示在所有节点数中左上角和右上角的元素节点数目，并用于从整体位移矢量U中提取元素位移矢量U e。
灵敏度分析将伴随着调用网格独立性过滤器(第26行)和最优准则优化器(28行)。
当前的合规以及其它参数在30-33行中体现出来，并绘制出最后结果的密度分布图(第35行)。
主循环会在改变值（改变值在30行确定）小于1%是终止，否则以上步骤将重复执行。
3.2 基于优化程序的最优化准则（37-48行）更新后的设计变量是通过优化程序找到的（37-48行）。
由于材料体积(sum(sum(xnew)))是一个拉格朗日乘子的单调递减函数，满足体积约束的拉格朗日乘子数值可以通过bi-sectioning algorithm算法得到（40-48行）。
bi-sectioning algorithm的初始化时通过假设一个小于l1和大于l2结合的拉格朗日乘数（39行）。
限定拉格朗日乘数的间隔大小会一直减半直到它的大小小于收敛判别准则（40行）。
3.3 网格独立性过滤技术（49-64行）Matlab中第49-64行是为了实现条件（5）。
注意，并不是在设计域里所有的元素都要被搜索来寻找那些在处于过滤半径之内的元素，而是在那些被考虑元素周围，2倍过滤半径大小的区域。
通过选择过滤半径rmin小于调用子程序后的半径，被过滤后的敏度将会和初始状态一样，并使得过滤器不起作用。
3.4有限元分析代码（65-90行）有限元分析代码位于65-90行中。
注意解决者是利用matlab中的稀疏选项。
整体刚度矩阵是由所有元素的循环建立的（70-77行）。
正如主程序中，变量n1和n2指示在所有节点数目中左上部分和右边节点数目并插入到整体刚度矩阵合适的位置中。
正如前面所说，节点和元素都是从左到右纵向排布。
而且，每个节点都有2个自由度（水平和垂直的），因此指令F(2,1)=-1.第79行应用了一个在左上角的垂直单元力。
消除线性方程中的固定自由度来实现支承结构。
Matlab可以很好的处理，由这一行其中freedofs表示不受约束的自由度。
通常来说，定义一个固定的自由度要简单的多，其后会被自动找到，通过使用matlab算子setdiff因无约束自由度与固定自由度的不同来找到无约束自由度（82行）。
元素的刚度矩阵在86-99行中进行计算。
这个针对4个节点的8*8矩阵可通过使用人工智能软件进行处理分析。
杨氏模量E和泊松比可以在88和89行进行改变。
4 拓展附录给出的matlab代码解决如图1中材料分配的优化问题以达到最小化原则。
而算法中许多的拓展和改变都值得思考，其中一部分将在下面进行阐述。
4.1 其他边界条件为解决其他优化问题，改变边界条件和支承条件非常容易。
为了解决如图2所示的短悬臂梁问题，只需要要79和80作相应改变，根据这些改变，输入行变为4.2 多重负载问题扩展算法来处理多重负载问题也很容易。
实际上，仅需额外增加3行，并对其他4行做微小调整即可。
就两个负载工况而言，力和适量位移必须定义为两方向向量，这就意味着第69行改变为现在目标函数为二维之和因此第20-22行即被以下行所替代为解决如图3的二重负载问题，一个右上方的负载要添加到底79行，如下4.3 无源元件在某些情况下,一些元素可能需要采取最低密度值(例如一个管上的洞)。
