88行代码

Efficient topology optimization in MATLAB using 88 lines of code

概述

在99行代码的优化效率、增加密度过滤

对有限元组装、计算柔度、过滤操作等部分进行向量化，去除 for 循环
增加灵敏度过滤和密度过滤两种方式

优化下面问题

\[ \min_{x} \quad c(x) = \mathbf{U}^\mathrm{T} \mathbf{K}(x) \mathbf{U} = \sum_{e=1}^{N} E_e(x_e) \cdot \mathbf{u}_e^\mathrm{T} \mathbf{k}_0 \mathbf{u}_e \]

\[ \text{s.t.} \quad \begin{cases} \sum x_e \leq \text{volfrac} \cdot N \\ 0 \leq x_e \leq 1 \\ \mathbf{K}(x)\mathbf{U} = \mathbf{F} \end{cases} \]

SIMP方法：

\[ E(x_e) = E_{\min} + x_e^p (E_0 - E_{\min}) \]

其中 $ p $ 是惩罚因子（通常为 3）
引入 $ E_{\min} $ 避免刚度矩阵奇异性

装配全局刚度矩阵K

通过向量操作加速全局刚度矩阵K的装配

A11 = [12  3 -6 -3;  3 12  3  0; -6  3 12 -3; -3  0 -3 12];
A12 = [-6 -3  0  3; -3 -6 -3 -6;  0 -3 -6  3;  3 -6  3 -6];
B11 = [-4  3 -2  9;  3 -4 -9  4; -2 -9 -4 -3;  9  4 -3 -4];
B12 = [ 2 -3  4 -9; -3  2  9 -2;  4  9  2  3; -9 -2  3  2];
KE = 1/(1-nu^2)/24*([A11 A12;A12' A11]+nu*[B11 B12;B12' B11]);%单元刚度阵
nodenrs = reshape(1:(1+nelx)*(1+nely),1+nely,1+nelx);                           % 单元编号
edofVec = reshape(2*nodenrs(1:end-1,1:end-1)+1,nelx*nely,1);                    % 每个单元左下角自由度索引
edofMat = repmat(edofVec,1,8)+repmat([0 1 2*nely+[2 3 0 1] -2 -1],nelx*nely,1); % 每个单元8个自由度索引
iK = reshape(kron(edofMat,ones(8,1))',64*nelx*nely,1);                          % 全局刚度矩阵行索引
jK = reshape(kron(edofMat,ones(1,8))',64*nelx*nely,1);                          % 全局刚度矩阵列索引

灵敏度过滤预处理

预处理一个邻接矩阵H，表示灵敏度过滤中每个单元受领域内其他单元影响，权重是距离的反比例函数

iH = ones(nelx*nely*(2*(ceil(rmin)-1)+1)^2,1);
jH = ones(size(iH));
sH = zeros(size(iH));
k = 0;
for i1 = 1:nelx
    for j1 = 1:nely
    e1 = (i1-1)*nely+j1;
    for i2 = max(i1-(ceil(rmin)-1),1):min(i1+(ceil(rmin)-1),nelx)
        for j2 = max(j1-(ceil(rmin)-1),1):min(j1+(ceil(rmin)-1),nely)
        e2 = (i2-1)*nely+j2;
        k = k+1;
        iH(k) = e1;   % e1单元受e2单元影响
        jH(k) = e2;   % e1单元受e2单元影响
        sH(k) = max(0,rmin-sqrt((i1-i2)^2+(j1-j2)^2)); % 存储权重
        end
    end
    end
end
H = sparse(iH,jH,sH);  % 整体处理每个单元与周围邻居的距离加权影响。n * n
Hs = sum(H,2);         % 按行求和，求每个单元受领域内其他单元的影响，后面用于归一化 n * 1

密度过滤以及灵敏度过滤

密度型过滤

\[ \tilde{dc}_i = \frac{ \displaystyle\sum_{j \in N(i)} H_{ij} \cdot x_j \cdot dc_j }{ x_i \cdot \displaystyle\sum_{j \in N(i)} H_{ij} } \]

$\tilde{dc}_i$：单元 $i$ 的滤波后灵敏度
$dc_j $：邻居单元 $j$ 的原始灵敏度
$x_j$：邻居单元的密度
$H_{ij}$：距离加权系数
$N(i)$：单元 $i$ 的邻居集合

灵敏度型过滤

\[ \tilde{dc}_i = \frac{ \displaystyle\sum_{j \in N(i)} H_{ij} \cdot dc_j }{ \displaystyle\sum_{j \in N(i)} H_{ij} } \]

$dc_i$：单元 $i$ 的滤波后灵敏度
$dc_j$：邻居单元 $j$ 的原始灵敏度
$H_{ij}$：距离加权系数
$N(i)$：单元 $i$ 的邻居集合