matlab结构拓扑优化问题求教

Dwayne.s

大家好，小弟做结构拓扑优化的，小弟接触matlab编程不久，目前在看徐斌老师的《MTALAB》有限元结构动力学分析与工程应用这本书，在看第八章的例8.4：连续体拓扑优化的时候遇到了麻烦，由于查找文献的局限性，未能找到这方面的文献，所以对以下三部分的程序不解。望懂的大神能指导下小弟么？小弟定当感激不尽第一块是function [xnew]=OC(nelx,nely,x,h,U5,q,G,freq,d,g,freq1,p,a,b,dU,G1,G2)
%--------------------------------------------------------------------------------------
%  OPTIMALITY CRITERIA UPDATE
%-------------------------------------------------------------------------------------
step1=1;
for ely = 1:nely
    for elx = 1:nelx
     xnew1(ely,elx)=real(x(ely,elx)*U5(ely,elx))^(1/(1+2*h));     %  update design variables based displacement sensitivity
     xnew2(ely,elx)=real((((freq1/(freq(1)))^d)*(((G1/G2)*(G(ely,elx)/((1/dU(ely,elx)))))^g))*x(ely,elx));  %  update design variables based frequency sensitivity
     TT1((elx-1)*nely+ely)=xnew1(ely,elx);
     TT2((elx-1)*nely+ely)=xnew2(ely,elx);
   end
end


for ely = 1:nely
    for elx = 1:nelx
     xnew2(ely,elx)=xnew2(ely,elx)*(median(TT1)/median(TT2));  %  tension of different update design variables based different sensitivities  
     xnew(ely,elx)=min(max(xnew1(ely,elx),xnew2(ely,elx)),(xnew1(ely,elx)+xnew2(ely,elx))/2);  %  tradeoff of update design variables
     T((elx-1)*nely+ely)=xnew(ely,elx);
   end
end

%  deletion of certain percentage of elements using Evolutionary Structural Optimization(ESO)
Mt=sort(T);
     for i=1:(nely*nelx)
       if((0.00005)<Mt(i)&Mt(i)<(0.009))
        Maa(step1)=Mt(i);
        step1=step1+1;
       end
     end
TT=2*floor(step1*0.005)
for ely = 1:nely
    for elx = 1:nelx
     if((xnew(ely,elx)>0.000009)&(xnew(ely,elx)<0.009)&TT~=0)
        if(xnew(ely,elx)<=Maa(TT))
         xnew(ely,elx)=0.000009;
        end
     end
     TTT((elx-1)*nely+ely)=xnew(ely,elx);   
    end
end        

%  division of design variable based average value
for ely = 1:nely
    for elx = 1:nelx
     w(ely,elx)=q*sign(xnew(ely,elx)-mean(TTT))*xnew(ely,elx);
     xnew(ely,elx)=xnew(ely,elx)+w(ely,elx);
     xnew(ely,elx)=min(max(xnew(ely,elx),0.000009),0.009);
     if((abs(x(ely,elx)-0.000009)<=1e-5)|(abs(x(ely,elx)-0.009)<=1e-5))
     xnew(ely,elx)=x(ely,elx);
     end
    end
end
第二块是：function [G,G1,G2]=refdistrib(nelx,nely,G,p,a,b,dU)
%-------------------------------------------------------------------------
% frequency sensitivity and sensitivity redistrubution
%-------------------------------------------------------------------------
G1=0;
G2=0;
for ely = 1:nely+1
    for elx = 1:nelx+1
        sum2=0;
        MM2=0;
        aa2=[elx ely;elx-1 ely-1;elx-1 ely;elx ely-1];
        for i=1:4
         if(aa2(i,:)>0&aa2(i,1)<=nelx&aa2(i,2)<=nely)
          sum2=sum2+G(aa2(i,2),aa2(i,1));
           MM2=MM2+2;
         end
        end
       G3(ely,elx)=sum2/MM2;
    end
end
for ely = 1:nely
    for elx = 1:nelx
     sum21=0;
     MM21=0;
     aa21=[elx ely;elx+1 ely+1;elx ely+1;elx+1 ely];
      for i=1:4
       if(aa21(i,:)>0&aa21(i,1)<=nelx+1&aa21(i,2)<=nely+1)
       sum21=sum21+G3(aa21(i,2),aa21(i,1));
       MM21=MM21+1;
       end
      end
     G(ely,elx)=sum21/MM21;
    end
end
for ely = 1:nely
    for elx = 1:nelx
    G1=G1+(1/dU(ely,elx))*G(ely,elx);
    G2=G2+G(ely,elx)^2;
    end
end
第三块是：function [U5]=reddistrib(nelx,nely,x,dU,lamda,p,a,b,h)
%------------------------------------------------------------------------
% displacement sensitivity and sensitivity redistrubution
%---------------------------------------------------------------------------
for ely = 1:nely
    for elx = 1:nelx
      U4(ely,elx)=((lamda/(4*p*a*b))*dU(ely,elx)*x(ely,elx))^h;
  end
end
for ely = 1:nely+1
    for elx = 1:nelx+1
        sum1=0;
        MM1=0;
        aa1=[elx ely;elx-1 ely-1;elx-1 ely;elx ely-1];
        for i=1:4
         if(aa1(i,:)>0&aa1(i,1)<=nelx&aa1(i,2)<=nely)
          sum1=sum1+U4(aa1(i,2),aa1(i,1));
           MM1=MM1+1;
         end
        end
       U6(ely,elx)=sum1/MM1;
    end
end
for ely = 1:nely
    for elx = 1:nelx
     sum11=0;
     MM11=0;
     aa11=[elx ely;elx+1 ely+1;elx ely+1;elx+1 ely];
      for i=1:4
       if(aa11(i,:)>0&aa11(i,1)<=nelx+1&aa11(i,2)<=nely+1)
       sum11=sum11+U6(aa11(i,2),aa11(i,1));
       MM11=MM11+1;
       end
      end
     U5(ely,elx)=sum11/MM11;
    end
end
望大神们能指导下小弟！谢谢！