Matlab top99 拓扑优化中经典算法.zip

function top99(nelx,nely,volfrac,penal,rmin) % 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(7) k(8) k(1) k(2) k(3) k(4) k(6) k(5) k(4) k(3) k(2) k(1) k(8) k(7) k(7) k(4) k(5) k(2) k(3) k(8) k(1) k(6) k(8) k(3) k(2) k(5) k(4) k(7) k(6) k(1)]; % 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