toph.zip

% a 99 line topology optimization code by Ole Sigmund,October 199 % toph(40,40,0.4,3.0,1.2) function toph(nelx,nely,volfrac,penal,rmin) % initialize x(1:nely,1:nelx)=volfrac; % x是设计变量 loop=0; % 存放迭代次数的变量 change=1.; % 每次迭代，x的改变值，用来判断何时收敛 % start ineration while change>0.01 %当两次连续目标函数迭代的差<=0.01时，且迭代次数小于50次，迭代结束 loop=loop+1; xold=x; %把前一次的设计变量付给xold % FE analysis [U]=FE(nelx,nely,x,penal); %有限元分析，得到位移矢量U % 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([n1;n2;n2+1;n1+1],1); %所示单元的自由度，左上，右上，右下，左下 c=c+(0.001+0.999*x(ely,elx)^penal)*Ue'*KE*Ue; % 计算目标函数的值（柔度） dc(ely,elx)=-0.999*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的改变量 abs绝对值，两个max指x行y列。 disp(['It.:' sprintf( '%4i',loop) ' Obj.:' sprintf(' %10.4f',c) ... %It 迭代，obj目标函数（柔度）， ' Vol.:' sprintf('%6.3f',sum(sum(x))/(nelx*nely)) ... %vol材料体积比 ' ch.:' sprintf('%6.3f',change)]) % plot densities 优化结果的图形显示 colormap(gray);imagesc(-x);axis equal;axis tight; axis off;pause(1e-6); % axis坐标轴控制函数 axis off 坐标系设为不可见，（1e-6）10的-6次方 end UU=reshape(U,nely,nelx); contourf(UU,8); colorbar;% imagesc 函数显示灰度图像 下面的代码是具有两个输入参数的 imagesc 函数显示一副灰度图像 imagesc(1,[0,1]); colormap(gray); imagesc 函数中的第二个参数确定灰度范围。 end % FE analysis 有限元求解子程序 function [U]=FE(nelx,nely,x,penal) %自定义函数，最后返回[U] [KE]=lk; % sparse 把一个全矩阵转化为一个稀疏矩阵，只存储每一个非零元素的三个值：元素值，元素的行号和列号 % 总体刚度矩阵的稀疏矩阵 % *2是因为x，y都有一个数 K=sparse((nelx+1)*(nely+1),(nelx+1)*(nely+1)); %k表示整体刚度矩阵 F=sparse((nely+1)*(nelx+1),1); %F表示力的向量 U=sparse((nely+1)*(nelx+1),1); %U表示整体变形 for elx=1:nelx for ely=1:nely %一列列的排序，y*x n1=(nely+1)*(elx-1)+ely; %左上 n2=(nely+1)*elx +ely; %右上 % 左上，右上，右下，左下 自由度 % 一个点有两个，所以要*2。第一个从1开始，所以*2之后要-1。 edof=[n1;n2;n2+1;n1+1]; K(edof,edof)=K(edof,edof)+(0.001+0.999*x(ely,elx)^penal)*KE; % 将单元刚度矩阵组装成总的刚度矩阵 end end % define loads and supports 定义荷载和支撑（半MBB梁） F(:,1)=0.01; %施加载荷 垂直单元力 %施加位移约束 fixeddofs =[nely/2+1-((nely/20)):2:nely/2+1+(nely/20)]; %固定节点 alldofs =[1:(nely+1)*(nelx+1)]; %固定所有节点 % 因无约束自由度与固定自由度的不同来找到无约束自由度 freedofs =setdiff(alldofs,fixeddofs); % solving 求解线性方程组，得到节点自由度的位移值 U(freedofs,:)=K(freedofs,freedofs)\F(freedofs,:); U(fixeddofs,:)=0; %受约束的节点自由度的位移值设为0 end % mesh-independency filter 独立网格过滤 敏度过滤技术子程序 function [dcn]=check(nelx,nely,rmin,x,dc) %dcn 清零，dcn用来保存更新的目标函数灵敏度 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) %floor（）函数时取整的意思 for l=max(j-floor(rmin),1):min(j+floor(rmin),nely) fac=rmin-sqrt((i-k)^2+(j-l)^2); %卷积算子（重量系数）Hf=过滤半径-控制函数（两元素中心距离） 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 end % Element stiffness matrix 单元刚度矩阵 function [KE]=lk KE=[2/3 -1/6 -1/3 -1/6 -1/6 2/3 -1/6 -1/3 -1/3 -1/6 2/3 -1/6 -1/6 -1/3 -1/6 2/3]; 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; else l2=lmid; end end end