99行代码的注释及拓扑优化入门详解

Sigmund在2001年在Structural and Multidisciplinary Optimization 上发表一篇名为 “A 99 line topology optimization code written in Matlab”论文。该论文后附带了一个Matlab拓扑优化程序。这个只有99行代码程序基于Matlab环境构建了一个完整的拓扑优化流程，其中包含：前处理(构建有限元仿真模型), 有限元模型分析计算,拓扑优化迭代和后处理(分析结果显示)。Sigmund在论文从拓扑优化角度对这个程序做了详细解释。本文仅从程序设计角度解析这段代码。

Sigmund的99行Matlab拓扑优化程序使用模块化方法设计，主要包含以下几个模块： 
- 程序主流程 
- 有限元模型求解模块 
- Filter模块 
- 单元刚度阵计算模块：计算平面四边形单元的刚度矩阵 
- 优化模块： 使用优化准则法更新设计变量 
- 目标函数计算和灵敏度分析模块

（若还有疑问，详情参见 https://blog.csdn.net/cocoonyang/article/details/80494678及北京理工大学，王路，拓扑优化优化学习报告）

99行拓扑优化代码

Sigmund的99行Matlab拓扑优化程序如下所示

% a 99 line topology optimization code by Ole Sigmund,October 1999 
clear 
nelx=60; 
nely=40; 
volfrac=0.5; 
penal=3.; rmin=1.5; % initialize x(1:nely,1:nelx)=volfrac; loop=0; change=1; % start ineration while change>0.01 loop=loop+1; xold=x; % FE analysis [U]=FE(nelx,nely,x,penal); % objective function and sensitivity analysis [KE]=lk;; c=0.; for ely=1:nely for elx=1:nelx n1=(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; end end % filtering of sensitivities [dc]=check(nelx,nely,rmin,x,dc); % design update by the optimality criteria method [x]=oc(nelx,nely,x,volfrac,dc); % print result change=max(max(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 densities colormap(gray);imagesc(-x);axis equal;axis tight; axis off;pause(1e-6); end % 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=sparse(2*(nely+1)*(nelx+1),1); for elx=1:nelx for ely=1:nely n1=(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; end end % define loads and supports ip=(nelx+1)*(nely+1); F(2*ip,1)=-1; fixeddofs =[1:2*(nely+1)]; alldofs =[1:2*(nely+1)*(nelx+1)]; freedofs =setdiff(alldofs,fixeddofs); % solving U(freedofs,:)=K(freedofs,freedofs)\F(freedofs,:); U(fixeddofs,:)=0; % mesh-independency 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 dcn(j,i)=dcn(j,i)/(x(j,i)*sum); end end % Element stiffness matrix function [KE]=lk E=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(

转载于:https://www.cnblogs.com/hyb221512/p/9212924.html