1. 柔性支撑结构拓扑优化基础

任务说明

本案例面向柔性支撑结构的初步设计，以经典简支梁为对象，建立拓扑优化模型。目标是在体积约束为50%的条件下，以结构最小柔度（即最大刚度）为优化目标，求解材料的最优分布，从而理解柔性支撑中“去除低效材料、保留传力路径”的基本思想。优化结果将呈现为清晰的桁架状材料布局，为后续柔性铰链或柔顺机构的设计提供初始构型参考。

问题定义与优化模型

设计域为长宽比2:1的矩形简支梁（尺寸2m×1m），左下端固定铰支，右下端滚动铰支，上边中点受竖直向下集中力F=1kN。材料为各向同性线弹性体，弹性模量E=1MPa（归一化），泊松比ν=0.3。采用SIMP（固体各向同性材料惩罚）插值模型，设计变量为每个单元的伪密度x∈[0,1]，惩罚因子p=3，体积约束为设计域总体积的50%。优化目标为最小化结构柔度c=U^TKU，即最大化结构整体刚度。

拓扑优化实现（88行经典SIMP代码）

% 简支梁拓扑优化 - 最小柔度, 50%体积约束
% 输入: nelx=60, nely=30 (单元数), volfrac=0.5, penal=3
function top_beam()
    nelx = 60; nely = 30; volfrac = 0.5; penal = 3;
    ft = 1; rmin = 1.5; % 灵敏度过滤半径

    % 有限元预处理
    [KE] = lk(); % 单元刚度矩阵
    [nodeN, eN] = size(KE); % 节点自由度数
    nnp = (nelx+1)*(nely+1); % 节点总数
    ndof = 2*nnp; % 总自由度
    F = sparse(ndof,1); U = zeros(ndof,1);
    % 载荷: 上边中点向下
    midNode = (nely+1)*(nelx/2) + nely/2 + 1;
    F(2*midNode-1,1) = 0; F(2*midNode,1) = -1000;
    % 边界条件: 左下固定铰支, 右下滚动铰支
    fixeddofs = [1, 2]; % 左下节点x,y
    fixeddofs = [fixeddofs, 2*(nnp) - 1]; % 右下节点x (滚动)
    alldofs = 1:ndof;
    freedofs = setdiff(alldofs, fixeddofs);

    % 初始化设计变量
    x = repmat(volfrac, nely, nelx);
    xPhys = x; % 物理密度(经滤波)
    loop = 0; change = 1;

    % 优化循环
    while change > 0.01 && loop < 200
        loop = loop + 1;
        % 有限元求解
        [U] = FE(nelx, nely, xPhys, penal, KE, F, freedofs, fixeddofs);
        % 目标函数和灵敏度
        [c, dc] = compliance(nelx, nely, xPhys, penal, KE, U);
        % 灵敏度过滤
        dc = check(nelx, nely, rmin, x, dc);
        % OC法更新设计变量
        xnew = OC(nelx, nely, x, volfrac, dc);
        change = max(abs(xnew(:) - x(:)));
        x = xnew; xPhys = x;
    end

    % 输出结果
    fprintf('迭代次数: %d, 最终柔度: %.4e\n', loop, c);
    colormap(gray); imagesc(-xPhys); axis equal; axis off;
end

% 单元刚度矩阵
function [KE] = lk()
    E = 1e6; 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

% 有限元组装与求解
function [U] = FE(nelx, nely, x, penal, KE, F, freedofs, fixeddofs)
    ndof = 2*(nelx+1)*(nely+1);
    K = sparse(ndof, ndof);
    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
    U = zeros(ndof,1);
    U(freedofs) = K(freedofs, freedofs) \ F(freedofs);
end

% 柔度与灵敏度
function [c, dc] = compliance(nelx, nely, x, penal, KE, U)
    c = 0; dc = zeros(nely, nelx);
    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];
            ue = U(edof);
            c = c + x(ely,elx)^penal * ue'*KE*ue;
            dc(ely,elx) = -penal * x(ely,elx)^(penal-1) * ue'*KE*ue;
        end
    end
end

% 灵敏度过滤
function [dcn] = check(nelx, nely, rmin, x, dc)
    dcn = zeros(nely, nelx);
    for i = 1:nelx
        for j = 1:nely
            sum = 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
end

% OC优化准则法
function [xnew] = OC(nelx, nely, x, volfrac, dc)
    l1 = 0; l2 = 1e9; move = 0.2;
    while (l2-l1)/(l1+l2) > 1e-3
        lmid = 0.5*(l2+l1);
        xnew = max(0, max(x-move, min(1, min(x+move, x.*sqrt(-dc./lmid)))));
        if sum(xnew(:)) - volfrac*nelx*nely > 0
            l1 = lmid;
        else
            l2 = lmid;
        end
    end
end

运行与验证

运行方式：在MATLAB中保存上述代码为top_beam.m，直接运行即可。需安装MATLAB基础环境，无需额外工具箱。

验证方法：观察优化迭代过程中柔度值单调下降并最终收敛。程序结束后显示灰度图像，其中白色区域（密度≈1）为保留材料，黑色区域（密度≈0）为去除材料。预期优化结果呈现对称的桁架结构：从加载点向两支座方向形成两条主要传力路径，中间区域材料被大量去除，整体呈“倒V”或拱形骨架分布。

判定标准：满足以下条件即判定案例通过：

迭代曲线收敛（change<0.01），最终体积分数严格为50%；
优化结果中无灰度单元（密度值接近0或1），无棋盘格现象；
材料分布关于梁中线左右对称，且载荷传递路径清晰连续。

结果解读与柔性支撑关联

优化得到的材料分布呈现典型的“倒V”形桁架结构，载荷通过斜撑直接传递至两端支座，中间大部分区域材料被去除。这正是柔性支撑设计的核心思想：在满足刚度需求的前提下，通过拓扑优化识别出不可省略的传力路径，将材料集中于此，其余区域留空或填充柔性铰链结构。该简支梁案例虽为刚性结构，但其优化流程与结果可直接迁移至柔性支撑的初始构型设计——将优化后的实体区域作为柔性铰链的连接节点，将去除区域作为运动释放空间，即可获得兼具承载能力与变形能力的柔性支撑雏形。

📚 课程目录