99_BESO matalb

拓扑优化程序99行MATAB代码，针对目标简支梁的半个对称模型，联合BESO算法，对柔顺度进行拓扑优化分析。

%固定两边中间的点，在最上方中点施加压力
function top（nelxnelyvolfracpenalrmin）;
nelx=80;
nely=50;
volfrac=0.5;
penal=3;
rmin=3;
% INITIALIZE
x（1:nely1:nelx） = volfrac; 
loop = 0; 
change = 1.;
% START ITERATION
while change > 0.01  
  loop = loop + 1;
  xold = x;
% FE-ANALYSIS
  [U]=FE（nelxnelyxpenal）;  
% 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（elyelx）^penal*Ue‘*KE*Ue;
      dc（elyelx） = -penal*x（elyelx）^（penal-1）*Ue‘*KE*Ue;
    end
  end
% FILTERING OF SENSITIVITIES
  [dc]   = check（nelxnelyrminxdc）;    
% DESIGN UPDATE BY THE OPTIMALITY CRITERIA METHOD
  [x]    = OC（nelxnelyxvolfracdc）; 
% PRINT RESULTS
  change = max（max（abs（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 
%%%%%%%%%% OPTIMALITY CRITERIA UPDATE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [xnew]=OC（nelxnelyxvolfracdc）  
l1 = 0; l2 = 100000; move = 0.2;
while （l2-l1 > 1e-4）
  lmid = 0.5*（l2+l1）;
  xnew = max（0.001max（x-movemin（1.min（x+movex.*sqrt（-dc./lmid）））））;
  if sum（sum（xnew）） - volfrac*nelx*nely > 0;
    l1 = lmid;
  else
    l2 = lmid;
  end
end
%%%%%%%%%% MESH-INDEPENDENCY FILTER %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [dcn]=check（nelxnelyrminxdc）
dcn=zeros（nelynelx）;
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（0fac）;
        dcn（ji） = dcn（ji） + max（0fac）*x（lk）*dc（lk）;
      end
    end
    dcn（ji） = dcn（ji）/（x（ji）*sum）;
  end
end
%%%%%%%%%% FE-ANALYSIS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [U]=FE（nelxnelyxpenal）
[KE] = lk; 
K = sparse（2*（nelx+1）*（nely+1） 2*（nelx+1）*（nely+1））;
F = sparse（2*（nely+1）*（nelx+1）1）; U = zeros（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（edofedof） = K（edofedof） + x（elyelx）^penal*KE;
  end
end
% DEFINE LOADS AND SUPPORTS （HALF MBB-BEAM）
%固定两边中间的点，在最上方中点施加压力
F（2*（nelx/2+1）*（nely+1）-2*nely1） = 1;
fixeddofs   = [2*（nely/2+1）2*nelx*（nely+1）+2*（nely/2+1）];
alldofs     = [1:2*（nely+1）*（nelx+1）];
freedofs    = setdiff（alldofsfixeddofs）;
U（freedofs:） = K（freedofsfreedofs） \ F（freedofs:）;      
U（fixeddofs:）= 0;
% SOLVING
%%U（freedofs :）= K（freedofsfreedofs） \ F（freedofs:??;      
%U（fixeddofs:??= 

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       3865  2015-10-15 11:21  top99____1.m

----------- ---------  ---------- -----  ----

                 3865                    1