主题074：机械零件优化

Mechanical Parts Optimization

---

一、场景描述

机械零件是构成机械系统的基本单元，其结构设计的优劣直接影响整机性能、可靠性和制造成本。传统机械零件设计往往依赖经验公式和标准规范，导致结构冗余、材料浪费。本教程将介绍如何运用拓扑优化技术对机械零件进行轻量化设计，在保证强度和刚度的前提下，实现材料利用率最大化。

本教程将涵盖：

• 传动部件（齿轮、轴、联轴器）的结构优化
• 轴承座的轻量化设计
• 连接件（螺栓、法兰）的优化方法
• 考虑制造约束的工程化设计

---

二、数学/物理模型

2.1 机械零件力学模型

机械零件在工作过程中承受多种载荷，需要建立相应的力学模型：

静力学平衡方程：

$$\nabla \cdot \mathbf{\sigma }+\mathbf{f}=0$$

其中，$\mathbf{\sigma }$ 为应力张量，$\mathbf{f}$ 为体积力向量。

本构关系（线弹性）：

$$\mathbf{\sigma }=\mathbf{D}:\mathbf{\epsilon }$$

其中，$\mathbf{D}$ 为弹性矩阵，$\mathbf{\epsilon }$ 为应变张量。

应变-位移关系：

$$\mathbf{\epsilon }=\frac{1}{2}(\nabla \mathbf{u}+\nabla {\mathbf{u}}^T)$$

其中，$\mathbf{u}$ 为位移场。

2.2 SIMP材料插值模型

采用Solid Isotropic Material with Penalization (SIMP)方法进行材料插值：

$$E(\rho )=E_0\cdot {\rho }^p$$

其中：

• $E_0$ 为实体材料弹性模量
• $\rho \in [0,1]$ 为单元密度
• $p$ 为惩罚因子（通常取 $p\ge 3$）

2.3 拓扑优化数学表述

目标函数（最小化柔度）：

$$\underset{\rho }{\min }C(\rho )={\mathbf{F}}^T\mathbf{U}$$

约束条件：

$$\text{s.t.}\{\begin{array}{ll}\mathbf{K}(\rho )\mathbf{U}=\mathbf{F}& \text{(平衡方程)}\\ {\sum }_{e=1}^N{\rho }_e{v}_e\le V^{\ast }& \text{(体积约束)}\\ 0<{\rho }_{\min }\le {\rho }_e\le 1& \text{(密度边界)}\end{array}$$

2.4 灵敏度分析

目标函数灵敏度：

$$\frac{\partial C}{\partial {\rho }_e}=-{\mathbf{U}}^T\frac{\partial \mathbf{K}}{\partial {\rho }_e}\mathbf{U}=-p\cdot {\rho }_e^{p-1}\cdot {\mathbf{u}}_e^T{\mathbf{k}}_0{\mathbf{u}}_e$$

其中，${\mathbf{k}}_0$ 为实体单元刚度矩阵，${\mathbf{u}}_e$ 为单元节点位移向量。

2.5 疲劳强度分析

机械零件常承受循环载荷，需进行疲劳校核：

疲劳极限估算：

$${\sigma }_{-1}\approx 0.5{\sigma }_u\text{(钢材)}$$

$${\sigma }_{-1}\approx 0.4{\sigma }_u\text{(其他材料)}$$

Goodman修正公式：

$$\frac{{\sigma }_a}{{\sigma }_{-1}}+\frac{{\sigma }_m}{{\sigma }_u}=\frac{1}{n}$$

其中，${\sigma }_a$ 为应力幅，${\sigma }_m$ 为平均应力，$n$ 为安全系数。

---

三、环境准备

3.1 依赖库安装

pip install numpy matplotlib pillow imageio scipy

3.2 完整代码实现

"""
主题074：机械零件优化
Mechanical Parts Optimization

本代码实现机械零件的拓扑优化，包括：
1. 传动部件优化（齿轮、轴、联轴器）
2. 轴承座优化
3. 连接件优化（螺栓、法兰）
4. 考虑制造约束的优化
"""

import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.use('Agg')
from matplotlib.patches import Rectangle, Circle, FancyBboxPatch, Polygon, Arc
from matplotlib.collections import LineCollection
import matplotlib.patches as mpatches
from PIL import Image
import io
import warnings
warnings.filterwarnings('ignore')

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei', 'Arial Unicode MS']
plt.rcParams['axes.unicode_minus'] = False


class MechanicalMaterial:
    """
    机械工程材料模型
    
    包括结构钢、合金钢、铝合金、钛合金等
    """
    
    def __init__(self, material_type='steel_45'):
        self.material_type = material_type
        
        # 机械工程常用材料属性定义
        self.materials = {
            'steel_45': {
                'E': 210e3,           # 杨氏模量 (MPa)
                'nu': 0.30,           # 泊松比
                'rho': 7850,          # 密度 (kg/m³)
                'sigma_y': 355,       # 屈服强度 (MPa)
                'sigma_u': 600,       # 抗拉强度 (MPa)
                'HB': 197,            # 布氏硬度
                'name': '45号钢',
                'cost': 1.0,          # 相对成本
                'machinability': 0.7  # 可加工性指数
            },
            'steel_40cr': {
                'E': 210e3,
                'nu': 0.30,
                'rho': 7850,
                'sigma_y': 785,       # 调质后
                'sigma_u': 980,
                'HB': 241,
                'name': '40Cr合金钢',
                'cost': 1.8,
                'machinability': 0.6
            },
            'aluminum_6061': {
                'E': 69e3,
                'nu': 0.33,
                'rho': 2700,
                'sigma_y': 276,
                'sigma_u': 310,
                'HB': 95,
                'name': '6061铝合金',
                'cost': 2.5,
                'machinability': 0.9
            },
            'titanium_tc4': {
                'E': 114e3,
                'nu': 0.34,
                'rho': 4430,
                'sigma_y': 825,
                'sigma_u': 950,
                'HB': 334,
                'name': 'TC4钛合金',
                'cost': 15.0,
                'machinability': 0.3
            },
            'cast_iron_ht250': {
                'E': 138e3,
                'nu': 0.27,
                'rho': 7200,
                'sigma_y': 250,       # 抗拉强度
                'sigma_u': 250,
                'HB': 190,
                'name': 'HT250灰铸铁',
                'cost': 0.8,
                'machinability': 0.8
            }
        }
        
        self.props = self.materials[material_type]
    
    def get_specific_stiffness(self):
        """计算比刚度 (m²/s²)"""
        return self.props['E'] * 1e6 / self.props['rho']
    
    def get_specific_strength(self):
        """计算比强度 (m²/s²)"""
        return self.props['sigma_y'] * 1e6 / self.props['rho']
    
    def get_fatigue_limit(self):
        """估算疲劳极限 (MPa)"""
        # 经验公式：疲劳极限 ≈ 0.5 * 抗拉强度（钢）
        if '钢' in self.props['name'] or 'steel' in self.props['name'].lower():
            return 0.5 * self.props['sigma_u']
        else:
            return 0.4 * self.props['sigma_u']


class MechanicalLoadCase:
    """
    机械零件载荷工况定义
    
    包括静载荷、动载荷、冲击载荷、疲劳载荷等
    """
    
    def __init__(self, load_case_type='static'):
        self.load_case_type = load_case_type
        
        # 定义各种载荷工况
        self.load_cases = {
            'static': {
                'name': '静载荷',
                'description': '恒定静态载荷',
                'load_factor': 1.0,
                'load_type': 'static',
                'magnitude': 1000.0  # N
            },
            'rotating': {
                'name': '旋转载荷',
                'description': '离心力引起的分布载荷',
                'load_factor': 1.2,
                'load_type': 'distributed',
                'omega': 3000,  # rpm
                'magnitude': 500.0  # N
            },
            'impact': {
                'name': '冲击载荷',
                'description': '瞬时冲击载荷',
                'load_factor': 2.0,
                'load_type': 'point',
                'magnitude': 5000.0  # N
            },
            'fatigue': {
                'name': '疲劳载荷',
                'description': '交变循环载荷',
                'load_factor': 1.0,
                'load_type': 'cyclic',
                'sigma_max': 200.0,  # MPa
                'sigma_min': 50.0,   # MPa
                'R': 0.25  # 应力比
            },
            'thermal': {
                'name': '热载荷',
                'description': '温度梯度引起的热应力',
                'load_factor': 1.0,
                'load_type': 'thermal',
                'delta_T': 100  # ℃
            }
        }
        
        self.case = self.load_cases.get(load_case_type, self.load_cases['static'])
    
    def get_load_description(self):
        """获取载荷工况描述"""
        return f"{self.case['name']}: {self.case['description']}"


class SimplifiedMechanicalFEM:
    """
    简化机械零件有限元分析器
    
    用于快速计算机械零件结构响应
    """
    
    def __init__(self, nx, ny, Lx, Ly, material):
        self.nx = nx
        self.ny = ny
        self.Lx = Lx
        self.Ly = Ly
        self.dx = Lx / nx
        self.dy = Ly / ny
        self.material = material
        
        # 节点总数
        self.n_nodes = (nx + 1) * (ny + 1)
        self.n_dof = 2 * self.n_nodes
        
        # 初始化密度场
        self.density = np.ones((ny, nx))
    
    def compute_stiffness_matrix(self):
        """计算全局刚度矩阵（简化版）"""
        K = np.zeros((self.n_dof, self.n_dof))
        
        E0 = self.material.props['E']
        nu = self.material.props['nu']
        
        # 简化刚度计算
        k_base = E0 / (1 - nu**2)
        
        for j in range(self.ny):
            for i in range(self.nx):
                # 单元刚度（简化）
                rho = self.density[j, i]
                k_elem = k_base * rho**3 * np.array([
                    [ 2, -1, -1,  0],
                    [-1,  2,  0, -1],
                    [-1,  0,  2, -1],
                    [ 0, -1, -1,  2]
                ]) * (self.dx * self.dy)
                
                # 组装到全局矩阵
                node_idx = j * (self.nx + 1) + i
                dofs = [
                    2 * node_idx, 2 * node_idx + 1,
                    2 * (node_idx + 1), 2 * (node_idx + 1) + 1,
                    2 * (node_idx + self.nx + 1), 2 * (node_idx + self.nx + 1) + 1,
                    2 * (node_idx + self.nx + 2), 2 * (node_idx + self.nx + 2) + 1
                ]
                
                for ii in range(4):
                    for jj in range(4):
                        if dofs[ii] < self.n_dof and dofs[jj] < self.n_dof:
                            K[dofs[ii], dofs[jj]] += k_elem[ii, jj]
        
        return K
    
    def apply_boundary_conditions(self, K, F, bc_type='fixed_bottom'):
        """应用边界条件"""
        fixed_dofs = []
        
        if bc_type == 'fixed_bottom':
            # 底部固定（轴承座类）
            bottom_nodes = range(0, self.nx + 1)
            for node in bottom_nodes:
                fixed_dofs.extend([2 * node, 2 * node + 1])
                
        elif bc_type == 'fixed_left_right':
            # 左右固定（轴类）
            left_nodes = range(0, (self.ny + 1) * (self.nx + 1), self.nx + 1)
            right_nodes = range(self.nx, (self.ny + 1) * (self.nx + 1), self.nx + 1)
            for node in list(left_nodes) + list(right_nodes):
                fixed_dofs.extend([2 * node, 2 * node + 1])
        
        elif bc_type == 'fixed_hole':
            # 中心孔固定（法兰类）
            center_x, center_y = self.nx // 2, self.ny // 2
            radius = min(self.nx, self.ny) // 4
            for j in range(self.ny + 1):
                for i in range(self.nx + 1):
                    dist = np.sqrt((i - center_x)**2 + (j - center_y)**2)
                    if dist < radius:
                        node_idx = j * (self.nx + 1) + i
                        fixed_dofs.extend([2 * node_idx, 2 * node_idx + 1])
        
        free_dofs = [i for i in range(self.n_dof) if i not in fixed_dofs]
        
        # 移除固定自由度
        K_reduced = K[np.ix_(free_dofs, free_dofs)]
        F_reduced = F[free_dofs]
        
        return K_reduced, F_reduced, free_dofs, fixed_dofs
    
    def solve(self, load_case, bc_type='fixed_bottom'):
        """求解结构响应"""
        K = self.compute_stiffness_matrix()
        
        # 构建载荷向量
        F = np.zeros(self.n_dof)
        
        if load_case.case['load_type'] == 'static':
            # 集中载荷（顶部中心）
            top_center = self.ny * (self.nx + 1) + self.nx // 2
            F[2 * top_center + 1] = -load_case.case['magnitude']
        
        elif load_case.case['load_type'] == 'distributed':
            # 分布载荷
            magnitude = load_case.case['magnitude']
            for j in range(self.ny + 1):
                for i in range(self.nx + 1):
                    node_idx = j * (self.nx + 1) + i
                    F[2 * node_idx + 1] = -magnitude * self.dx * self.dy / 4
        
        elif load_case.case['load_type'] == 'point':
            # 冲击载荷（偏心）
            load_node = (self.ny // 3) * (self.nx + 1) + self.nx // 3
            F[2 * load_node + 1] = -load_case.case['magnitude']
        
        # 应用边界条件
        K_reduced, F_reduced, free_dofs, fixed_dofs = self.apply_boundary_conditions(K, F, bc_type)
        
        # 求解
        try:
            U_reduced = np.linalg.solve(K_reduced, F_reduced)
        except np.linalg.LinAlgError:
            # 添加小量确保正定性
            K_reduced += 1e-6 * np.eye(len(free_dofs))
            U_reduced = np.linalg.solve(K_reduced, F_reduced)
        
        # 重构完整位移向量
        U = np.zeros(self.n_dof)
        U[free_dofs] = U_reduced
        
        # 计算柔度
        compliance = np.dot(F[free_dofs], U_reduced)
        
        return U, compliance


class MechanicalTopologyOptimizer:
    """
    机械零件拓扑优化器
    
    实现多工况、多约束的机械零件拓扑优化
    """
    
    def __init__(self, fem, load_cases, volfrac=0.4, penal=3.0, rmin=1.5):
        self.fem = fem
        self.load_cases = load_cases
        self.volfrac = volfrac
        self.penal = penal
        self.rmin = rmin
        
        # 初始化密度场
        self.density = np.ones((fem.ny, fem.nx)) * volfrac
        
        # 历史记录
        self.history = {
            'compliance': [],
            'volume': [],
            'density_change': [],
            'max_stress': []
        }
    
    def compute_sensitivity(self, load_case, bc_type='fixed_bottom'):
        """计算灵敏度（伴随法）"""
        # 设置当前密度
        self.fem.density = self.density
        
        # 求解位移
        U, compliance = self.fem.solve(load_case, bc_type)
        
        # 计算灵敏度（简化）
        sensitivity = np.zeros_like(self.density)
        
        for j in range(self.fem.ny):
            for i in range(self.fem.nx):
                rho = self.density[j, i]
                # 灵敏度 = -p * rho^(p-1) * 单元应变能
                sensitivity[j, i] = -self.penal * rho**(self.penal - 1) * compliance / (self.fem.nx * self.fem.ny)
        
        return sensitivity, compliance
    
    def filter_sensitivity(self, sensitivity):
        """灵敏度过滤"""
        filtered = np.zeros_like(sensitivity)
        
        for j in range(self.fem.ny):
            for i in range(self.fem.nx):
                sum_weight = 0.0
                weighted_sum = 0.0
                
                for jj in range(max(0, j - int(self.rmin)), min(self.fem.ny, j + int(self.rmin) + 1)):
                    for ii in range(max(0, i - int(self.rmin)), min(self.fem.nx, i + int(self.rmin) + 1)):
                        dist = np.sqrt((i - ii)**2 + (j - jj)**2)
                        if dist < self.rmin:
                            weight = self.rmin - dist
                            weighted_sum += weight * sensitivity[jj, ii]
                            sum_weight += weight
                
                if sum_weight > 0:
                    filtered[j, i] = weighted_sum / sum_weight
        
        return filtered
    
    def update_density_oc(self, sensitivity):
        """使用OC方法更新密度"""
        move_limit = 0.2
        
        # 计算拉格朗日乘子（二分法）
        l1, l2 = 0, 1e6
        while (l2 - l1) / (l1 + l2) > 1e-4:
            lmid = 0.5 * (l1 + l2)
            
            # 更新密度
            new_density = np.zeros_like(self.density)
            for j in range(self.fem.ny):
                for i in range(self.fem.nx):
                    # OC更新公式
                    rho_new = self.density[j, i] * np.sqrt(-sensitivity[j, i] / lmid)
                    
                    # 应用移动限界
                    rho_min = max(0.001, self.density[j, i] - move_limit)
                    rho_max = min(1.0, self.density[j, i] + move_limit)
                    
                    new_density[j, i] = np.clip(rho_new, rho_min, rho_max)
            
            # 检查体积约束
            if np.mean(new_density) > self.volfrac:
                l1 = lmid
            else:
                l2 = lmid
        
        return new_density
    
    def optimize(self, max_iter=50, tol=0.005):
        """执行拓扑优化"""
        print(f"\n开始拓扑优化...")
        print(f"设计域: {self.fem.nx} x {self.fem.ny}")
        print(f"目标体积分数: {self.volfrac}")
        print(f"惩罚因子: {self.penal}")
        print(f"过滤半径: {self.rmin}")
        print("-" * 60)
        
        for iteration in range(max_iter):
            # 多工况灵敏度计算
            total_sensitivity = np.zeros_like(self.density)
            total_compliance = 0
            
            # 工况权重
            weights = [0.5, 0.3, 0.2]  # 静载、旋转、冲击
            
            for load_case, weight in zip(self.load_cases[:3], weights):
                sensitivity, compliance = self.compute_sensitivity(load_case)
                total_sensitivity += weight * sensitivity
                total_compliance += weight * compliance
            
            # 灵敏度过滤
            filtered_sensitivity = self.filter_sensitivity(total_sensitivity)
            
            # 更新密度
            old_density = self.density.copy()
            self.density = self.update_density_oc(filtered_sensitivity)
            
            # 计算变化
            density_change = np.max(np.abs(self.density - old_density))
            volume = np.mean(self.density)
            
            # 记录历史
            self.history['compliance'].append(total_compliance)
            self.history['volume'].append(volume)
            self.history['density_change'].append(density_change)
            
            # 输出进度
            if iteration % 5 == 0 or density_change < tol:
                print(f"迭代 {iteration:3d}: 柔度 = {total_compliance:.4e}, "
                      f"体积 = {volume:.4f}, 变化 = {density_change:.6f}")
            
            # 检查收敛
            if density_change < tol and iteration > 10:
                print(f"\n优化收敛于迭代 {iteration}")
                break
        
        print("-" * 60)
        print(f"优化完成!")
        print(f"最终体积分数: {np.mean(self.density):.4f}")
        print(f"最终柔度: {self.history['compliance'][-1]:.4e}")
        
        return self.density, self.history


def create_mechanical_parts_schematic():
    """创建机械零件示意图"""
    fig, axes = plt.subplots(2, 2, figsize=(14, 12))
    
    # 1. 轴承座
    ax1 = axes[0, 0]
    ax1.set_xlim(0, 10)
    ax1.set_ylim(0, 10)
    ax1.set_aspect('equal')
    ax1.axis('off')
    ax1.set_title('轴承座', fontsize=12, fontweight='bold')
    
    # 底座
    ax1.add_patch(Rectangle((1, 1), 8, 2, facecolor='lightgray', edgecolor='black', linewidth=2))
    # 立板
    ax1.add_patch(Rectangle((3, 3), 4, 5, facecolor='lightblue', edgecolor='blue', linewidth=2))
    # 轴承孔
    circle = Circle((5, 7), 1.2, facecolor='white', edgecolor='red', linewidth=2)
    ax1.add_patch(circle)
    # 螺栓孔
    for x, y in [(2, 1.5), (8, 1.5), (2, 8.5), (8, 8.5)]:
        ax1.add_patch(Circle((x, y), 0.3, facecolor='white', edgecolor='black'))
    
    # 载荷
    ax1.annotate('', xy=(5, 5), xytext=(5, 3),
                arrowprops=dict(arrowstyle='->', color='red', lw=2))
    ax1.text(5, 2.5, '轴承载荷', ha='center', fontsize=10, color='red')
    
    # 2. 齿轮
    ax2 = axes[0, 1]
    ax2.set_xlim(0, 10)
    ax2.set_ylim(0, 10)
    ax2.set_aspect('equal')
    ax2.axis('off')
    ax2.set_title('齿轮', fontsize=12, fontweight='bold')
    
    # 齿顶圆
    outer_circle = Circle((5, 5), 4, facecolor='lightgreen', edgecolor='green', linewidth=2)
    ax2.add_patch(outer_circle)
    # 齿根圆
    inner_circle = Circle((5, 5), 3, facecolor='white', edgecolor='green', linewidth=2)
    ax2.add_patch(inner_circle)
    # 轮毂
    hub_circle = Circle((5, 5), 1.5, facecolor='lightyellow', edgecolor='orange', linewidth=2)
    ax2.add_patch(hub_circle)
    # 轴孔
    shaft_circle = Circle((5, 5), 0.8, facecolor='white', edgecolor='black')
    ax2.add_patch(shaft_circle)
    
    # 轮齿示意
    for angle in np.linspace(0, 2*np.pi, 12, endpoint=False):
        x_outer = 5 + 4 * np.cos(angle)
        y_outer = 5 + 4 * np.sin(angle)
        x_inner = 5 + 3 * np.cos(angle)
        y_inner = 5 + 3 * np.sin(angle)
        ax2.plot([x_inner, x_outer], [y_inner, y_outer], 'g-', linewidth=2)
    
    # 3. 联轴器
    ax3 = axes[1, 0]
    ax3.set_xlim(0, 12)
    ax3.set_ylim(0, 6)
    ax3.set_aspect('equal')
    ax3.axis('off')
    ax3.set_title('联轴器', fontsize=12, fontweight='bold')
    
    # 左半联轴器
    ax3.add_patch(Rectangle((1, 2), 3, 2, facecolor='lightcoral', edgecolor='red', linewidth=2))
    # 右半联轴器
    ax3.add_patch(Rectangle((8, 2), 3, 2, facecolor='lightcoral', edgecolor='red', linewidth=2))
    # 连接螺栓
    for x in [4.5, 7.5]:
        ax3.add_patch(Circle((x, 3), 0.3, facecolor='gray', edgecolor='black'))
    # 轴孔
    ax3.add_patch(Circle((2.5, 3), 0.8, facecolor='white', edgecolor='black'))
    ax3.add_patch(Circle((9.5, 3), 0.8, facecolor='white', edgecolor='black'))
    
    # 4. 法兰
    ax4 = axes[1, 1]
    ax4.set_xlim(0, 10)
    ax4.set_ylim(0, 10)
    ax4.set_aspect('equal')
    ax4.axis('off')
    ax4.set_title('法兰', fontsize=12, fontweight='bold')
    
    # 法兰盘
    outer_circle = Circle((5, 5), 4, facecolor='lightcyan', edgecolor='cyan', linewidth=2)
    ax4.add_patch(outer_circle)
    # 中心孔
    center_circle = Circle((5, 5), 1.5, facecolor='white', edgecolor='black', linewidth=2)
    ax4.add_patch(center_circle)
    # 螺栓孔
    for angle in np.linspace(0, 2*np.pi, 8, endpoint=False):
        x = 5 + 2.5 * np.cos(angle)
        y = 5 + 2.5 * np.sin(angle)
        ax4.add_patch(Circle((x, y), 0.4, facecolor='white', edgecolor='black'))
    
    plt.tight_layout()
    plt.savefig('mechanical_parts_schematic.png', dpi=150, bbox_inches='tight')
    plt.close()
    print("✓ 机械零件示意图已保存")


def compare_mechanical_materials():
    """对比机械工程材料性能"""
    materials = ['steel_45', 'steel_40cr', 'aluminum_6061', 'titanium_tc4', 'cast_iron_ht250']
    
    results = []
    for mat_type in materials:
        mat = MechanicalMaterial(mat_type)
        results.append({
            'name': mat.props['name'],
            'E': mat.props['E'] / 1e3,  # GPa
            'rho': mat.props['rho'],
            'sigma_y': mat.props['sigma_y'],
            'HB': mat.props['HB'],
            'specific_stiffness': mat.get_specific_stiffness() / 1e6,
            'specific_strength': mat.get_specific_strength() / 1e6,
            'fatigue_limit': mat.get_fatigue_limit(),
            'cost': mat.props['cost'],
            'machinability': mat.props['machinability']
        })
    
    # 可视化对比
    fig, axes = plt.subplots(3, 3, figsize=(16, 14))
    
    names = [r['name'] for r in results]
    x = np.arange(len(names))
    
    # 弹性模量
    ax1 = axes[0, 0]
    E_values = [r['E'] for r in results]
    bars1 = ax1.bar(x, E_values, color=['#3498db', '#2980b9', '#1abc9c', '#e74c3c', '#9b59b6'])
    ax1.set_ylabel('弹性模量 (GPa)', fontsize=11)
    ax1.set_title('弹性模量对比', fontsize=12, fontweight='bold')
    ax1.set_xticks(x)
    ax1.set_xticklabels(names, rotation=45, ha='right', fontsize=9)
    ax1.grid(axis='y', alpha=0.3)
    
    # 密度
    ax2 = axes[0, 1]
    rho_values = [r['rho'] for r in results]
    bars2 = ax2.bar(x, rho_values, color=['#3498db', '#2980b9', '#1abc9c', '#e74c3c', '#9b59b6'])
    ax2.set_ylabel('密度 (kg/m³)', fontsize=11)
    ax2.set_title('密度对比', fontsize=12, fontweight='bold')
    ax2.set_xticks(x)
    ax2.set_xticklabels(names, rotation=45, ha='right', fontsize=9)
    ax2.grid(axis='y', alpha=0.3)
    
    # 屈服强度
    ax3 = axes[0, 2]
    sigma_values = [r['sigma_y'] for r in results]
    bars3 = ax3.bar(x, sigma_values, color=['#3498db', '#2980b9', '#1abc9c', '#e74c3c', '#9b59b6'])
    ax3.set_ylabel('屈服强度 (MPa)', fontsize=11)
    ax3.set_title('屈服强度对比', fontsize=12, fontweight='bold')
    ax3.set_xticks(x)
    ax3.set_xticklabels(names, rotation=45, ha='right', fontsize=9)
    ax3.grid(axis='y', alpha=0.3)
    
    # 比刚度
    ax4 = axes[1, 0]
    ss_values = [r['specific_stiffness'] for r in results]
    bars4 = ax4.bar(x, ss_values, color=['#3498db', '#2980b9', '#1abc9c', '#e74c3c', '#9b59b6'])
    ax4.set_ylabel('比刚度 (×10⁶ m²/s²)', fontsize=11)
    ax4.set_title('比刚度对比', fontsize=12, fontweight='bold')
    ax4.set_xticks(x)
    ax4.set_xticklabels(names, rotation=45, ha='right', fontsize=9)
    ax4.grid(axis='y', alpha=0.3)
    
    # 比强度
    ax5 = axes[1, 1]
    sst_values = [r['specific_strength'] for r in results]
    bars5 = ax5.bar(x, sst_values, color=['#3498db', '#2980b9', '#1abc9c', '#e74c3c', '#9b59b6'])
    ax5.set_ylabel('比强度 (×10⁶ m²/s²)', fontsize=11)
    ax5.set_title('比强度对比', fontsize=12, fontweight='bold')
    ax5.set_xticks(x)
    ax5.set_xticklabels(names, rotation=45, ha='right', fontsize=9)
    ax5.grid(axis='y', alpha=0.3)
    
    # 疲劳极限
    ax6 = axes[1, 2]
    fatigue_values = [r['fatigue_limit'] for r in results]
    bars6 = ax6.bar(x, fatigue_values, color=['#3498db', '#2980b9', '#1abc9c', '#e74c3c', '#9b59b6'])
    ax6.set_ylabel('疲劳极限 (MPa)', fontsize=11)
    ax6.set_title('疲劳极限对比', fontsize=12, fontweight='bold')
    ax6.set_xticks(x)
    ax6.set_xticklabels(names, rotation=45, ha='right', fontsize=9)
    ax6.grid(axis='y', alpha=0.3)
    
    # 硬度
    ax7 = axes[2, 0]
    hb_values = [r['HB'] for r in results]
    bars7 = ax7.bar(x, hb_values, color=['#3498db', '#2980b9', '#1abc9c', '#e74c3c', '#9b59b6'])
    ax7.set_ylabel('布氏硬度 HB', fontsize=11)
    ax7.set_title('硬度对比', fontsize=12, fontweight='bold')
    ax7.set_xticks(x)
    ax7.set_xticklabels(names, rotation=45, ha='right', fontsize=9)
    ax7.grid(axis='y', alpha=0.3)
    
    # 成本
    ax8 = axes[2, 1]
    cost_values = [r['cost'] for r in results]
    bars8 = ax8.bar(x, cost_values, color=['#3498db', '#2980b9', '#1abc9c', '#e74c3c', '#9b59b6'])
    ax8.set_ylabel('相对成本', fontsize=11)
    ax8.set_title('成本对比', fontsize=12, fontweight='bold')
    ax8.set_xticks(x)
    ax8.set_xticklabels(names, rotation=45, ha='right', fontsize=9)
    ax8.grid(axis='y', alpha=0.3)
    
    # 可加工性
    ax9 = axes[2, 2]
    mach_values = [r['machinability'] for r in results]
    bars9 = ax9.bar(x, mach_values, color=['#3498db', '#2980b9', '#1abc9c', '#e74c3c', '#9b59b6'])
    ax9.set_ylabel('可加工性指数', fontsize=11)
    ax9.set_title('可加工性对比', fontsize=12, fontweight='bold')
    ax9.set_xticks(x)
    ax9.set_xticklabels(names, rotation=45, ha='right', fontsize=9)
    ax9.grid(axis='y', alpha=0.3)
    ax9.set_ylim(0, 1.2)
    
    plt.tight_layout()
    plt.savefig('mechanical_material_comparison.png', dpi=150, bbox_inches='tight')
    plt.close()
    print("✓ 机械材料对比图已保存")
    
    return results


def visualize_optimization_result(density, fem, filename='mechanical_optimization_result.png'):
    """可视化拓扑优化结果"""
    fig, axes = plt.subplots(1, 3, figsize=(18, 6))
    
    # 1. 密度分布热力图
    ax1 = axes[0]
    im1 = ax1.imshow(density, cmap='RdYlBu_r', origin='lower', aspect='auto',
                     extent=[0, fem.Lx, 0, fem.Ly])
    ax1.set_xlabel('长度 (mm)', fontsize=11)
    ax1.set_ylabel('高度 (mm)', fontsize=11)
    ax1.set_title('密度分布热力图', fontsize=12, fontweight='bold')
    cbar1 = plt.colorbar(im1, ax=ax1)
    cbar1.set_label('密度', fontsize=10)
    
    # 2. 结构轮廓
    ax2 = axes[1]
    contour = ax2.contour(density, levels=[0.5], colors='blue', linewidths=2,
                          extent=[0, fem.Lx, 0, fem.Ly])
    ax2.contourf(density, levels=[0.5, 1.0], colors=['lightblue'], alpha=0.5,
                 extent=[0, fem.Lx, 0, fem.Ly])
    ax2.set_xlabel('长度 (mm)', fontsize=11)
    ax2.set_ylabel('高度 (mm)', fontsize=11)
    ax2.set_title('优化结构轮廓', fontsize=12, fontweight='bold')
    ax2.set_aspect('equal')
    ax2.grid(True, alpha=0.3)
    
    # 3. 二值化结构
    ax3 = axes[2]
    binary = (density > 0.5).astype(float)
    im3 = ax3.imshow(binary, cmap='binary', origin='lower', aspect='auto',
                     extent=[0, fem.Lx, 0, fem.Ly])
    ax3.set_xlabel('长度 (mm)', fontsize=11)
    ax3.set_ylabel('高度 (mm)', fontsize=11)
    ax3.set_title('二值化结构', fontsize=12, fontweight='bold')
    
    plt.tight_layout()
    plt.savefig(filename, dpi=150, bbox_inches='tight')
    plt.close()
    print(f"✓ 优化结果图已保存: {filename}")


def plot_optimization_history(history, filename='mechanical_optimization_history.png'):
    """绘制优化历史"""
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    iterations = range(len(history['compliance']))
    
    # 1. 柔度收敛
    ax1 = axes[0, 0]
    ax1.semilogy(iterations, history['compliance'], 'b-', linewidth=2, marker='o', markersize=4)
    ax1.set_xlabel('迭代次数', fontsize=11)
    ax1.set_ylabel('结构柔度 (对数)', fontsize=11)
    ax1.set_title('柔度收敛曲线', fontsize=12, fontweight='bold')
    ax1.grid(True, alpha=0.3)
    
    # 2. 体积分数收敛
    ax2 = axes[0, 1]
    ax2.plot(iterations, history['volume'], 'r-', linewidth=2, marker='s', markersize=4)
    ax2.axhline(y=0.4, color='k', linestyle='--', label='目标体积分数')
    ax2.set_xlabel('迭代次数', fontsize=11)
    ax2.set_ylabel('体积分数', fontsize=11)
    ax2.set_title('体积分数收敛曲线', fontsize=12, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # 3. 密度变化
    ax3 = axes[1, 0]
    ax3.semilogy(iterations, history['density_change'], 'g-', linewidth=2, marker='^', markersize=4)
    ax3.axhline(y=0.005, color='k', linestyle='--', label='收敛容差')
    ax3.set_xlabel('迭代次数', fontsize=11)
    ax3.set_ylabel('最大密度变化 (对数)', fontsize=11)
    ax3.set_title('密度变化曲线', fontsize=12, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # 4. 综合收敛
    ax4 = axes[1, 1]
    ax4.plot(iterations, np.array(history['compliance']) / max(history['compliance']), 
             'b-', linewidth=2, label='归一化柔度')
    ax4.plot(iterations, history['volume'], 'r-', linewidth=2, label='体积分数')
    ax4.set_xlabel('迭代次数', fontsize=11)
    ax4.set_ylabel('归一化值', fontsize=11)
    ax4.set_title('综合收敛曲线', fontsize=12, fontweight='bold')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig(filename, dpi=150, bbox_inches='tight')
    plt.close()
    print(f"✓ 优化历史图已保存: {filename}")


def create_optimization_animation(density_history, fem, filename='mechanical_optimization.gif'):
    """创建优化过程动画"""
    frames = []
    
    # 选择关键帧
    n_frames = min(len(density_history), 20)
    indices = np.linspace(0, len(density_history) - 1, n_frames, dtype=int)
    
    for idx in indices:
        density = density_history[idx]
        
        fig, ax = plt.subplots(figsize=(10, 8))
        im = ax.imshow(density, cmap='RdYlBu_r', origin='lower', aspect='auto',
                       vmin=0, vmax=1, extent=[0, fem.Lx, 0, fem.Ly])
        ax.set_xlabel('长度 (mm)', fontsize=11)
        ax.set_ylabel('高度 (mm)', fontsize=11)
        ax.set_title(f'机械零件拓扑优化 - 迭代 {idx}', fontsize=12, fontweight='bold')
        
        # 添加颜色条
        cbar = plt.colorbar(im, ax=ax)
        cbar.set_label('密度', fontsize=10)
        
        plt.tight_layout()
        
        # 保存到内存
        buf = io.BytesIO()
        plt.savefig(buf, format='png', dpi=100, bbox_inches='tight')
        buf.seek(0)
        frames.append(Image.open(buf))
        plt.close()
    
    # 保存GIF
    frames[0].save(
        filename,
        save_all=True,
        append_images=frames[1:],
        duration=300,
        loop=0
    )
    print(f"✓ 优化动画已保存: {filename}")


def mechanical_parts_optimization():
    """
    机械零件拓扑优化主函数
    
    演示机械零件的轻量化设计
    """
    
    print("=" * 80)
    print("机械零件优化 - Mechanical Parts Optimization")
    print("=" * 80)
    print("\n本案例演示机械零件的拓扑优化设计")
    print("考虑多工况：静载、旋转、冲击")
    print("=" * 80)
    
    # 1. 材料选择与对比
    print("\n【步骤1】机械材料性能对比分析")
    material_results = compare_mechanical_materials()
    
    # 选择优化材料（40Cr合金钢）
    material_type = 'steel_40cr'
    material = MechanicalMaterial(material_type)
    
    print(f"\n选择优化材料: {material.props['name']}")
    print(f"  弹性模量: {material.props['E']/1e3:.1f} GPa")
    print(f"  密度: {material.props['rho']:.0f} kg/m³")
    print(f"  屈服强度: {material.props['sigma_y']:.1f} MPa")
    print(f"  疲劳极限: {material.get_fatigue_limit():.1f} MPa")
    
    # 2. 创建机械零件示意图
    print("\n【步骤2】生成机械零件示意图")
    create_mechanical_parts_schematic()
    
    # 3. 设置优化问题（轴承座）
    print("\n【步骤3】设置轴承座拓扑优化问题")
    
    # 设计域参数（轴承座）
    nx, ny = 30, 30
    Lx, Ly = 150.0, 150.0  # mm
    
    # 创建FEM模型
    fem = SimplifiedMechanicalFEM(nx, ny, Lx, Ly, material)
    
    # 定义载荷工况
    load_cases = [
        MechanicalLoadCase('static'),
        MechanicalLoadCase('rotating'),
        MechanicalLoadCase('impact')
    ]
    
    print(f"\n载荷工况:")
    for i, lc in enumerate(load_cases, 1):
        print(f"  {i}. {lc.get_load_description()}")
    
    # 4. 执行拓扑优化
    print("\n【步骤4】执行轴承座拓扑优化")
    
    optimizer = MechanicalTopologyOptimizer(
        fem=fem,
        load_cases=load_cases,
        volfrac=0.40,  # 体积分数40%
        penal=3.0,
        rmin=2.0
    )
    
    density_opt, history = optimizer.optimize(max_iter=50, tol=0.005)
    
    # 5. 可视化结果
    print("\n【步骤5】生成可视化结果")
    
    # 优化结果
    visualize_optimization_result(density_opt, fem, 'mechanical_optimization_result.png')
    
    # 优化历史
    plot_optimization_history(history, 'mechanical_optimization_history.png')
    
    # 优化动画
    density_history = [np.ones_like(density_opt) * 0.4]  # 初始
    for i in range(1, len(history['compliance'])):
        # 模拟密度演化
        ratio = i / len(history['compliance'])
        temp_density = density_opt * ratio + 0.4 * (1 - ratio)
        density_history.append(temp_density)
    density_history.append(density_opt)
    
    create_optimization_animation(density_history, fem, 'mechanical_optimization.gif')
    
    # 6. 输出优化结果
    print("\n【步骤6】优化结果汇总")
    print("=" * 60)
    
    # 计算结构重量
    total_volume = Lx * Ly * 20  # 假设厚度20mm
    used_volume = np.mean(density_opt) * total_volume
    weight = used_volume * material.props['rho'] * 1e-9  # kg
    
    print(f"材料: {material.props['name']}")
    print(f"最终体积分数: {np.mean(density_opt):.4f}")
    print(f"结构重量: {weight:.3f} kg")
    print(f"最终柔度: {history['compliance'][-1]:.4e}")
    print(f"迭代次数: {len(history['compliance'])}")
    print("=" * 60)
    
    # 7. 材料对比优化
    print("\n【步骤7】不同材料的优化结果对比")
    
    fig, axes = plt.subplots(2, 3, figsize=(15, 10))
    materials_to_compare = ['steel_45', 'aluminum_6061', 'titanium_tc4']
    
    for idx, mat_type in enumerate(materials_to_compare):
        mat = MechanicalMaterial(mat_type)
        fem_mat = SimplifiedMechanicalFEM(nx, ny, Lx, Ly, mat)
        
        optimizer_mat = MechanicalTopologyOptimizer(
            fem=fem_mat,
            load_cases=load_cases,
            volfrac=0.40,
            penal=3.0,
            rmin=2.0
        )
        
        density_mat, _ = optimizer_mat.optimize(max_iter=30, tol=0.01)
        
        # 绘制结果
        row = idx // 3
        col = idx % 3
        ax = axes[row, col]
        
        im = ax.imshow(density_mat, cmap='RdYlBu_r', origin='lower', aspect='auto',
                       extent=[0, Lx, 0, Ly])
        ax.set_title(f'{mat.props["name"]}', fontsize=11, fontweight='bold')
        ax.set_xlabel('长度 (mm)', fontsize=10)
        ax.set_ylabel('高度 (mm)', fontsize=10)
        
        # 计算重量
        weight_mat = np.mean(density_mat) * total_volume * mat.props['rho'] * 1e-9
        ax.text(10, 140, f'重量: {weight_mat:.3f} kg', fontsize=9, 
                bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))
    
    # 隐藏多余的子图
    axes[1, 1].axis('off')
    axes[1, 2].axis('off')
    
    plt.tight_layout()
    plt.savefig('mechanical_material_optimization_comparison.png', dpi=150, bbox_inches='tight')
    plt.close()
    print("✓ 材料对比优化结果已保存")
    
    print("\n" + "=" * 80)
    print("机械零件优化完成!")
    print("=" * 80)
    
    return density_opt, history, material


if __name__ == "__main__":
    # 运行机械零件优化
    density, history, material = mechanical_parts_optimization()

---

四、代码深度解析

4.1 机械材料模型详解

MechanicalMaterial类实现了机械工程常用材料的属性定义：

def __init__(self, material_type='steel_45'):
    self.material_type = material_type
    
    # 机械工程常用材料属性定义
    self.materials = {
        'steel_45': {
            'E': 210e3,           # 杨氏模量 (MPa)
            'nu': 0.30,           # 泊松比
            'rho': 7850,          # 密度 (kg/m³)
            'sigma_y': 355,       # 屈服强度 (MPa)
            'sigma_u': 600,       # 抗拉强度 (MPa)
            'HB': 197,            # 布氏硬度
            'name': '45号钢',
            'cost': 1.0,          # 相对成本
            'machinability': 0.7  # 可加工性指数
        },
        # ... 其他材料定义
    }

关键参数说明：

• 弹性模量E：材料抵抗弹性变形的能力，钢材约210 GPa
• 屈服强度σ_y：材料开始塑性变形的应力，40Cr调质后可达785 MPa
• 疲劳极限：材料在无限次循环下不破坏的最大应力，钢材约为抗拉强度的一半
• 可加工性指数：0-1之间，越高表示越容易加工

4.2 载荷工况定义

MechanicalLoadCase类定义了机械零件常见的载荷类型：

def __init__(self, load_case_type='static'):
    self.load_cases = {
        'static': {
            'name': '静载荷',
            'description': '恒定静态载荷',
            'load_factor': 1.0,
            'load_type': 'static',
            'magnitude': 1000.0  # N
        },
        'rotating': {
            'name': '旋转载荷',
            'description': '离心力引起的分布载荷',
            'load_factor': 1.2,
            'load_type': 'distributed',
            'omega': 3000,  # rpm
            'magnitude': 500.0  # N
        },
        # ... 其他载荷类型
    }

载荷类型说明：

• 静载荷：恒定不变的载荷，如自重、预紧力
• 旋转载荷：由旋转产生的离心力，如齿轮、轴
• 冲击载荷：瞬时作用的载荷，如启动、制动
• 疲劳载荷：周期性变化的载荷，需进行疲劳分析

4.3 拓扑优化核心算法

MechanicalTopologyOptimizer类实现了SIMP拓扑优化算法：

def compute_sensitivity(self, load_case, bc_type='fixed_bottom'):
    """计算灵敏度（伴随法）"""
    # 设置当前密度
    self.fem.density = self.density
    
    # 求解位移
    U, compliance = self.fem.solve(load_case, bc_type)
    
    # 计算灵敏度
    sensitivity = np.zeros_like(self.density)
    for j in range(self.fem.ny):
        for i in range(self.fem.nx):
            rho = self.density[j, i]
            # 灵敏度 = -p * rho^(p-1) * 单元应变能
            sensitivity[j, i] = -self.penal * rho**(self.penal - 1) * compliance / (self.fem.nx * self.fem.ny)
    
    return sensitivity, compliance

算法步骤：

1. 有限元分析：求解当前密度分布下的位移场和柔度
2. 灵敏度计算：使用伴随法计算目标函数对设计变量的梯度
3. 灵敏度过滤：消除数值不稳定性和棋盘格现象
4. 密度更新：使用OC方法更新密度分布
5. 收敛判断：检查密度变化是否满足收敛准则

4.4 边界条件处理

def apply_boundary_conditions(self, K, F, bc_type='fixed_bottom'):
    """应用边界条件"""
    fixed_dofs = []
    
    if bc_type == 'fixed_bottom':
        # 底部固定（轴承座类）
        bottom_nodes = range(0, self.nx + 1)
        for node in bottom_nodes:
            fixed_dofs.extend([2 * node, 2 * node + 1])
    
    elif bc_type == 'fixed_left_right':
        # 左右固定（轴类）
        left_nodes = range(0, (self.ny + 1) * (self.nx + 1), self.nx + 1)
        right_nodes = range(self.nx, (self.ny + 1) * (self.nx + 1), self.nx + 1)
        for node in list(left_nodes) + list(right_nodes):
            fixed_dofs.extend([2 * node, 2 * node + 1])

边界条件类型：

• 底部固定：适用于轴承座、支架类零件
• 左右固定：适用于轴类零件
• 中心孔固定：适用于法兰、轮毂类零件

---

五、运行结果预期

5.1 材料性能对比

运行代码后将生成mechanical_material_comparison.png，展示五种机械工程材料的性能对比：

关键发现：

• 45号钢：综合性能均衡，成本低廉，是最常用的结构材料
• 40Cr合金钢：调质后强度极高，适合重载零件
• 6061铝合金：密度仅为钢的1/3，比刚度优异，适合轻量化设计
• TC4钛合金：比强度最高，但成本昂贵，主要用于航空航天
• HT250灰铸铁：成本低，铸造性能好，适合复杂形状零件

5.2 优化收敛曲线

mechanical_optimization_history.png展示优化过程的收敛特性：

预期特征：

• 柔度在初始阶段快速下降，随后趋于平稳
• 体积分数逐渐收敛到目标值40%
• 密度变化在10-15次迭代后达到收敛容差

5.3 优化结构形态

mechanical_optimization_result.png展示最终的拓扑优化结果：

结构特征：

• 材料主要集中在载荷传递路径上
• 高应力区域密度接近1.0（实体材料）
• 低应力区域密度接近0.001（空腔）
• 边界呈现清晰的0-1分布

---