经典 99 行 MATLAB 拓扑优化程序（可运行代码&&超详细注释）

原创已于 2025-11-12 10:34:32 修改 · 474 阅读

6 ·

CC 4.0 BY-SA版权

文章标签：

#算法 #matlab #拓扑学 #性能优化

于 2025-09-11 13:14:06 首次发布

这段代码是做什么的？

SIMP（Solid Isotropic Material with Penalization/固体各向同性材料惩罚模型）方法 求解 MBB 梁（半长） 的最优材料分布（即“拓扑优化”），目标是在给定体积分数下 最小化结构柔度（最大化刚度）。
说人话：“给定 60×30 网格、50 % 体积分数，自动找出最硬、最省材料的梁形状。”

什么是拓扑优化？

是一种通过数学方法对材料分布进行优化的设计技术，旨在在给定的设计空间内，通过调整材料布局来满足特定的性能目标（如刚度最大化、重量最轻等）。它广泛应用于结构设计、机械工程、航空航天、汽车工业等领域，能够显著提升产品的性能并减少材料浪费。
核心原理：

核心原理

拓扑优化的核心是通过有限元分析（FEA）和优化算法（如SIMP法、水平集方法等）迭代计算，逐步去除对结构性能贡献较小的材料，同时保留关键受力路径。最终生成轻量化且高强度的结构设计。

代码前面要输入参数才可以运行

单元位移矩阵

详情看这个up的视频https://www.bilibili.com/video/BV1h7421N7K3/?spm_id_from=333.788.videopod.sections&vd_source=02ffd836310dcb034a8bb8e545cdce63

可运行代码


function top99(nelx,nely,volfrac,penal,rmin)
% 输入
nelx=60 ;        %水平方向单元数（离散单元数）  矩阵列数
nely=20 ;        %竖直方向单元数               矩阵行数
volfrac=0.3;     %允许材料体积分数  将结构减重到原体积的百分之多少
penal=3.0 ;      %惩罚系数
rmin=1.5  ;      %过滤半径（网格无关性），防止出现棋盘格现象

% INITIALIZE   初始化  下面是将单元进行离散化
x(1:nely,1:nelx) = volfrac;      % 初始密度场 创建一个nely行nelx列的矩阵，用来存储个单元的伪密度
loop = 0;                        % 迭代计数
change = 1.;                     % 两次迭代最大变化量  要是浮点数后面要带点
                                 % 存储每次迭代后目标函数的改变值，用来判断是否收敛

% START ITERATION  iteration  开始迭代
while change > 0.01
    loop = loop + 1;
    xold = x;    %x用来存储每次迭代时的结果

% FE-ANALYSIS  
    [U] = FE(nelx,nely,x,penal);  % 调用有限元分析  （65——85）

    % OBJECTIVE FUNCTION AND SENSITIVITY ANALYSIS  目标函数和敏度分析
    [KE] = lk;   %调用单元刚度矩阵子程序 （86——99）
    c = 0.;      %c是目标函数，即最小化柔度
    for ely = 1:nely
        for elx = 1:nelx                %这里是两层for循环 相当于矩阵的遍历
            n1 = (nely+1)*(elx-1)+ely;  %n1 代表矩阵中一个网  格的左上角
            n2 = (nely+1)*elx+ely;      %n2 代表矩阵中一个网格的右上角   
            %这里不懂的话可以去b站参考视频  https://www.bilibili.com/video/BV1h7421N7K3/?spm_id_from=333.788.videopod.sections&vd_source=02ffd836310dcb034a8bb8e545cdce63
            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(ely,elx)^penal*Ue'*KE*Ue;  %最小化柔度 即 目标函数表达式  这个是公式
            dc(ely,elx) = -penal*x(ely,elx)^(penal-1)*Ue'*KE*Ue;  %灵敏度分析表达式  对上面式子求一次微分
        end
    end

    % FILTERING OF SENSITIVITIES  灵敏度过滤
    [dc] = check(nelx,nely,rmin,x,dc);     %调用（网格依赖性过滤）子程序（49——64）

    % DESIGN UPDATE BY THE OPTIMALITY CRITERIA METHOD
    [x] = OC(nelx,nely,x,volfrac,dc);      %调用（优化准则法）子程序（37——48）

    % PRINT RESULTS   打印结果
    change = max(max(abs(x-xold)));        %设计变量改变值的绝对值的最大值   双max是因为二维
    %disp为数据显示函数
    %...就是接上行的意思
    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 plot densiwies  图的显示
    colormap(gray); imagesc(-x); axis equal; axis tight; axis off; pause(1e-6);
    %colormap(gray);  （灰度） 色图
    %imagesc(-x);  使用缩放颜色显示图像
    %axis equal；  坐标轴单位长度相同
    %axis tight；  坐标轴的范围为数据范围
    %axis off；    关闭所有坐标轴线，刻度标记，标签
    %pause(1e-6);  暂停延时（1e-6）
    
end
%基于优化准则法的求解器（37——48）
%%%%%%%%%% OPTIMALITY CRITERIA UPDATE %%%%%%%%%
function [xnew]=OC(nelx,nely,x,volfrac,dc)  
l1 = 0; l2 = 100000; move = 0.2; % l1：二分法的上边界  l2：二分法的下边界  move：步长
while (l2-l1) > 1e-4    %循环终止条件：两个边界的插值小于1e-4
    lmid = 0.5*(l2+l1); %二分法发的中间值
    xnew = max(0.001,max(x-move,min(1.,min(x+move,x.*sqrt(-dc./lmid)))));%这句要是不明白 也可以去b站看那个视频
    %x-move 是变化的下限     x+move 是变化的下限
    %0.001 是相对密度的下限  1 是相对密度的上限
    %此句子是在判断x.*sqrt(-dc./lmid)是否位于x-move 和 x+move 之间，否则就取上下限
    %这样一个数轴表示   0.001  x-move  x+move  1 
    if sum(sum(xnew)) - volfrac*nelx*nely > 0   %说明区间的值变大了  说明lmid小了 （分母变小整体增大）
        l1 = lmid;                              %所以要把lmid增大
    else
        l2 = lmid;
    end
end
%网格依赖性过滤 49-64
%%%%%%%%%% MESH-INDEPENDENCY FILTER %%%%%%%%%%%
    function [dcn]=check(nelx,nely,rmin,x,dc) %rmin 过滤半径  x 相对密度  dc 敏度
dcn = zeros(nely,nelx);  %清零处理，保存更新的目标函数灵敏度
for i = 1:nelx
    for j = 1:nely
        sum = 0.0;
        for k = max(i-round(rmin),1):min(i+round(rmin),nelx)
            %round(rmin)：将过滤半径舍入到最近的整数
            for l = max(j-round(rmin),1):min(j+round(rmin),nely)
                %主循环已经确定了以某个单元坐标（j,i）为圆心，遍历附近以rmin为半径的单元，圆内以及圆上包含其他单元的中心坐标，将其添加到加权叠加灵敏度计算
                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
%end
%有限元分析  65-99
%%%%%%%%%% FE-ANALYSIS %%%%%%%%%%%%
function [U]=FE(nelx,nely,x,penal)  %节点位移矩阵
[KE] = lk;                          %计算单元刚度矩阵
%spare函数用于将矩阵的任意0元素去除，非0元素以及下标组成新的矩阵
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); %力矩阵和位移矩阵的稀疏矩阵，F=KU

for elx = 1:nelx
    for ely = 1:nely
        n1 = (nely+1)*(elx-1)+ely;
        n2 = (nely+1)*elx+ely;      %同上
        %将单元刚度矩阵组装成总体刚度矩阵K
        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
%end

% DEFINE LOADS AND SUPPORTS (HALF MBB-BEAM)
F(2,1) = -1;        %施加载荷  左上角竖直向下单位力
%union函数用于找到所有元素，去除掉重复元素之后升序排列
fixeddofs = union([1:2:2*(nely+1)], [2*(nelx+1)*(nely+1)]); %施加约束

alldofs = [1:2*(nely+1)*(nelx+1)];  %求所有自由度
freedofs = setdiff(alldofs,fixeddofs);  %返回第一个数组存在，但是第二个数组不存在的元素（求没有施加约束的自由度）

% SOLVING
U(freedofs,:) = K(freedofs,freedofs) \ F(freedofs,:);%求解线性方程组 得到各个节点的位移
U(fixeddofs,:) = 0;   %受约束节点的自由度位移量为0
%有限元分析  65-99
%%%%%%%%%% 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) ];
%end
%end