Topology Optimization 99

[MATLAB]拓扑优化99行代码

% a 99 line topology optimization code by Ole Sigmund,October 1999 
%example 60,20,0.5,3,1.5
function top99(nelx,nely,volfrac,penal,rmin)
%% initialize 初始化
x(1:nely,1:nelx)=volfrac; 
loop=0; 
change=1; 
%% start ineration 开始迭代
while change>0.01 
    loop=loop+1; 
    xold=x;
    % FE analysis 有限元分析
    [U]=FE(nelx,nely,x,penal); % 有限元模型分析，将计算得到的各个节点位移值保存在数组U中 
    % objective function and sensitivity analysis 目标函数和灵敏度分析
    [KE]=lk;  % 计算单元刚度阵 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); % 从数组 U 中根据节点*度编号 提取该单元的 位移向量 
            c=c+x(ely,elx)^penal*Ue'*KE*Ue; %  使用SIMP(Solid Isotropic Material with Penalization)材料插值模型
            dc(ely,elx)=-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)) %计算x的改变量
    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    
%% FE analysis 有限元求解子程序
function [U]=FE(nelx,nely,x,penal) 
[KE]=lk; 
K=sparse(2*(nelx+1)*(nely+1),2*(nelx+1)*(nely+1)); 
F=sparse(2*(nely+1)*(nelx+1),1);
U=sparse(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(edof,edof)=K(edof,edof)+x(ely,elx)^penal*KE; 
    end 
end 
%% define loads and supports 定义负载和支撑
F = sparse(2*(nely+1)*(nelx+1),1,-1,2*(nely+1)*(nelx+1),1);
fixeddofs = union([1:1:2*(nely+1)],[1]);
alldofs   =[1:2*(nely+1)*(nelx+1)]; 
freedofs  =setdiff(alldofs,fixeddofs); 
% solving 求解
U(freedofs,:)=K(freedofs,freedofs)\F(freedofs,:); 
U(fixeddofs,:)=0;

%% 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

%% Element stiffness matrix 单元刚度矩阵的子程序
function [KE]=lk 
E=1.;   %杨氏模量
nu=0.3; %泊松比 
k=[1/2-nu/6 1/8+nu/8 -1/4-nu/12 -1/8+3*nu/8 ... 
    -1/4+nu/12 -1/8-nu/8 nu/6   1/8-3*nu/8]; 
KE=E/(1-nu^2)*[ k(1)  k(2)  k(3)  k(4)  k(5)  k(6)  k(7)  k(8)   
                k(2)  k(1)  k(8)  k(7)  k(6)  k(5)  k(4)  k(3) 
                k(3)  k(8)  k(1)  k(6)  k(7)  k(4)  k(5)  k(2) 
                k(4)  k(7)  k(6)  k(1)  k(8)  k(3)  k(2)  k(5) 
                k(5)  k(6)  k(7)  k(8)  k(1)  k(2)  k(3)  k(4) 
                k(6)  k(5)  k(4)  k(3)  k(2)  k(1)  k(8)  k(7) 
                k(7)  k(4)  k(5)  k(2)  k(3)  k(8)  k(1)  k(6) 
                k(8)  k(3)  k(2)  k(5)  k(4)  k(7)  k(6)  k(1)];

%% 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))))); %这里是 OC 算法的核心所在
    if sum(sum(xnew))-volfrac*nelx*nely>0 
        l1=lmid; 
    else 
        l2=lmid; 
    end 
end