99行代码

A 99 line topology optimization code written in Matlab(2001)

概述

最基础的拓扑优化代码

FE代码

有限元分析：计算单元刚度矩阵，组装全局刚度矩阵，求解单元位移

function [U]=FE(nelx,nely,x,alfa,xfar)
[KE] = lk; %单元刚度矩阵
K = sparse(2*(nelx+1)*(nely+1), 2*(nelx+1)*(nely+1));
    %总体刚度矩阵的稀疏矩阵 注意节点数和单元数的关系
F = sparse(2*(nely+1)*(nelx+1),1); %载荷向量 有这么多节点(nely+1)*(nelx+1)
U = zeros(2*(nely+1)*(nelx+1),1); %位移向量 有这么多节点(nely+1)*(nelx+1)
% 总刚组装
for elx = 1:nelx   % 第elx行的单元
  for ely = 1:nely  %第ely列的单元
    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;  
       K(edof,edof) = K(edof,edof) +  (1+tanh(alfa*( x(ely,elx)-xfar)))*KE;  %将单元刚度矩阵组装成总的刚度矩阵（实际是取出edof构成的8×8矩阵，然后填入）
%             K(edof,edof) = K(edof,edof) + (1/(1+exp(-alfa*(x(ely,elx)-0.5))))*KE;  %将单元刚度矩阵组装成总的刚度矩阵
  end
end
% DEFINE LOADS AND SUPPORTS (HALF MBB-BEAM)
F(2,1) = -1; % 其他元素都是0

fixeddofs   = union([1:2:2*(nely+1)],[2*(nelx+1)*(nely+1)]); % 固定左边一列+右下角顶点x坐标
alldofs     = [1:2*(nely+1)*(nelx+1)];                       %所有结点 
freedofs    = setdiff(alldofs,fixeddofs);  % setdiff返回A中有，B中没有的值自由节点
% SOLVING
U(freedofs,:) = K(freedofs,freedofs) \ F(freedofs,:); %求解   
U(fixeddofs,:)= 0;

灵敏度过滤

去除棋盘格，先设置一个过滤半径，对过滤半径内的单元与当前单元距离计算权重（与距离成反比的权重函数），对密度乘以权重然后累加，最后除以加权总和，更新灵敏度

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

OC法

正常情况下，对于有解析解的形式同样是构造拉格朗日函数

\[ L(x,\lambda)=c(x)+\lambda(\sum x_i-V_0) \]

对其求极值并通过乘法形式更新，由于现在没有解析解，故采用二分查找 \(\lambda\)

\[ x_i^{new}=x_i\cdot\sqrt{-\frac{\partial c/\partial x_i}{\lambda}} \]

function [xnew]=OC(nelx,nely,x,volfrac,dc)  
L1 = 0; L2 = 100000; move = 0.2;
while (L2-L1 > 1e-4) %l 是小写的L 不是12 11  

  Lmid = 0.5*(L2+L1);
  xnew = max(0.001,max(x-move,min(1.,min(x+move,x.*sqrt(-dc./Lmid)))));
  if sum(sum(xnew)) - volfrac*nelx*nely > 0;
    L1 = Lmid;
  else
    L2 = Lmid;
  end
end

论文

论文地址