top88r_88行热传导程序_

%%%% AN 88 LINE TOPOLOGY OPTIMIZATION CODE %%%% function topr2(nelx,nely,volfrac,penal,rmin,ft) %ft指定应使用灵敏度过滤（ft=1）还是密度过滤（ft=2） %% MATERIAL PROPERTIES E0 = 1; %材料的刚度 Emin = 1e-3; %分配给空穴区域的非常小的刚度 %%PERPARE FINITE ELEMENT ANALYSIS 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]; %计算单元刚度矩阵 nodenrs = reshape(1:(1+nelx)*(1+nely),1+nely,1+nelx); %reshape：将A 的行列排列成m行n列，按照列取数据的 edofVec = reshape(nodenrs(1:end-1,1:end-1)+1,nelx*nely,1); %比如nodenrs有m行n列，nodenrs(1:end-1,1:end-1)表示取矩阵nodenrs的前m-1行前n-1列 edofMat = repmat(edofVec,1,4)+repmat([0 nely+[1 0] -1],nelx*nely,1); %B=remat(A,M,N)是以A的内容堆叠在（MxN）的矩阵B中 iK = reshape(kron(edofMat,ones(4,1))',16*nelx*nely,1); %kron(A,B):A每个元素都乘以矩阵B,然后堆叠。 jK = reshape(kron(edofMat,ones(1,4))',16*nelx*nely,1); % DEFINE LOADS AND SUPPORTS (HALF MBB-BEAM) F = sparse(1:(nely+1)*(nelx+1),1,0.01,(nely+1)*(nelx+1),1); U = zeros((nely+1)*(nelx+1),1); % fixeddofs = union([1:2:(nely+1)],[(nelx+1)*(nely+1)]); fixeddofs = [nely/2+1-((nely/20)):2:nely/2+1+(nely/20)]; alldofs=[1:(nely+1)*(nelx+1)]; freedofs=setdiff(alldofs,fixeddofs); %% PERPARE FILTER iH = ones(nelx*nely*(2*(ceil(rmin)-1)+1)^2,1); %w=ceil(z)函数将输入z中的元素取整，值w为不小于本身的最小整数。 jH = ones(size(iH)); sH = zeros(size(iH)); k=0; for i1 = 1:nelx for j1 = 1:nely e1 = (i1-1)*nely+j1; for i2 = max(i1-(ceil(rmin)-1),1):min(i1+(ceil(rmin)-1),nelx) for j2=max(j1-(ceil(rmin)-1),1):min(j1+(ceil(rmin)-1),nely) e2 =(i2-1)*nely+j2; k = k+1; iH(k) = e1; jH(k) = e2; sH(k) = max(0,rmin-sqrt((i1-i2)^2+(j1-j2)^2)); end end end end H = sparse(iH,jH,sH); Hs = sum(H,2); %% INITIULIZE ITERATION x = repmat(volfrac,nely,nelx); xPhys = x; %物理密度 loop = 0; change = 1; %% START ITERATION while change>0.01 loop = loop+1; %%FE-ANALYSIS sK = reshape(KE(:)*(Emin+xPhys(:)'.^penal*(E0-Emin)),16*nelx*nely,1); K = sparse(iK,jK,sK);K = (K+K')/2; %第二个语句确保刚度矩阵完全对称 U(freedofs)= K(freedofs,freedofs)\F(freedofs); %% OBJECTIVE FUNCTION AND SENSITIVITY ANALYSIS ce = reshape(sum((U(edofMat)*KE).*U(edofMat),2),nely,nelx); c = sum(sum((Emin+xPhys.^penal.*(E0-Emin)).*ce)); dc = -penal*(E0-Emin)*xPhys.^(penal-1).*ce; dv = ones(nely,nelx); %% FILTERING/MODIFICATION OF SENSITIVITIES if ft == 1 dc(:)=H*(x(:).*dc(:))./Hs./max(1e-3,x(:)); elseif ft == 2 dc(:) = H*(dc(:)./Hs); dv(:) = H*(dv(:)./Hs); end %% OPTIMALITY CRITERIA UPDATE OF DESIN VARIABLES AND PHYSICAL DENSITIES l1=0; l2=1e9; move=0.2; while(l2-l1)/(l1+l2)> 1e-3 % 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./dv/lmid))))); if ft == 1 xPhys = xnew; elseif ft == 2 xPhys(:)= (H*xnew(:))./Hs; end if sum(xPhys(:))>volfrac*nelx*nely, l1=lmid; else l2=lmid; end end change = max(abs(xnew(:)-x(:))); x=xnew; %% PRINT RESULT fprintf(' It.:%4i Obj.:%10.4f Vol.:%6.3f ch.:%6.3f\n',loop,c,... mean(xPhys(:)),change); %% PLOT DENSITIES colormap(gray); imagesc(1-xPhys); caxis([0 1]); axis equal; axis off; drawnow; end