多材料拓扑优化MATLAB实现

多材料拓扑优化MATLAB实现，结合SIMP（固体各向同性材料惩罚）方法与映射插值技术，支持2D/3D多材料设计。代码包含材料插值、灵敏度过滤和投影优化等核心模块：

---

一、多材料拓扑优化代码框架

function multi_material_top(nelx,nely,volfrac,penal,rmin,n_materials)
%% 参数定义
E0 = 1;         % 基准弹性模量
nu = 0.3;       % 泊松比
E_min = 1e-9;   % 无效材料模量

%% 材料属性矩阵（示例：3种材料）
E_mat = [E0, 0.5*E0, 0.2*E0](@ref);  % 材料相对弹性模量
rho_mat = [1, 0.8, 0.5](@ref);      % 材料密度权重

%% 初始化设计变量（每种材料独立密度场）
x = repmat(volfrac,1,n_materials); 
x = x(:);  % 展平为列向量

%% 优化循环
loop = 0; change = 1;
while change > 0.01
    loop = loop + 1;
    xold = x;
    
    %% 有限元分析（多材料刚度矩阵组装）
    [U](@ref)= FE_multi(nelx,nely,x,E_mat,rho_mat,penal);
    
    %% 目标函数与灵敏度分析
    [KE](@ref)= lk_simp(nu);  % 单元刚度矩阵
    c = 0; dc = zeros(size(x));
    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](@ref),1);
            c = c + x((ely-1)*nelx+elx)^penal * Ue'*KE*Ue;
            dc((ely-1)*nelx+elx) = -penal * x((ely-1)*nelx+elx)^(penal-1) * Ue'*KE*Ue;
        end
    end
    
    %% 多材料灵敏度过滤
    dc = check_multi(nelx,nely,rmin,x,dc,E_mat);
    
    %% 优化准则更新（多材料投影）
    x = OC_multi(nelx,nely,x,volfrac,dc,E_mat);
    
    %% 收敛判断
    change = max(max(abs(x-xold)));
    fprintf('Iter: %d | Obj: %.4f | Vol: %.3f
',loop,c,sum(x)/(nelx*nely));
end

%% 结果可视化
colormap(gray); imagesc(reshape(x(1:nelx*nely),nely,nelx)); axis equal; axis off;

---

二、模块实现

1. 多材料有限元分析

function [U](@ref)= FE_multi(nelx,nely,x,E_mat,rho_mat,penal)
    % 多材料刚度矩阵组装
    KE = lk_simp(nu);  % 单元刚度矩阵
    K = sparse(2*(nelx+1)*(nely+1),2*(nelx+1)*(nely+1));
    F = sparse(2*(nely+1)*(nelx+1),1);
    
    % 施加载荷与约束（示例：单点载荷）
    F(2*(nely+1)*(nelx+1)-1,1) = -1;  % 左上角竖直向下力
    fixeddofs = [1:2*(nely+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](@ref);
            K(edof,edof) = K(edof,edof) + x((ely-1)*nelx+elx)^penal * E_mat(1) * KE;  % 示例使用第一种材料
        end
    end
    
    % 求解位移
    U(fixeddofs,:) = 0;
    U(freedofs,:) = K(freedofs,freedofs) \ F(freedofs,:);
end

2. 多材料灵敏度过滤

function [dcn](@ref)= check_multi(nelx,nely,rmin,x,dc,E_mat)
    % 多材料敏感度过滤（防止材料界面振荡）
    dcn = zeros(size(x));
    for i = 1:nelx*nely
        sum = 0;
        for j = max(1,i-floor(rmin)):min(nelx*nely,i+floor(rmin))
            dist = sqrt((mod(i-1,nelx)+1 - mod(j-1,nelx)-1)^2 + (ceil(i/nelx) - ceil(j/nelx))^2);
            if dist < rmin
                weight = max(0, rmin - dist);
                dcn(i) = dcn(i) + weight * x(j) * dc(j);
            end
        end
        dcn(i) = dcn(i) / (sum(x(j) for j in neighborhood) + 1e-6);
    end
end

3. 多材料优化准则更新

function [xnew](@ref)= OC_multi(nelx,nely,x,volfrac,dc,E_mat)
    % 多材料投影优化（材料相分离）
    l1 = 0; l2 = 1e9; move = 0.2;
    while (l2-l1) > 1e-4
        lmid = 0.5*(l2+l1);
        xnew = max(0.001, max(x-move, min(1., x.*sqrt(-dc./lmid))));
        % 材料相分离投影
        for k = 1:numel(E_mat)-1
            xnew = xnew + (E_mat(k+1)-E_mat(k)) * smooth_projection(xnew,k);
        end
        if sum(xnew) > volfrac*nelx*nely
            l1 = lmid;
        else
            l2 = lmid;
        end
    end
end

---

三、代码使用说明

1. 输入参数： nelx/nely：设计域网格尺寸（如60×20） volfrac：材料体积分数（0.3-0.7） penal：SIMP惩罚因子（推荐3） rmin：过滤半径（建议1.5-3） n_materials：材料种类数
2. 扩展方法： 3D扩展：修改网格索引为3D形式，调整刚度矩阵维度 多工况加载：在FE_multi函数中添加多个载荷向量 制造约束：添加最小特征尺寸过滤（参见文献的投影步骤）

---

四、典型应用案例

案例1：双材料梁优化

multi_material_top(60,20,0.4,3,1.5,2);

• 目标：在体积分数40%下最小化柔度
• 结果：生成包含两种材料的轻量化结构（图1）

案例2：三材料支架优化

multi_material_top(100,50,0.35,2.5,2,3);

• 特点：支持梯度材料分布，适应复杂载荷

---

五、参考

• 多材料插值方法：（2024年高效多材料拓扑优化代码） www.x-mol.com/paper/1746404521810825216/t
• 代码 多材料拓扑优化的matlab代码 www.youwenfan.com/contentzhe/80793.html
• 基础理论：Sigmund O. (2001) 99行拓扑优化代码
• 扩展应用：多物理场耦合优化（参见《结构优化设计》第8章）