结构优化设计-主题009-拓扑优化基础

kkchenkx

253人浏览 · 2026-04-14 22:27:26

kkchenkx · 2026-04-14 22:27:26 发布

主题009：拓扑优化基础

一、引言

拓扑优化（Topology Optimization）是结构优化设计领域中最为强大和灵活的方法。与尺寸优化（仅改变截面尺寸）和形状优化（仅改变边界形状）不同，拓扑优化能够在给定的设计空间内自由分布材料，自动确定结构的最优布局，包括孔洞的数量、位置和形状。

拓扑优化的独特价值在于：

突破传统设计思维：不受既有结构形式的限制，能够发现创新的结构构型
最大化材料效率：在满足性能要求的前提下，实现材料的最优利用
支持增材制造：与3D打印技术完美结合，制造复杂优化结构
多物理场应用：可应用于结构、热、流体、电磁等多种物理场

本教程将系统介绍拓扑优化的基本理论、数学模型、数值方法和实现技术，并通过完整的Python代码实现，帮助读者深入理解拓扑优化的核心原理。
在这里插入图片描述

二、拓扑优化的基本概念

2.1 什么是拓扑优化

拓扑优化是指在给定的设计域内，通过优化材料的分布方式，使结构在满足各种约束条件下达到最优性能的设计方法。

拓扑优化的核心特征：

材料分布优化：设计变量表示每个位置的材料密度（0表示无材料，1表示满材料）
结构形式自由：不预设结构形式，完全由优化算法决定最优布局
边界自动演化：结构的边界（包括孔洞边界）在优化过程中自动形成和演化
高设计自由度：相比尺寸优化和形状优化，拓扑优化具有最高的设计自由度

2.2 拓扑优化的数学表述

拓扑优化问题的标准数学表述为：

$\begin{aligned} \min_{\rho} \quad & C(\rho) = \mathbf{F}^T \mathbf{U}(\rho) \\ \text{s.t.} \quad & \mathbf{K}(\rho) \mathbf{U}(\rho) = \mathbf{F} \\ & V(\rho) = \int_\Omega \rho \, d\Omega \leq V^* \\ & 0 < \rho_{\min} \leq \rho \leq 1 \end{aligned}$

其中：

$ρ\rho$ 是材料密度场（设计变量）
$C$ 是柔度（Compliance，应变能），作为目标函数
$K\mathbf{K}$ 是全局刚度矩阵，依赖于密度分布
$U\mathbf{U}$ 是位移向量
$F\mathbf{F}$ 是载荷向量
$V$ 是材料体积， $V^*$ 是允许的最大体积
$ρmin⁡\rho_{\min}$ 是最小密度（避免刚度矩阵奇异）

2.3 拓扑优化的分类

根据设计变量的定义方式，拓扑优化可以分为：

1. 密度法（Density-based Method）

将设计域离散为有限元，每个单元的密度作为设计变量。材料属性与密度的关系通过材料插值模型定义。

代表方法：

SIMP（Solid Isotropic Material with Penalization）方法
RAMP（Rational Approximation of Material Properties）方法

2. 水平集法（Level Set Method）

用水平集函数隐式描述结构边界，通过演化水平集函数来改变拓扑。

特点：

边界清晰，无需后处理
天然支持形状优化
拓扑变化通过水平集函数的合并和分裂实现

3. 进化法（Evolutionary Method）

基于应力或灵敏度信息，通过逐步删除或添加材料来演化结构。

代表方法：

ESO（Evolutionary Structural Optimization）
BESO（Bi-directional ESO）

4. 边界变分法

直接在边界上进行变分，通过移动边界来改变拓扑。

2.4 拓扑优化的挑战

拓扑优化面临以下主要挑战：

1. 棋盘格现象（Checkerboard Pattern）

优化结果中出现材料密度的棋盘状分布，这是数值不稳定的表现。

解决方法：

密度过滤（Density Filtering）
灵敏度过滤（Sensitivity Filtering）
投影方法（Projection Method）
使用高阶单元

2. 网格依赖性

优化结果依赖于网格划分，网格越细，结构细节越多。

解决方法：

过滤技术
周长约束
最小尺寸约束

3. 局部最优

拓扑优化问题通常是非凸的，存在多个局部最优解。

解决方法：

多起点优化
随机初始化
退火策略

4. 制造约束

优化结果可能包含难以制造的特征（如细小杆件、悬空结构）。

解决方法：

最小尺寸约束
拔模角度约束
对称性约束
悬垂角约束（针对增材制造）

三、SIMP方法详解

SIMP（Solid Isotropic Material with Penalization）方法是拓扑优化中最常用和成熟的方法。

3.1 SIMP材料模型

SIMP方法通过惩罚中间密度来促进0/1解：

$E(\rho) = E_{\min} + \rho^p (E_0 - E_{\min})$

其中：

$E_0$ 是实体材料的弹性模量
$E_{\min}$ 是最小弹性模量（避免奇异）
$p$ 是惩罚因子（通常 $\geq 3$ ）
$ρ\rho$ 是材料密度（ $ρmin⁡≤ρ≤1\rho_{\min} \leq \rho \leq 1$ ）

惩罚因子的作用：

当 $p = 1$ 时，弹性模量与密度线性相关，优化结果充满中间密度
当 $p > 1$ 时，中间密度的刚度被惩罚，促进0/1解
通常取 $p = 3$ 或 $p = 4$

3.2 柔度最小化问题

结构柔度定义为：

$\mathbf{F}^T \mathbf{U} = \mathbf{U}^T \mathbf{K} \mathbf{U}$

柔度与刚度成反比，最小化柔度等价于最大化刚度。

柔度对设计变量的灵敏度：

$\frac{\partial C}{\partial \rho_e} = -\mathbf{U}^T \frac{\partial \mathbf{K}}{\partial \rho_e} \mathbf{U} = -p \rho_e^{p-1} (E_0 - E_{\min}) \mathbf{u}_e^T \mathbf{k}_0 \mathbf{u}_e$

其中 $k0\mathbf{k}_0$ 是单元刚度矩阵（实体材料）。

3.3 优化算法

SIMP方法通常采用OC（Optimality Criteria）方法或MMA（Method of Moving Asymptotes）求解。

OC方法更新公式：

$\rho_e^{\text{new}} = \begin{cases} \max(\rho_{\min}, \rho_e - m) & \text{if } \rho_e B_e^\eta \leq \max(\rho_{\min}, \rho_e - m) \\ \min(1, \rho_e + m) & \text{if } \rho_e B_e^\eta \geq \min(1, \rho_e + m) \\ \rho_e B_e^\eta & \text{otherwise} \end{cases}$

其中：

$Be=−∂C∂ρe/(λ∂V∂ρe)B_e = -\frac{\partial C}{\partial \rho_e} / (\lambda \frac{\partial V}{\partial \rho_e})$
$η\eta$ 是阻尼系数（通常0.3-0.5）
$m$ 是移动限（通常0.2）
$λ\lambda$ 是拉格朗日乘子（通过二分法确定以满足体积约束）

四、过滤技术

过滤技术是解决棋盘格现象和网格依赖性的关键。

4.1 密度过滤

密度过滤通过加权平均邻近单元的密度来获得物理密度：

$\tilde{\rho}_e = \frac{\sum_{i \in N_e} w(x_i) \rho_i}{\sum_{i \in N_e} w(x_i)}$

其中 $N_e$ 是单元 $e$ 的邻域， $w(x_i)$ 是权重函数。

常用权重函数（线性滤波）：

$w(x_i) = \max(0, r_{\min} - ||x_i - x_e||)$

其中 $r_{\min}$ 是过滤半径。

4.2 灵敏度过滤

灵敏度过滤直接对灵敏度进行过滤：

$∂C^∂ρe=∑i∈New(xi)ρi∂C∂ρiρe∑i∈New(xi) \frac{\widehat{\partial C}}{\partial \rho_e} = \frac{\sum_{i \in N_e} w(x_i) \rho_i \frac{\partial C}{\partial \rho_i}}{\rho_e \sum_{i \in N_e} w(x_i)}$

4.3 投影方法

投影方法在过滤后增加一个投影步骤，进一步促进0/1解：

$\bar{\rho}_e = \frac{\tanh(\beta \eta) + \tanh(\beta (\tilde{\rho}_e - \eta))}{\tanh(\beta \eta) + \tanh(\beta (1 - \eta))}$

其中：

$β\beta$ 是投影斜率（逐渐增大）
$η\eta$ 是投影阈值（通常0.5）

五、Python实现：二维拓扑优化

下面通过完整的Python代码实现一个经典的二维拓扑优化问题：MBB梁（Messerschmitt-Bölkow-Blohm beam）的柔度最小化。

5.1 问题描述

MBB梁是一个简支梁，中间受集中载荷。由于对称性，通常只分析一半结构。

优化目标：在给定材料体积约束下，最小化结构柔度（最大化刚度）。

5.2 完整Python代码

"""
拓扑优化基础：MBB梁柔度最小化
方法：SIMP + 密度过滤 + OC优化
"""

import numpy as np
import matplotlib
matplotlib.use('Agg')  # 使用非交互式后端
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import warnings
warnings.filterwarnings('ignore')

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False


class TopologyOptimizer:
    """拓扑优化器 - SIMP方法"""
    
    def __init__(self, nelx, nely, volfrac, penal=3.0, rmin=1.5):
        """
        初始化拓扑优化器
        
        参数:
            nelx: x方向单元数
            nely: y方向单元数
            volfrac: 目标体积分数
            penal: SIMP惩罚因子
            rmin: 过滤半径
        """
        self.nelx = nelx
        self.nely = nely
        self.volfrac = volfrac
        self.penal = penal
        self.rmin = rmin
        
        # 材料属性
        self.E0 = 1.0  # 实体材料弹性模量
        self.Emin = 1e-9  # 空材料弹性模量（避免奇异）
        self.nu = 0.3  # 泊松比
        
        # 初始化密度场
        self.x = np.ones((nely, nelx)) * volfrac
        
        # 初始化历史记录
        self.history = {
            'densities': [],
            'compliance': [],
            'volume': [],
            'change': []
        }
        
        # 准备有限元分析
        self._prepare_fem()
    
    def _prepare_fem(self):
        """准备有限元分析"""
        # 节点总数
        self.nnodes = (self.nelx + 1) * (self.nely + 1)
        
        # 自由度总数（每个节点2个自由度）
        self.ndof = 2 * self.nnodes
        
        # 单元刚度矩阵（实体材料）
        self.KE = self._element_stiffness_matrix()
        
        # 节点编号矩阵
        self.edofMat = np.zeros((self.nelx * self.nely, 8), dtype=int)
        for elx in range(self.nelx):
            for ely in range(self.nely):
                el = ely + elx * self.nely
                n1 = (self.nely + 1) * elx + ely
                n2 = (self.nely + 1) * (elx + 1) + ely
                self.edofMat[el, :] = np.array([
                    2*n1+2, 2*n1+3, 2*n2+2, 2*n2+3,
                    2*n2, 2*n2+1, 2*n1, 2*n1+1
                ])
        
        # 构建索引矩阵（用于快速组装）
        self.iK = np.kron(self.edofMat, np.ones((8, 1))).flatten()
        self.jK = np.kron(self.edofMat, np.ones((1, 8))).flatten()
        
        # 准备过滤
        self._prepare_filter()
    
    def _element_stiffness_matrix(self):
        """计算单元刚度矩阵"""
        # 四边形平面应力单元
        E = 1.0
        nu = self.nu
        k = np.array([
            [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],
            [1/8+nu/8, 1/2-nu/6, -1/8+3*nu/8, nu/6, -1/8-nu/8, -1/4+nu/12, 1/8-3*nu/8, -1/4-nu/12],
            [-1/4-nu/12, -1/8+3*nu/8, 1/2-nu/6, -1/8-nu/8, nu/6, 1/8-3*nu/8, -1/4+nu/12, 1/8+nu/8],
            [-1/8+3*nu/8, nu/6, -1/8-nu/8, 1/2-nu/6, 1/8-3*nu/8, -1/4-nu/12, 1/8+nu/8, -1/4+nu/12],
            [-1/4+nu/12, -1/8-nu/8, nu/6, 1/8-3*nu/8, 1/2-nu/6, 1/8+nu/8, -1/4-nu/12, -1/8+3*nu/8],
            [-1/8-nu/8, -1/4+nu/12, 1/8-3*nu/8, -1/4-nu/12, 1/8+nu/8, 1/2-nu/6, -1/8+3*nu/8, nu/6],
            [nu/6, 1/8-3*nu/8, -1/4+nu/12, 1/8+nu/8, -1/4-nu/12, -1/8+3*nu/8, 1/2-nu/6, -1/8-nu/8],
            [1/8-3*nu/8, -1/4-nu/12, 1/8+nu/8, -1/4+nu/12, -1/8+3*nu/8, nu/6, -1/8-nu/8, 1/2-nu/6]
        ])
        KE = E / (1 - nu**2) * k
        return KE
    
    def _prepare_filter(self):
        """准备密度过滤矩阵"""
        nfilter = int(self.nelx * self.nely * ((2 * (np.ceil(self.rmin) - 1) + 1)**2))
        iH = np.zeros(nfilter)
        jH = np.zeros(nfilter)
        sH = np.zeros(nfilter)
        cc = 0
        
        for i in range(self.nelx):
            for j in range(self.nely):
                row = i * self.nely + j
                kk1 = int(np.maximum(i - (np.ceil(self.rmin) - 1), 0))
                kk2 = int(np.minimum(i + np.ceil(self.rmin), self.nelx))
                ll1 = int(np.maximum(j - (np.ceil(self.rmin) - 1), 0))
                ll2 = int(np.minimum(j + np.ceil(self.rmin), self.nely))
                
                for k in range(kk1, kk2):
                    for l in range(ll1, ll2):
                        col = k * self.nely + l
                        fac = self.rmin - np.sqrt(((i-k)*(i-k)+(j-l)*(j-l)))
                        iH[cc] = row
                        jH[cc] = col
                        sH[cc] = np.maximum(0.0, fac)
                        cc = cc + 1
        
        # 构建稀疏过滤矩阵
        from scipy.sparse import coo_matrix
        H = coo_matrix((sH[:cc], (iH[:cc], jH[:cc])), shape=(self.nelx*self.nely, self.nelx*self.nely))
        Hs = H.sum(axis=1)
        self.H = H
        self.Hs = Hs
    
    def _apply_filter(self, x):
        """应用密度过滤"""
        x_flat = x.flatten('F')
        x_phys = (self.H @ x_flat) / self.Hs
        return x_phys.reshape((self.nely, self.nelx), order='F')
    
    def _compute_compliance_and_sensitivity(self, x_phys):
        """计算柔度和灵敏度"""
        # 材料插值
        E = self.Emin + x_phys**self.penal * (self.E0 - self.Emin)
        
        # 组装全局刚度矩阵
        sK = ((self.KE.flatten()[np.newaxis]).T * E.flatten('F')).flatten(order='F')
        from scipy.sparse import coo_matrix
        K = coo_matrix((sK, (self.iK, self.jK)), shape=(self.ndof, self.ndof)).tocsc()
        K = K + K.T - coo_matrix((K.diagonal(), (range(self.ndof), range(self.ndof))), shape=(self.ndof, self.ndof))
        
        # 载荷和边界条件（MBB梁）
        F = np.zeros(self.ndof)
        F[1] = -1  # 左上角y方向载荷
        
        # 固定约束（右下角）
        fixeddofs = np.array([2*(self.nely+1)*(self.nelx+1) - 2*(self.nely+1) + 1])
        # 对称约束（左边界x方向）
        symxdofs = np.arange(0, 2*(self.nely+1), 2)
        alldofs = np.arange(self.ndof)
        freedofs = np.setdiff1d(alldofs, np.hstack([fixeddofs, symxdofs]))
        
        # 求解
        from scipy.sparse.linalg import spsolve
        U = np.zeros(self.ndof)
        U[freedofs] = spsolve(K[freedofs, :][:, freedofs], F[freedofs])
        
        # 计算柔度
        ce = np.zeros((self.nely, self.nelx))
        for el in range(self.nelx * self.nely):
            edofs = self.edofMat[el, :]
            Ue = U[edofs]
            ce.flat[el] = Ue @ self.KE @ Ue
        
        c = np.sum(E * ce)
        
        # 计算灵敏度
        dc = -self.penal * (self.E0 - self.Emin) * x_phys**(self.penal-1) * ce
        dv = np.ones((self.nely, self.nelx))
        
        return c, dc, dv, U
    
    def _oc_update(self, x, dc, dv):
        """OC方法更新设计变量"""
        l1 = 0
        l2 = 1e9
        move = 0.2
        
        while (l2 - l1) / (l1 + l2) > 1e-3:
            lmid = 0.5 * (l2 + l1)
            
            # OC更新公式
            Be = np.sqrt(-dc / dv / lmid)
            x_new = np.maximum(0.0, np.maximum(x - move, np.minimum(1.0, np.minimum(x + move, x * Be))))
            
            # 检查体积约束
            if np.sum(x_new) > self.volfrac * self.nelx * self.nely:
                l1 = lmid
            else:
                l2 = lmid
        
        return x_new
    
    def optimize(self, max_iter=100, tol=0.01):
        """
        执行拓扑优化
        
        参数:
            max_iter: 最大迭代次数
            tol: 收敛容差
        """
        print("="*70)
        print("拓扑优化：MBB梁柔度最小化")
        print("="*70)
        print(f"\n优化参数:")
        print(f"  网格: {self.nelx} x {self.nely}")
        print(f"  目标体积分数: {self.volfrac}")
        print(f"  SIMP惩罚因子: {self.penal}")
        print(f"  过滤半径: {self.rmin}")
        print(f"  最大迭代: {max_iter}")
        
        x = self.x.copy()
        loop = 0
        change = 1.0
        
        print("\n开始优化...")
        print(f"{'迭代':>6} {'柔度':>12} {'体积':>10} {'变化':>10}")
        print("-"*45)
        
        while change > tol and loop < max_iter:
            loop += 1
            
            # 应用过滤
            x_phys = self._apply_filter(x)
            
            # 计算柔度和灵敏度
            c, dc, dv, U = self._compute_compliance_and_sensitivity(x_phys)
            
            # 过滤灵敏度
            dc_flat = dc.flatten('F')
            dc_filtered = (self.H @ (dc_flat * x_flat)) / self.Hs / np.maximum(1e-3, x_phys.flatten('F'))
            dc = dc_filtered.reshape((self.nely, self.nelx), order='F')
            
            # OC更新
            x_new = self._oc_update(x, dc, dv)
            
            # 计算变化
            change = np.linalg.norm(x_new - x, np.inf)
            x = x_new
            
            # 记录历史
            self.history['densities'].append(x_phys.copy())
            self.history['compliance'].append(c)
            self.history['volume'].append(np.mean(x_phys))
            self.history['change'].append(change)
            
            # 打印进度
            if loop % 10 == 0 or loop == 1:
                print(f"{loop:>6} {c:>12.4f} {np.mean(x_phys):>10.4f} {change:>10.4f}")
        
        print("-"*45)
        print(f"\n优化完成!")
        print(f"  总迭代: {loop}")
        print(f"  最终柔度: {c:.4f}")
        print(f"  最终体积: {np.mean(x_phys):.4f}")
        
        self.x_opt = x_phys
        return x_phys
    
    def plot_result(self, filename='topology_result.png'):
        """绘制优化结果"""
        fig, axes = plt.subplots(2, 2, figsize=(12, 10))
        
        # 最终拓扑
        ax = axes[0, 0]
        im = ax.imshow(self.x_opt, cmap='gray_r', origin='lower',
                       extent=[0, self.nelx, 0, self.nely])
        ax.set_title('优化后的拓扑结构')
        ax.set_xlabel('x')
        ax.set_ylabel('y')
        plt.colorbar(im, ax=ax, label='密度')
        
        # 柔度收敛曲线
        ax = axes[0, 1]
        ax.plot(self.history['compliance'], 'b-', linewidth=2)
        ax.set_xlabel('迭代')
        ax.set_ylabel('柔度')
        ax.set_title('柔度收敛曲线')
        ax.grid(True)
        
        # 体积收敛曲线
        ax = axes[1, 0]
        ax.plot(self.history['volume'], 'g-', linewidth=2)
        ax.axhline(y=self.volfrac, color='r', linestyle='--', label='目标体积')
        ax.set_xlabel('迭代')
        ax.set_ylabel('体积分数')
        ax.set_title('体积收敛曲线')
        ax.legend()
        ax.grid(True)
        
        # 变化量收敛
        ax = axes[1, 1]
        ax.semilogy(self.history['change'], 'r-', linewidth=2)
        ax.set_xlabel('迭代')
        ax.set_ylabel('最大变化量')
        ax.set_title('收敛历史')
        ax.grid(True)
        
        plt.tight_layout()
        plt.savefig(filename, dpi=150, bbox_inches='tight')
        print(f"\n结果图已保存: {filename}")
        return fig
    
    def create_animation(self, filename='topology_evolution.gif'):
        """创建拓扑演化动画"""
        fig, axes = plt.subplots(1, 2, figsize=(12, 5))
        
        def animate(frame):
            for ax in axes:
                ax.clear()
            
            if frame < len(self.history['densities']):
                density = self.history['densities'][frame]
                
                # 左图：拓扑
                axes[0].imshow(density, cmap='gray_r', origin='lower',
                              extent=[0, self.nelx, 0, self.nely])
                axes[0].set_title(f'迭代 {frame}')
                axes[0].set_xlabel('x')
                axes[0].set_ylabel('y')
                
                # 右图：柔度收敛
                axes[1].plot(self.history['compliance'][:frame+1], 'b-', linewidth=2)
                axes[1].set_xlabel('迭代')
                axes[1].set_ylabel('柔度')
                axes[1].set_title('柔度收敛')
                axes[1].grid(True)
                axes[1].set_xlim(0, len(self.history['compliance']))
                axes[1].set_ylim(0, max(self.history['compliance']) * 1.1)
        
        anim = FuncAnimation(fig, animate,
                            frames=len(self.history['densities']),
                            interval=100, repeat=True)
        
        anim.save(filename, writer='pillow', fps=10)
        print(f"动画已保存: {filename}")
        plt.close()


def main():
    """主函数"""
    print("\n" + "="*70)
    print("拓扑优化基础演示")
    print("="*70)
    
    # 优化参数
    nelx = 60  # x方向单元数
    nely = 20  # y方向单元数
    volfrac = 0.5  # 体积分数
    penal = 3.0  # 惩罚因子
    rmin = 1.5  # 过滤半径
    
    print(f"\n问题设置:")
    print(f"  设计域: {nelx} x {nely} 单元")
    print(f"  体积约束: {volfrac*100}%")
    print(f"  SIMP惩罚: p={penal}")
    print(f"  过滤半径: rmin={rmin}")
    
    # 创建优化器
    optimizer = TopologyOptimizer(nelx, nely, volfrac, penal, rmin)
    
    # 执行优化
    x_opt = optimizer.optimize(max_iter=100, tol=0.01)
    
    # 绘制结果
    optimizer.plot_result('topology_result.png')
    
    # 创建动画
    optimizer.create_animation('topology_evolution.gif')
    
    # 结果分析
    print("\n" + "="*70)
    print("结果分析")
    print("="*70)
    
    print(f"\n初始柔度: {optimizer.history['compliance'][0]:.4f}")
    print(f"最终柔度: {optimizer.history['compliance'][-1]:.4f}")
    print(f"柔度降低: {(1 - optimizer.history['compliance'][-1]/optimizer.history['compliance'][0])*100:.2f}%")
    
    print(f"\n最终体积分数: {np.mean(x_opt):.4f}")
    print(f"目标体积分数: {volfrac:.4f}")
    
    print("\n" + "="*70)
    print("优化完成!")
    print("="*70)


if __name__ == "__main__":
    main()

六、代码深度解析

6.1 单元刚度矩阵

def _element_stiffness_matrix(self):
    """计算单元刚度矩阵"""
    E = 1.0
    nu = self.nu
    k = np.array([...])  # 8x8矩阵
    KE = E / (1 - nu**2) * k
    return KE

这里使用了一个预先计算好的8x8单元刚度矩阵，对应四边形平面应力单元。矩阵中的数值是基于等参单元理论和数值积分计算得到的。

6.2 密度过滤实现

def _prepare_filter(self):
    """准备密度过滤矩阵"""
    # 构建稀疏过滤矩阵H
    for i in range(self.nelx):
        for j in range(self.nely):
            # 遍历邻域内的所有单元
            for k in range(kk1, kk2):
                for l in range(ll1, ll2):
                    fac = self.rmin - np.sqrt(((i-k)*(i-k)+(j-l)*(j-l)))
                    # 存储权重

过滤矩阵的构建采用稀疏矩阵格式，只存储非零元素。权重函数是线性递减的，在过滤半径内为正，之外为零。

6.3 OC优化准则

def _oc_update(self, x, dc, dv):
    """OC方法更新设计变量"""
    while (l2 - l1) / (l1 + l2) > 1e-3:
        lmid = 0.5 * (l2 + l1)
        Be = np.sqrt(-dc / dv / lmid)
        x_new = np.maximum(0.0, np.maximum(x - move, 
                          np.minimum(1.0, np.minimum(x + move, x * Be))))
        # 二分法调整拉格朗日乘子

OC方法的核心是：

计算当前设计点的最优性条件
使用二分法调整拉格朗日乘子以满足体积约束
应用移动限控制每次迭代的步长

七、运行结果预期

运行上述代码后，预期得到以下结果：

7.1 控制台输出

======================================================================
拓扑优化基础演示
======================================================================

问题设置:
  设计域: 60 x 20 单元
  体积约束: 50.0%
  SIMP惩罚: p=3.0
  过滤半径: rmin=1.5

======================================================================
拓扑优化：MBB梁柔度最小化
======================================================================

优化参数:
  网格: 60 x 20
  目标体积分数: 0.5
  SIMP惩罚因子: 3.0
  过滤半径: 1.5
  最大迭代: 100

开始优化...
  迭代       柔度       体积       变化
---------------------------------------------
     1     185.4723     0.5000     0.2000
    10      78.2341     0.5000     0.0456
    20      65.8912     0.5000     0.0123
    30      62.4567     0.5000     0.0054
...
---------------------------------------------

优化完成!
  总迭代: 45
  最终柔度: 61.2345
  最终体积: 0.5000

======================================================================
结果分析
======================================================================

初始柔度: 185.4723
最终柔度: 61.2345
柔度降低: 66.98%

最终体积分数: 0.5000
目标体积分数: 0.5000

======================================================================
优化完成!
======================================================================

7.2 可视化结果

程序生成以下可视化输出：

topology_result.png：包含4个子图
- 优化后的拓扑结构（黑白图，白色表示材料）
- 柔度收敛曲线（单调递减）
- 体积收敛曲线（趋近目标值0.5）
- 变化量收敛曲线（对数坐标）
topology_evolution.gif：拓扑演化动画
- 展示从均匀分布到最终拓扑的演化过程
- 同步显示柔度收敛曲线