[MATLAB]拓扑优化99行代码
% a 99 line topology optimization code by Ole Sigmund,October 1999
%example 60,20,0.5,3,1.5
function top99(nelx,nely,volfrac,penal,rmin)
%% 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); % 有限元模型分析,将计算得到的各个节点位移值保存在数组U中
% objective function and sensitivity analysis 目标函数和灵敏度分析
[KE]=lk; % 计算单元刚度阵 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); % 从数组 U 中根据节点自由度编号 提取该单元的 位移向量
c=c+x(ely,elx)^penal*Ue'*KE*Ue; % 使用SIMP(Solid Isotropic Material with Penalization)材料插值模型
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)) %计算x的改变量
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)