91行散热（传热）结构拓扑优化matlab程序

% a 91 line topology optimization code by Ole Sigmund,October 1999 clear nelx=40; nely=40; volfrac=0.4; penal=3.; rmin=1.2; % initialize x(1:nely,1:nelx)=volfrac; loop=0; change=1; % start ineration while change>0.001 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([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)) ; 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 figure(1); colormap(gray);imagesc(-x);axis equal;axis tight; axis off;pause(1e-6); figure(2); U2=full(U); T=reshape(U2,nelx+1,nely+1); pcolor(T(1:(nelx+1),1:(nely+1)));xlabel('x');ylabel('y');title('Temperature distribution');colorbar; end % FE analysis function [U]=FE(nelx,nely,x,penal) [KE]=lk; K = sparse((nelx +1) *(nely+1), (nelx+1) *(nely+1)); F = sparse((nely +1) *(nelx+1),1); U = sparse((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 = [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 F(:,1)=0.01; fixeddofs = [nely/2+1-(nely/20):1:nely/2+1+(nely/20)]; % fixeddofs = [1,2,nely+2,2*nely+3]; alldofs=[1:(nely+1)*(nelx+1)]; freedofs =setdiff(alldofs,fixeddofs); % solving % U(freedofs,:)=K(freedofs,freedofs)\F(freedofs,:)+273; U(freedofs,:)=333+exp(-K(freedofs,freedofs)\F(freedofs,:)); U(fixeddofs,:)=273; end % 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 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