利用双向渐进结构优化法（BESO ）优化连续体的结构拓扑

利用双向渐进结构优化法（BESO ）优化连续体的结构拓扑

%%%% A SOFT-KILL BESO CODE BY X. HUANG and Y.M. XIE %%%%
%%%% A further review of ESO type methods for topology optimization %%%%%
%%%% SBESO 方法 体积约束下的最小柔度问题%%%%%%%%%%%%
function sbeso（nelxnelyvolfracerrmin）
nelx=80;
nely=40;
volfrac=0.4;
penal = 3;  
rmin=2;
er=0.02;
%INITIALIZE
x（1:nely1:nelx）=1.;vol=1.;i=0;change=1.;penal=3.;
%START iTH ITERATION
while change>0.001
    i=i+1; 
    vol= max（vol*（1-er）volfrac）;
    if i>1;olddc=dc;
    end
 %FE-ANALYSIS
 [U]=FE（nelxnelyxpenal）;
 %objectIVE FUNCTION AND SENSITIVITY ANALYSIS
 [KE]=lk;
 c（i）=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（i）=c（i）+0.5*x（elyelx）^penal*Ue‘*KE*Ue;
         dc（elyelx）=0.5*x（elyelx）^（penal-1）*Ue‘*KE*Ue;
     end
 end
 %FILTERING OF SENSITIVITIES
 [dc]=check（nelxnelyrminxdc）;
 %STABILIZATION OF EVOLUTIONARY PROCESS
 if i>1;
     dc=（dc+olddc）/2.;
 end
 %BESO DESIGN UPDATE
 [x]=ADDDEL（nelxnelyvoldcx）;
 %PRINT RESULTS
 if i>10;
     change=abs（sum（c（i-9:i-5））-sum（c（i-4:i）））/sum（c（i-4:i））;
 end
 disp（[‘ It.: ‘ sprintf（‘%4i‘i） ‘ Obj.: ‘ sprintf（‘%10.4f‘c（i）） ...
       ‘ 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

%%%%%%%%%% OPTIMIZATION CRITERA UPDATE %%%%%%%%%%%%%%%%%%%%%%%%%%
function [x]=ADDDEL（nelxnelyvolfradcx）;
l1=min（min（dc））;l2=max（max（dc））;
while （（l2-l1）/l2>1.0e-5）
    th=（l1+l2）/2.0;
    x=max（0.001sign（dc-th））;
    if sum（sum（x））-volfra*（nelx*nely）>0;
        l1=th;
    else 
        l2=th;
    end
end

%%%%%%%%%% MESH-INDEPENDENCY FILTER %%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [dcf]=check（nelxnelyrminxdc）
dcf=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-1）^2）;
                sum=sum+max（0fac）;
                dcf（ji）=dcf（ji）+max（0fac）*dc（lk）;
            end
        end
        dcf（ji）=dcf（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 （Cantilever）
ny=ceil（nely/2）;
        F（2*（（nely+1）*nelx+（ny+1））1） = -1000;
fixeddofs=[1:2*（nely+1）];
alldofs  =[1:2*（nely+1）*（nelx+1）];
freedofs =setdiff（alldofsfixeddofs）;
% SOLVIN

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     目录           0  2014-03-27 07:35  基于软杀死的BESO拓扑优化\
     文件        3837  2014-03-26 20:31  基于软杀死的BESO拓扑优化\sbeso.m