拓扑优化99行代码详解注释版top99.m

本代码是作者为了方便拓扑优化初学者尽快入门而制作的，在原99行代码的基础上详细注释了每行代码的程序翻译，并与《A 99 line topology optimization code written in Matlab》O. Sigmund 原文公式进行对照标注。

%%%% A 99 LINE TOPOLOGY OPTIMIZATION CODE BY OLE SIGMUND OCTOBER 1999 %%% 
function top（nelxnelyvolfracpenalrmin）;%nelx和nely分别是水平和垂直方向上的元素数量volfrac容积率%penal惩罚因子一般为3.0，rmin过滤半径
nelx=60;nely=30;volfrac=0.5;penal=3;rmin=1.2;
% INITIALIZE 
x（1:nely1:nelx） = volfrac;%均匀分布材料
loop = 0; 
change = 1.; 
% START ITERATION 
while change > 0.01 
loop = loop + 1; 
xold = x; 
% FE-ANALYSIS 
[U]=FE（nelxnelyxpenal）;%有限元子程序 
% 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（elyelx）^penal*Ue‘*KE*Ue; %全局柔度最小，刚度最大，目标函数
     dc（elyelx） = -penal*x（elyelx）^（penal-1）* Ue‘*KE*Ue;%目标函数的灵敏度
   end 
end 
% FILTERING OF SENSITIVITIES 网格独立过滤器
[dc] = check（nelxnelyrminxdc）;
% DESIGN UPDATE BY THE OPTIMALITY CRITERIA METHOD 最优标准优化器
[x] = OC（nelxnelyxvolfracdc）; %%定义优化准则函数
% PRINT RESULTS 
% A=sum（sum（x））
change = max（max（abs（x-xold）））;%如果设计变量的变化（第30行中确定的变化）小于1％，则终止主循环，abs表示取绝对值
disp（[‘It.:‘ sprintf（‘%4i‘loop）   ‘Obj.:‘  sprintf（‘%10.4f‘c） ... 
‘ Vol.: ‘ sprintf（‘%6.3f‘sum（sum（x））/ （nelx*nely）） ... 
‘ ch.: ‘ sprintf（‘%6.3f‘change ）]） %sum（x）表示[x]的列总和sum（sum（x））/ （nelx*nely）表示体积分数，It是迭代次数，c是全局刚度，ch是设计变量的变化
% PLOT DENSITIES最后结果的密度分布图
colormap（gray）; imagesc（-x）; axis equal; axis tight; axis off;pause（1e-6）; 
end 
%%%%%%%%%% OPTIMALITY CRITERIA UPDATE %%%%%%%%% 基于优化程序的最优化准则oc法
function [xnew]=OC（nelxnelyxvolfracdc） 
l1 = 0; l2 = 100000; move = 0.2; 
while （l2-l1 > 1e-4）% bi-sectioning算法求出满足体积约束的拉格朗日乘数数值（41-49）小于l1，大于l2的这个范围，限定拉格朗日乘数的间隔大小会一直减半直到他的大小小于收敛判别准则
lmid = 0.5*（l2+l1）; 
% A=x+move
xnew = max（0.001  max（x-move  min（1. min（x+movex.*sqrt（-dc./lmid）））））; %开根号 1/2次方代表数值阻尼系数η x.*sqrt（-dc./lmid）为论文中的XeBeη
if sum（sum（xnew）） - volfrac*nelx*nely > 0 
l1 = lmid; 
else
l2 = lmid; 
end 
end 
%%%%%%%%%% MESH-INDEPENDENCY FILTER %%%%%%%%%%% 网格独立性过滤技术
function [dcn]=check（nelxnelyrminxdc） 
dcn=zeros（nelynelx）;
for i = 1:nelx 
    for j = 1:nely 
     sum=