topmma.zip

% a 99 line topology optimization code by Ole Sigmund,October 199 %整块板均匀加热，左边界给定温度T=0，其余边界为绝热边界 function topmma(nelx,nely,volfrac,penal,rmin) % initialize x(1:nely,1:nelx)=volfrac; % x是设计变量 loop=0; % 存放迭代次数的变量 change=1.; % 每次迭代，x的改变值，用来判断何时收敛 % start ineration xval=reshape(x,nely*nelx,1); xold1=xval; xold2=xold1; low=[]; upp=[]; while change>0.04 %当两次连续目标函数迭代的差<=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; % 单位刚度矩阵 obj=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); %所示单元的自由度，左上，右上，右下，左下 obj=obj+(0.001+0.999*penal*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 最佳优化方式设计更新 % DEFINE INITIAL CONDITIONS m=1; n=nelx*nely; iter=loop; xmin=0.001*ones(n,1); xmax=ones(n,1); f0val=obj; df0dx=reshape(dc,n,1); df0dx2=0*df0dx; % 定义约束条件 fval=sum(sum(x))/(volfrac*n)-1; dfdx=ones(1,n)/(volfrac*n); dfdx2=0*dfdx; a0=1; a=0*ones(m,1); c=10000*ones(m,1); d=0*ones(m,1); % RUN MMA [xmma,ymma,zmma,lam,xsi,eta,mu,zet,s,low,upp] = ... mmasub(m,n,iter,xval,xmin,xmax,xold1,xold2, ... f0val,df0dx,df0dx2,fval,dfdx,dfdx2,low,upp,a0,a,c,d); % END UPDATE BY MMA METHOD xold2=xold1; xold1=xval; xval=xmma; x=reshape(xval,nely,nelx); % PRINT RESULTS change=max(max(abs(x-xold))); disp(['It.:' sprintf('%4i',loop) 'Obj.:' sprintf('%10.4f',obj) '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 % axis坐标轴控制函数 axis off 坐标系设为不可见，（1e-6）10的-6次方 end % imagesc 函数显示灰度图像 下面的代码是具有两个输入参数的 imagesc 函数显示一副灰度图像 imagesc(1,[0,1]); colormap(gray); imagesc 函数中的第二个参数确定灰度范围。 % 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 % This is the file mmasub.m % function [xmma,ymma,zmma,lam,xsi,eta,mu,zet,s,low,upp] = ... mmasub(m,n,iter,xval,xmin,xmax,xold1,xold2, ... f0val,df0dx,df0dx2,fval,dfdx,dfdx2,low,upp,a0,a,c,d); % % Written in May 1999 by % Krister Svanberg <krille@math.kth.se> % Department of Mathematics % SE-10044 Stockholm, Sweden. % % Modified ("spdiags" instead of "diag") April 2002 % % % This function mmasub performs one MMA-iteration, aimed at % solving the nonlinear programming problem:（这个函数mmasub执行一个mma迭代，目标是求解非线性规划问题:） % % Minimize f_0(x) + a_0*z + sum( c_i*y_i + 0.5*d_i*(y_i)^2 ) % subject to f_i(x) - a_i*z - y_i <= 0, i = 1,...,m % xmin_j <= x_j <= xmax_j, j = 1,...,n % z >= 0, y_i >= 0, i = 1,...,m %*** INPUT: % % m = The number of general constraints.（一般约束的数量） % n = The number of variables x_j.（变量的个数xj） % iter = Current iteration number ( =1 the first time mmasub is called).（当前迭代数(第一次调用mmasub时=1)。） % xval = Column vector with the current values of the variables x_j. （ xval =带有变量x_j的当前值的列向量。） % xmin = Column vector with the lower bounds for the variables x_j. （xmin =具有变量x_j下界的列向量。） % xmax = Column vector with the upper bounds for the variables x_j. （ xmax =变量x_j的上界列向量。） % xold1 = xval, one iteration ago (provided that iter>1). （ xold1 = xval，一次迭代前(假设iter>1)。） % xold2 = xval, two iterations ago (provided that iter>2). （xold2 = xval，两次迭代前(假设iter>2)。） % f0val = The value of the objective function f_0 at xval. （f0val =目标函数在xval处的值f_0。） % df0dx = Column vector with the derivatives of the objective function （目标函数导数的列向量f_0对于变量x_j，在xval处计算。） % f_0 with respect to the variables x_j, calculated at xval. % df0dx2 = Column vector with the non-mixed second derivatives of the （非混合二阶导数的列向量，目标函数f_0对于变量x_j按xval计算。df0dx2(j) =二阶导数 % objective function f_0 with respect to the variables x_j,f_0相对于x_j的%(2次)。重要提示:如果有二阶导数，）让df0dx2 = 0*df0dx。） % calculated at xval. df0dx2(j) = the second derivative % of f_0 with respect to x_j (twice). % Important note: If second derivatives are not available, % simply let df0dx2 = 0*df0dx. % fval = Column vector with the values of the constraint functions f_i, （约束函数f_i的值的列向量，按xval计算。） % calculated at xval. % dfdx = (m x n)-matrix