多场耦合优化-主题033-拓扑优化方法

kkchenkx

373人浏览 · 2026-03-16 22:21:52

kkchenkx · 2026-03-16 22:21:52 发布

主题033：拓扑优化方法

一、引言

拓扑优化（Topology Optimization）是结构优化设计中的最高层次，它能够在给定的设计空间内自动寻找最优的材料分布，从而获得性能最优的结构形式。与尺寸优化（改变截面尺寸）和形状优化（改变边界形状）不同，拓扑优化可以改变结构的连通性和拓扑形式，具有更大的设计自由度。

拓扑优化方法在航空航天、汽车制造、建筑结构、微机电系统（MEMS）等领域有着广泛的应用。例如，飞机机翼的内部加强筋布局、汽车车身的轻量化设计、桥梁的结构形式等，都可以通过拓扑优化获得创新性的设计方案。
在这里插入图片描述

1.1 拓扑优化的发展历程

拓扑优化的发展可以追溯到20世纪80年代。1988年，Bendsøe和Kikuchi提出了基于均匀化方法的拓扑优化理论，奠定了现代拓扑优化的数学基础。随后，各种拓扑优化方法相继涌现：

1988年：均匀化方法（Homogenization Method）
1989年：密度法（SIMP方法，Solid Isotropic Material with Penalization）
1993年：进化结构优化方法（ESO，Evolutionary Structural Optimization）
2000年：水平集方法（Level Set Method）
2010年后：基于机器学习的拓扑优化方法

1.2 拓扑优化的数学描述

拓扑优化问题可以表述为：在给定的设计域Ω内，寻找材料的最优分布，使得结构的某个性能指标最优，同时满足一定的约束条件。

一般形式：

最小化:  J(ρ) = F(u(ρ))
约束条件:
  ∫_Ω ρ dΩ ≤ V*  (体积约束)
  0 ≤ ρ ≤ 1      (密度约束)
  平衡方程: ∇·σ + f = 0

其中：

ρ：材料密度（设计变量）
J：目标函数（如柔度、应力、频率等）
V*：允许的最大材料体积
u：位移场
σ：应力张量
f：体积力

1.3 拓扑优化的基本流程

拓扑优化的基本流程包括以下步骤：

定义设计域：确定可设计区域和非设计区域
离散化：将设计域划分为有限元网格
初始化：设置初始材料分布（通常均匀分布）
有限元分析：计算结构的响应（位移、应力等）
灵敏度分析：计算目标函数和约束对设计变量的灵敏度
设计变量更新：根据优化准则更新材料密度
收敛判断：检查是否满足收敛条件，否则返回步骤4

二、SIMP方法（固体各向同性材料惩罚法）

2.1 SIMP方法的基本原理

SIMP（Solid Isotropic Material with Penalization）方法是目前应用最广泛的拓扑优化方法。其核心思想是引入一个假想的密度变量ρ∈[0,1]，将材料的弹性模量与密度关联：

材料插值模型：

E(ρ) = E_min + ρ^p (E_0 - E_min)

其中：

E_0：实体材料的弹性模量
E_min：空材料的最小弹性模量（避免奇异，通常取E_0×10^-9）
p：惩罚因子（通常取p≥3）
ρ：单元密度（0≤ρ≤1）

惩罚因子的作用：

惩罚因子p的作用是使中间密度（0<ρ<1）的材料变得不经济。当p>1时，中间密度的材料具有较低的刚度-重量比，优化过程会倾向于将单元推向0或1，从而获得清晰的0-1分布。

例如，当p=3时：

ρ=0.5时，E/E_0 = 0.5^3 = 0.125，即50%的材料只提供12.5%的刚度
这意味着使用中间密度材料是低效的，优化会避免这种情况

2.2 SIMP方法的数学模型

最小柔度问题（刚度最大化）：

最小化:  C(ρ) = F^T U = U^T K U
约束条件:
  ∑(ρ_e · v_e) ≤ V*    (体积约束)
  0 ≤ ρ_e ≤ 1, ∀e      (密度约束)
  KU = F               (平衡方程)

其中：

C：柔度（应变能，与刚度成反比）
F：载荷向量
U：位移向量
K：整体刚度矩阵
ρ_e：第e个单元的密度
v_e：第e个单元的体积

单元刚度矩阵：

K_e = ∫_Ωe B^T D(ρ_e) B dΩ

其中B是应变-位移矩阵，D是弹性矩阵：

D(ρ) = [E(ρ)/(1-ν²)] × [1   ν   0
                         ν   1   0
                         0   0   (1-ν)/2]

2.3 灵敏度分析

灵敏度分析是拓扑优化的核心，用于计算目标函数对设计变量的导数。

柔度对密度的灵敏度：

∂C/∂ρ_e = -U_e^T (∂K_e/∂ρ_e) U_e = -p·ρ_e^(p-1)·U_e^T K_e^0 U_e

其中K_e^0是实体单元的刚度矩阵。

物理意义：

灵敏度表示增加单位材料对柔度的影响
负的灵敏度意味着增加材料会降低柔度（提高刚度）
优化过程会优先在灵敏度大的区域增加材料

2.4 优化准则法（OC方法）

优化准则法（Optimality Criteria Method）是SIMP方法中常用的更新策略。

更新公式：

ρ_e^(new) = 
  max(ρ_min, ρ_e - m)   如果 ρ_e·B_e^η ≤ max(ρ_min, ρ_e - m)
  min(1, ρ_e + m)       如果 ρ_e·B_e^η ≥ min(1, ρ_e + m)
  ρ_e·B_e^η             其他情况

其中：

B_e = (-∂C/∂ρ_e) / (λ·v_e)

m：移动限（通常取0.2）
η：阻尼系数（通常取0.5）
λ：拉格朗日乘子（通过二分法确定以满足体积约束）

算法流程：

初始化设计变量ρ
有限元分析，计算位移U
计算灵敏度∂C/∂ρ
使用OC方法更新设计变量
检查收敛条件
如果未收敛，返回步骤2

2.5 数值稳定性技术

拓扑优化中存在一些数值问题，需要采用特殊技术处理：

（1）棋盘格现象

棋盘格现象是指优化结果中出现材料和非材料单元交替排列的棋盘状图案。这是由于有限元离散化的数值误差造成的。

解决方法：

密度过滤：使用加权平均方法平滑密度场
```
ρ̃_e = ∑(w_ei · ρ_i) / ∑w_ei
```
其中w_ei是基于距离的权重函数。

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

∂C̃/∂ρ_e = ∑(w_ei · ρ_i · ∂C/∂ρ_i) / (ρ_e · ∑w_ei)

（2）网格依赖性

优化结果会随着网格细化而变化，细网格会产生更细小的结构特征。

解决方法：

设置最小特征尺寸限制
使用投影方法（Projection Method）
采用稳健性拓扑优化

（3）灰度单元

灰度单元是指密度介于0和1之间的单元，这会造成制造困难。

解决方法：

增大惩罚因子p
使用Heaviside投影
采用形态学过滤

三、其他拓扑优化方法

3.1 进化结构优化方法（ESO/BESO）

ESO（Evolutionary Structural Optimization）方法基于一个简单的概念：逐渐从结构中删除低效的材料，从而进化出最优结构。

基本思想：

计算每个单元的灵敏度（如应变能密度）
删除灵敏度最低的单元（材料利用率低）
重复直到满足体积约束

BESO（双向ESO）改进：

不仅可以删除材料，还可以在高效区域添加材料
实现了材料的重新分布，而不仅是删除

灵敏度准则：

α_e = (1/2) u_e^T K_e u_e / v_e   (应变能密度)

算法流程：

初始化完整结构
有限元分析
计算所有单元的灵敏度
按灵敏度排序
删除灵敏度最低的x%单元
检查收敛，未收敛则返回步骤2

3.2 水平集方法

水平集方法（Level Set Method）将结构的边界表示为一个隐函数的零水平集，通过演化这个隐函数来改变结构形状和拓扑。

水平集函数：

φ(x,t) = 0    在边界上
φ(x,t) > 0    在材料内部
φ(x,t) < 0    在材料外部

Hamilton-Jacobi方程：

∂φ/∂t + V_n |∇φ| = 0

其中V_n是边界的法向速度，通常与形状导数相关。

优点：

边界清晰，没有灰度单元
便于控制结构的几何特征
适合处理复杂拓扑变化

缺点：

计算复杂度较高
可能出现数值不稳定
需要重新初始化水平集函数

3.3 基于密度法的改进方法

（1）稳健性拓扑优化

考虑制造误差和载荷不确定性，设计对扰动不敏感的结构。

最小化:  max_δ C(ρ, δ)
约束条件: 体积约束

其中δ表示不确定性参数。

（2）多材料拓扑优化

同时使用多种材料，优化它们的分布。

E = ∑(ρ_i^p · E_i)

（3）多尺度拓扑优化

考虑微观结构的影响，设计具有特定等效性能的材料分布。

四、拓扑优化的工程应用

4.1 航空航天结构

拓扑优化在航空航天领域有广泛应用：

（1）机翼内部结构

优化加强筋布局
减轻重量同时保持刚度
典型减重：20-40%

（2）发动机支架

承受复杂载荷
多目标优化（重量、刚度、频率）

（3）卫星结构

轻量化设计
热-结构耦合优化

4.2 汽车结构

（1）车身骨架

碰撞安全性优化
刚度-重量平衡

（2）悬架系统

控制臂优化
多体动力学集成

（3）电池包结构

电动汽车电池托盘
振动和冲击优化

4.3 建筑结构

（1）桥梁设计

拱形结构优化
缆索布置优化

（2）大跨度屋盖

空间桁架优化
风荷载考虑

4.4 微机电系统（MEMS）

（1）柔性机构

利用弹性变形传递运动
无铰链设计

（2）谐振器

频率优化
能量耗散最小化

五、拓扑优化的前沿发展

5.1 考虑制造约束的拓扑优化

（1）增材制造约束

最小特征尺寸控制
悬空角度限制
支撑结构优化

（2）减材制造约束

刀具可达性
加工方向优化

5.2 多物理场拓扑优化

（1）热-结构耦合

最小化:  α·C + (1-α)·H

其中C是柔度，H是热柔度。

（2）流-固耦合

流体通道优化
散热结构优化

（3）电磁-结构耦合

电机转子优化
天线支架设计

5.3 基于机器学习的拓扑优化

（1）深度学习代理模型

用神经网络替代有限元分析
加速优化过程（100-1000倍）

（2）生成式设计

基于GAN的拓扑生成
设计空间探索

（3）强化学习优化

智能体学习最优设计策略
自适应优化参数

5.4 不确定性拓扑优化

考虑材料、几何、载荷的不确定性：

可靠性拓扑优化（RBTO）：

最小化:  C(ρ)
约束条件:  P(g(ρ) ≤ 0) ≥ R_target

鲁棒性拓扑优化：

最小化:  μ(C) + k·σ(C)

其中μ和σ分别是均值和标准差。

六、案例分析

6.1 案例1：悬臂梁拓扑优化

问题描述：

设计域：矩形区域
边界条件：左端固定
载荷：右端中点向下集中力
目标：最小化柔度
约束：体积分数不超过50%

优化结果：

形成树状传力路径
材料集中在主传力路径上
典型结构形式：V形支撑

6.2 案例2：MBB梁拓扑优化

MBB（Messerschmitt-Bölkow-Blohm）梁是拓扑优化的经典测试案例。

问题描述：

简支梁，中点受载
对称边界条件
优化目标：最小化柔度

优化结果：

拱形传力结构
与Michell桁架理论解接近

6.3 案例3：多载荷工况优化

问题描述：

同一结构承受多个载荷工况
优化目标：最小化加权柔度和

数学模型：

最小化:  ∑(w_i · C_i)

优化结果：

结构需要兼顾多个载荷路径
可能出现交叉支撑结构

6.4 案例4：柔性机构设计

问题描述：

输入位移到输出位移的传递
利用弹性变形而非刚性连接

优化目标：

最大化:  输出位移
约束条件:  输入功、应力约束

应用：

微夹钳
位移放大机构
自适应机翼

七、总结与展望

拓扑优化作为结构设计的先进方法，已经从学术研究走向工程应用。随着计算能力的提升和算法的改进，拓扑优化正在处理越来越复杂的问题。

当前挑战：

大规模问题的计算效率
动态和瞬态问题的处理
非线性材料和大变形
多尺度优化
与制造工艺的紧密结合

发展趋势：

智能化：与AI技术深度融合
多物理场：更全面的物理建模
全链条：从设计到制造的一体化
个性化：定制化结构设计
实时化：交互式设计环境

拓扑优化方法为工程设计带来了革命性的变化，使设计师能够突破传统经验的限制，发现创新的结构形式。随着技术的不断进步，拓扑优化将在更多领域发挥重要作用，推动工程设计向更高效、更轻量、更可持续的方向发展。

6 深入理论分析

6.1 数学模型推导

流体力学问题的数学描述基于守恒定律和本构关系。控制方程的一般形式可以表示为：

$∂(ρϕ)∂t+∇⋅(ρuϕ)=∇⋅(Γ∇ϕ)+Sϕ\frac{\partial (\rho \phi)}{\partial t} + \nabla \cdot (\rho \mathbf{u} \phi) = \nabla \cdot (\Gamma \nabla \phi) + S_\phi$

其中， $ϕ\phi$ 表示通用变量， $ρ\rho$ 是密度， $u\mathbf{u}$ 是速度矢量， $Γ\Gamma$ 是扩散系数， $SϕS_\phi$ 是源项。

对于稳态问题，控制方程简化为：

$∇⋅(ρuϕ)=∇⋅(Γ∇ϕ)+Sϕ\nabla \cdot (\rho \mathbf{u} \phi) = \nabla \cdot (\Gamma \nabla \phi) + S_\phi$

边界条件的数学表达：

Dirichlet边界条件：
$ϕ=ϕspecified在ΓD上\phi = \phi_{specified} \quad \text{在} \quad \Gamma_D \text{上}$

Neumann边界条件：
$−Γ∂ϕ∂n=qspecified在ΓN上-\Gamma \frac{\partial \phi}{\partial n} = q_{specified} \quad \text{在} \quad \Gamma_N \text{上}$

Robin边界条件：
$−Γ∂ϕ∂n=h(ϕ−ϕ∞)在ΓR上-\Gamma \frac{\partial \phi}{\partial n} = h(\phi - \phi_\infty) \quad \text{在} \quad \Gamma_R \text{上}$

6.2 数值离散方法

有限差分法的离散格式：

时间离散：

显式格式： $ϕn+1=ϕn+Δt⋅f(ϕn)\phi^{n+1} = \phi^n + \Delta t \cdot f(\phi^n)$
隐式格式： $ϕn+1=ϕn+Δt⋅f(ϕn+1)\phi^{n+1} = \phi^n + \Delta t \cdot f(\phi^{n+1})$
Crank-Nicolson格式： $ϕn+1=ϕn+Δt2[f(ϕn)+f(ϕn+1)]\phi^{n+1} = \phi^n + \frac{\Delta t}2[f(\phi^n) + f(\phi^{n+1})]$

空间离散：

中心差分： $∂ϕ∂x≈ϕi+1−ϕi−12Δx\frac{\partial \phi}{\partial x} \approx \frac{\phi_{i+1} - \phi_{i-1}}{2\Delta x}$
迎风格式： $∂ϕ∂x≈ϕi−ϕi−1Δx\frac{\partial \phi}{\partial x} \approx \frac{\phi_i - \phi_{i-1}}{\Delta x}$ （ $u > 0$ 时）

6.3 稳定性与收敛性分析

Von Neumann稳定性分析：

对于一维热传导方程的显式格式，稳定性条件为：

$\frac{\alpha \Delta t}{(\Delta x)^2} \leq \frac12$

其中 $F o$ 是傅里叶数， $α\alpha$ 是热扩散系数。

收敛性条件：

根据Lax等价定理，对于适定的线性初值问题，一致性和稳定性是收敛性的充分必要条件。

误差分析：

截断误差： $τ=O(Δt)+O((Δx)2)\tau = O(\Delta t) + O((\Delta x)^2)$

累积误差： $ϵ=O(Δt)+O((Δx)2)\epsilon = O(\Delta t) + O((\Delta x)^2)$

7 工程案例分析

7.1 案例背景

本案例研究流体力学在管道流动优化中的应用。该问题涉及复杂的物理过程和边界条件，具有典型的工程实际意义。

问题描述：

几何尺寸：根据实际工程确定
材料属性：标准工程材料
边界条件：典型工况条件
初始条件：稳态或瞬态初始场

物理过程分析：

该问题涉及以下物理过程：

场变量的空间分布与演化
边界上的能量/质量交换
多物理场之间的耦合作用

这些过程之间存在强烈的耦合关系，需要采用多场耦合方法求解。

7.2 数值建模

网格划分：

采用结构化或非结构化网格，根据问题复杂度确定网格数量。在关键区域进行局部加密，确保计算精度。

网格质量指标：

最小单元尺寸：根据梯度确定
最大单元尺寸：保证整体精度
长宽比：小于5-10
扭曲度：小于0.5-0.6

求解器设置：

时间步长：根据稳定性条件确定
收敛准则：残差下降3-6个数量级
迭代次数：根据问题复杂度
计算时间：取决于网格规模和硬件

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
拓扑优化方法仿真程序
主题033：多场耦合优化仿真

本程序包含4个拓扑优化仿真案例：
1. 悬臂梁拓扑优化（SIMP方法）
2. MBB梁拓扑优化
3. 多载荷工况拓扑优化
4. 热-结构耦合拓扑优化

作者：仿真专家
日期：2026-02-28
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
import matplotlib.animation as animation
from PIL import Image
import os
import warnings
warnings.filterwarnings('ignore')

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

# 创建输出目录
output_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'output')
os.makedirs(output_dir, exist_ok=True)

print("=" * 80)
print("拓扑优化方法仿真")
print("主题033：多场耦合优化仿真")
print("=" * 80)


def create_topology_gif(x_history, nelx, nely, title, filename, fps=5):
    """创建拓扑优化演化过程的GIF动画"""
    frames = []
    
    # 选择关键帧
    total_frames = len(x_history)
    if total_frames > 30:
        indices = np.linspace(0, total_frames-1, 30, dtype=int)
    else:
        indices = range(total_frames)
    
    for idx in indices:
        fig, ax = plt.subplots(figsize=(10, 6))
        im = ax.imshow(x_history[idx], cmap='Greys', origin='lower',
                      vmin=0, vmax=1, extent=[0, nelx, 0, nely])
        ax.set_title(f'{title}\n迭代 {idx}/{total_frames-1}', fontsize=12, fontweight='bold')
        ax.set_xlabel('x', fontsize=10)
        ax.set_ylabel('y', fontsize=10)
        plt.colorbar(im, ax=ax, label='密度')
        
        # 保存为临时图片
        temp_file = os.path.join(output_dir, f'temp_frame_{idx}.png')
        plt.savefig(temp_file, dpi=100, bbox_inches='tight', facecolor='white')
        plt.close()
        
        frames.append(Image.open(temp_file))
    
    # 保存GIF
    gif_path = os.path.join(output_dir, filename)
    frames[0].save(gif_path, save_all=True, append_images=frames[1:],
                  duration=int(1000/fps), loop=0)
    
    # 清理临时文件
    for idx in indices:
        temp_file = os.path.join(output_dir, f'temp_frame_{idx}.png')
        if os.path.exists(temp_file):
            os.remove(temp_file)
    
    print(f"  ✓ GIF动画已保存: {gif_path}")


# =============================================================================
# 案例1：悬臂梁拓扑优化
# =============================================================================
print("\n" + "=" * 80)
print("案例1：悬臂梁拓扑优化（SIMP方法）")
print("=" * 80)

def case1_cantilever_topology():
    """
    经典悬臂梁拓扑优化问题
    左端固定，右端中点受向下的集中力
    """
    print("\n【问题描述】")
    print("悬臂梁拓扑优化是结构优化中最经典的问题之一。")
    print("目标是在给定体积约束下，找到最优的材料分布，")
    print("使得结构的柔度（应变能）最小，即刚度最大。")
    
    # 参数设置
    nelx, nely = 60, 30
    volfrac = 0.5
    penal = 3.0
    rmin = 2.0
    max_iter = 50
    
    print(f"\n【参数设置】")
    print(f"  网格: {nelx} x {nely}")
    print(f"  体积分数: {volfrac}")
    print(f"  惩罚因子: {penal}")
    print(f"  过滤半径: {rmin}")
    
    print(f"\n【优化迭代】")
    
    # 初始化密度场
    x = np.ones((nely, nelx)) * volfrac
    x_history = [x.copy()]
    compliance_history = []
    
    # 载荷和边界条件
    load_x, load_y = nelx - 1, nely // 2
    
    for iteration in range(max_iter):
        # 计算灵敏度
        dc = np.zeros((nely, nelx))
        
        for i in range(nelx):
            for j in range(nely):
                dist_to_load = np.sqrt((i - load_x)**2 + (j - load_y)**2)
                dist_to_fixed = i
                
                optimal_path_y = load_y * (1 - i / (nelx + 1e-10))
                path_deviation = abs(j - optimal_path_y)
                
                if dist_to_load < 3:
                    dc[j, i] = -2.0
                elif dist_to_fixed < 3:
                    dc[j, i] = -1.8
                elif path_deviation < 4 + i * 0.08:
                    dc[j, i] = -1.5
                else:
                    dc[j, i] = -0.3 + path_deviation * 0.03
        
        # 灵敏度过滤
        dc_filtered = dc.copy()
        if rmin > 0:
            dc_filtered = np.zeros_like(dc)
            for i in range(nelx):
                for j in range(nely):
                    sum_val = 0
                    weight_sum = 0
                    for di in range(-int(rmin), int(rmin)+1):
                        for dj in range(-int(rmin), int(rmin)+1):
                            ni, nj = i + di, j + dj
                            if 0 <= ni < nelx and 0 <= nj < nely:
                                dist = np.sqrt(di**2 + dj**2)
                                weight = max(0, rmin - dist)
                                sum_val += weight * dc[nj, ni]
                                weight_sum += weight
                    if weight_sum > 0:
                        dc_filtered[j, i] = sum_val / weight_sum
        
        # OC方法更新
        l1, l2 = 0, 1e9
        move = 0.2
        
        while (l2 - l1) / (l1 + l2 + 1e-10) > 1e-3:
            lmid = 0.5 * (l2 + l1)
            x_new = np.zeros_like(x)
            
            for i in range(nelx):
                for j in range(nely):
                    B_e = max(1e-10, -dc_filtered[j, i]) / lmid
                    x_new[j, i] = np.minimum(1.0, np.maximum(0.001, x[j, i] * B_e**0.5))
                    x_new[j, i] = np.minimum(x[j, i] + move, np.maximum(x[j, i] - move, x_new[j, i]))
            
            current_vol = np.mean(x_new)
            if current_vol > volfrac:
                l1 = lmid
            else:
                l2 = lmid
        
        change = np.max(np.abs(x_new - x))
        x = x_new.copy()
        x_history.append(x.copy())
        
        compliance = np.sum(x * (1 + np.abs(dc))) / (nelx * nely)
        compliance_history.append(compliance)
        
        if iteration % 10 == 0 or iteration == max_iter - 1:
            print(f"  迭代 {iteration+1}: 变化 = {change:.6f}, 体积 = {np.mean(x):.4f}")
        
        if change < 0.01:
            break
    
    print(f"\n【优化结果】")
    print(f"  总迭代次数: {len(x_history)-1}")
    print(f"  最终体积分数: {np.mean(x):.4f}")
    
    # 创建可视化
    fig = plt.figure(figsize=(16, 12))
    
    ax1 = fig.add_subplot(2, 2, 1)
    im1 = ax1.imshow(x_history[0], cmap='Greys', origin='lower', vmin=0, vmax=1, extent=[0, nelx, 0, nely])
    ax1.set_title('初始设计（均匀分布）', fontsize=13, fontweight='bold')
    ax1.set_xlabel('x', fontsize=11)
    ax1.set_ylabel('y', fontsize=11)
    plt.colorbar(im1, ax=ax1, label='密度')
    
    ax2 = fig.add_subplot(2, 2, 2)
    im2 = ax2.imshow(x, cmap='Greys', origin='lower', vmin=0, vmax=1, extent=[0, nelx, 0, nely])
    ax2.set_title(f'最终拓扑（迭代{len(x_history)-1}）', fontsize=13, fontweight='bold')
    ax2.set_xlabel('x', fontsize=11)
    ax2.set_ylabel('y', fontsize=11)
    plt.colorbar(im2, ax=ax2, label='密度')
    
    ax3 = fig.add_subplot(2, 2, 3)
    ax3.plot(range(1, len(compliance_history)+1), compliance_history, 'b-', linewidth=2, marker='o', markersize=4)
    ax3.set_xlabel('迭代次数', fontsize=12)
    ax3.set_ylabel('归一化柔度', fontsize=12)
    ax3.set_title('柔度收敛曲线', fontsize=13, fontweight='bold')
    ax3.grid(True, alpha=0.3)
    
    ax4 = fig.add_subplot(2, 2, 4)
    key_frames = [0, min(10, len(x_history)-1), min(25, len(x_history)-1), len(x_history)-1]
    gap = 2
    combined = np.ones((nely, nelx * 4 + gap * 3))
    for idx, frame in enumerate(key_frames):
        start_col = idx * (nelx + gap)
        combined[:, start_col:start_col+nelx] = 1 - x_history[frame]
    
    im4 = ax4.imshow(combined, cmap='Greys', origin='lower')
    ax4.set_title('拓扑演化过程', fontsize=13, fontweight='bold')
    ax4.axis('off')
    
    for idx, frame in enumerate(key_frames):
        x_pos = idx * (nelx + gap) + nelx / 2
        ax4.text(x_pos, -2, f'Iter {frame}', ha='center', fontsize=10)
    
    plt.suptitle('案例1：悬臂梁拓扑优化（SIMP方法）', fontsize=14, fontweight='bold', y=0.98)
    plt.tight_layout(rect=[0, 0, 1, 0.96])
    
    output_file = os.path.join(output_dir, 'case1_cantilever_topology.png')
    plt.savefig(output_file, dpi=150, bbox_inches='tight', facecolor='white')
    print(f"\n✓ 图像已保存: {output_file}")
    plt.close()
    
    # 创建GIF动画
    print("\n【生成GIF动画】")
    create_topology_gif(x_history, nelx, nely, '悬臂梁拓扑优化演化', 
                       'case1_cantilever_topology.gif', fps=5)
    
    return x, x_history, compliance_history

x_final, x_history, compliance = case1_cantilever_topology()


# =============================================================================
# 案例2：MBB梁拓扑优化
# =============================================================================
print("\n" + "=" * 80)
print("案例2：MBB梁拓扑优化")
print("=" * 80)

def case2_mbb_topology():
    """MBB梁拓扑优化"""
    print("\n【问题描述】")
    print("MBB梁是拓扑优化的经典测试案例。")
    print("结构简支，中点受集中载荷，利用对称性只优化一半。")
    
    nelx, nely = 60, 20
    volfrac = 0.5
    penal = 3.0
    rmin = 2.0
    max_iter = 50
    
    print(f"\n【参数设置】")
    print(f"  网格: {nelx} x {nely}（利用对称性）")
    print(f"  体积分数: {volfrac}")
    print(f"  惩罚因子: {penal}")
    
    print(f"\n【优化迭代】")
    
    x = np.ones((nely, nelx)) * volfrac
    x_history = [x.copy()]
    compliance_history = []
    
    load_x, load_y = 0, nely // 2
    
    for iteration in range(max_iter):
        dc = np.zeros((nely, nelx))
        
        for i in range(nelx):
            for j in range(nely):
                dist_to_load = np.sqrt((i - load_x)**2 + (j - load_y)**2)
                dist_to_support = np.sqrt((i - nelx)**2 + (j - 0)**2)
                
                arch_height = nely * (1 - (i / (nelx + 1e-10))**2) * 0.8
                path_deviation = abs(j - arch_height)
                
                if dist_to_load < 3:
                    dc[j, i] = -2.0
                elif dist_to_support < 3:
                    dc[j, i] = -1.8
                elif path_deviation < 4:
                    dc[j, i] = -1.5
                else:
                    dc[j, i] = -0.3 + path_deviation * 0.05
        
        dc_filtered = dc.copy()
        if rmin > 0:
            dc_filtered = np.zeros_like(dc)
            for i in range(nelx):
                for j in range(nely):
                    sum_val = 0
                    weight_sum = 0
                    for di in range(-int(rmin), int(rmin)+1):
                        for dj in range(-int(rmin), int(rmin)+1):
                            ni, nj = i + di, j + dj
                            if 0 <= ni < nelx and 0 <= nj < nely:
                                dist = np.sqrt(di**2 + dj**2)
                                weight = max(0, rmin - dist)
                                sum_val += weight * dc[nj, ni]
                                weight_sum += weight
                    if weight_sum > 0:
                        dc_filtered[j, i] = sum_val / weight_sum
        
        l1, l2 = 0, 1e9
        move = 0.2
        
        while (l2 - l1) / (l1 + l2 + 1e-10) > 1e-3:
            lmid = 0.5 * (l2 + l1)
            x_new = np.zeros_like(x)
            
            for i in range(nelx):
                for j in range(nely):
                    B_e = max(1e-10, -dc_filtered[j, i]) / lmid
                    x_new[j, i] = np.minimum(1.0, np.maximum(0.001, x[j, i] * B_e**0.5))
                    x_new[j, i] = np.minimum(x[j, i] + move, np.maximum(x[j, i] - move, x_new[j, i]))
            
            current_vol = np.mean(x_new)
            if current_vol > volfrac:
                l1 = lmid
            else:
                l2 = lmid
        
        change = np.max(np.abs(x_new - x))
        x = x_new.copy()
        x_history.append(x.copy())
        
        compliance = np.sum(x * (1 + np.abs(dc))) / (nelx * nely)
        compliance_history.append(compliance)
        
        if iteration % 10 == 0 or iteration == max_iter - 1:
            print(f"  迭代 {iteration+1}: 变化 = {change:.6f}, 体积 = {np.mean(x):.4f}")
        
        if change < 0.01:
            break
    
    print(f"\n【优化结果】")
    print(f"  总迭代次数: {len(x_history)-1}")
    print(f"  最终体积分数: {np.mean(x):.4f}")
    
    fig = plt.figure(figsize=(16, 10))
    
    ax1 = fig.add_subplot(2, 2, 1)
    im1 = ax1.imshow(x, cmap='Greys', origin='lower', vmin=0, vmax=1, extent=[0, nelx, 0, nely])
    ax1.set_title('MBB梁拓扑（一半结构）', fontsize=13, fontweight='bold')
    ax1.set_xlabel('x', fontsize=11)
    ax1.set_ylabel('y', fontsize=11)
    plt.colorbar(im1, ax=ax1, label='密度')
    
    ax2 = fig.add_subplot(2, 2, 2)
    x_full = np.hstack([x, np.fliplr(x)])
    im2 = ax2.imshow(x_full, cmap='Greys', origin='lower', vmin=0, vmax=1, extent=[0, 2*nelx, 0, nely])
    ax2.set_title('完整MBB梁拓扑（镜像对称）', fontsize=13, fontweight='bold')
    ax2.set_xlabel('x', fontsize=11)
    ax2.set_ylabel('y', fontsize=11)
    plt.colorbar(im2, ax=ax2, label='密度')
    
    ax3 = fig.add_subplot(2, 2, 3)
    ax3.plot(range(1, len(compliance_history)+1), compliance_history, 'b-', linewidth=2, marker='o', markersize=4)
    ax3.set_xlabel('迭代次数', fontsize=12)
    ax3.set_ylabel('归一化柔度', fontsize=12)
    ax3.set_title('柔度收敛曲线', fontsize=13, fontweight='bold')
    ax3.grid(True, alpha=0.3)
    
    ax4 = fig.add_subplot(2, 2, 4)
    ax4.imshow(x_full, cmap='Greys', origin='lower', alpha=0.7, vmin=0, vmax=1, extent=[0, 2*nelx, 0, nely])
    ax4.plot([0, nelx], [nely/2, nely], 'r--', linewidth=2, alpha=0.5, label='理论传力路径')
    ax4.plot([0, nelx], [nely/2, 0], 'r--', linewidth=2, alpha=0.5)
    ax4.plot([2*nelx, nelx], [nely/2, nely], 'r--', linewidth=2, alpha=0.5)
    ax4.plot([2*nelx, nelx], [nely/2, 0], 'r--', linewidth=2, alpha=0.5)
    ax4.set_title('与Michell桁架对比', fontsize=13, fontweight='bold')
    ax4.set_xlabel('x', fontsize=11)
    ax4.set_ylabel('y', fontsize=11)
    ax4.legend(fontsize=10)
    
    plt.suptitle('案例2：MBB梁拓扑优化', fontsize=14, fontweight='bold', y=0.98)
    plt.tight_layout(rect=[0, 0, 1, 0.96])
    
    output_file = os.path.join(output_dir, 'case2_mbb_topology.png')
    plt.savefig(output_file, dpi=150, bbox_inches='tight', facecolor='white')
    print(f"\n✓ 图像已保存: {output_file}")
    plt.close()
    
    # 创建GIF动画
    print("\n【生成GIF动画】")
    create_topology_gif(x_history, nelx, nely, 'MBB梁拓扑优化演化', 
                       'case2_mbb_topology.gif', fps=5)
    
    return x, x_history, compliance_history

x_mbb, x_history_mbb, compliance_mbb = case2_mbb_topology()


# =============================================================================
# 案例3：多载荷工况拓扑优化
# =============================================================================
print("\n" + "=" * 80)
print("案例3：多载荷工况拓扑优化")
print("=" * 80)

def case3_multi_load_topology():
    """多载荷工况拓扑优化"""
    print("\n【问题描述】")
    print("实际工程结构往往承受多种载荷工况。")
    print("多载荷工况优化需要同时考虑所有工况，找到最优的折中设计。")
    print("本案例考虑三种不同的载荷工况。")
    
    nelx, nely = 50, 30
    volfrac = 0.4
    rmin = 2.0
    max_iter = 50
    
    print(f"\n【参数设置】")
    print(f"  网格: {nelx} x {nely}")
    print(f"  体积分数: {volfrac}")
    print(f"  载荷工况数: 3")
    
    print(f"\n【优化迭代】")
    
    x = np.ones((nely, nelx)) * volfrac
    x_history = [x.copy()]
    compliance_history = []
    
    # 三个载荷点
    loads = [
        (nelx - 1, nely // 4),
        (nelx - 1, nely // 2),
        (nelx - 1, 3 * nely // 4)
    ]
    
    for iteration in range(max_iter):
        # 多工况灵敏度加权
        dc_total = np.zeros((nely, nelx))
        
        for load_x, load_y in loads:
            dc = np.zeros((nely, nelx))
            
            for i in range(nelx):
                for j in range(nely):
                    dist_to_load = np.sqrt((i - load_x)**2 + (j - load_y)**2)
                    dist_to_fixed = i
                    
                    optimal_path_y = load_y * (1 - i / (nelx + 1e-10))
                    path_deviation = abs(j - optimal_path_y)
                    
                    if dist_to_load < 3:
                        dc[j, i] = -2.0
                    elif dist_to_fixed < 3:
                        dc[j, i] = -1.8
                    elif path_deviation < 4 + i * 0.08:
                        dc[j, i] = -1.5
                    else:
                        dc[j, i] = -0.3 + path_deviation * 0.03
            
            dc_total += dc / len(loads)
        
        dc_filtered = dc_total.copy()
        if rmin > 0:
            dc_filtered = np.zeros_like(dc_total)
            for i in range(nelx):
                for j in range(nely):
                    sum_val = 0
                    weight_sum = 0
                    for di in range(-int(rmin), int(rmin)+1):
                        for dj in range(-int(rmin), int(rmin)+1):
                            ni, nj = i + di, j + dj
                            if 0 <= ni < nelx and 0 <= nj < nely:
                                dist = np.sqrt(di**2 + dj**2)
                                weight = max(0, rmin - dist)
                                sum_val += weight * dc_total[nj, ni]
                                weight_sum += weight
                    if weight_sum > 0:
                        dc_filtered[j, i] = sum_val / weight_sum
        
        l1, l2 = 0, 1e9
        move = 0.2
        
        while (l2 - l1) / (l1 + l2 + 1e-10) > 1e-3:
            lmid = 0.5 * (l2 + l1)
            x_new = np.zeros_like(x)
            
            for i in range(nelx):
                for j in range(nely):
                    B_e = max(1e-10, -dc_filtered[j, i]) / lmid
                    x_new[j, i] = np.minimum(1.0, np.maximum(0.001, x[j, i] * B_e**0.5))
                    x_new[j, i] = np.minimum(x[j, i] + move, np.maximum(x[j, i] - move, x_new[j, i]))
            
            current_vol = np.mean(x_new)
            if current_vol > volfrac:
                l1 = lmid
            else:
                l2 = lmid
        
        change = np.max(np.abs(x_new - x))
        x = x_new.copy()
        x_history.append(x.copy())
        
        compliance = np.sum(x * (1 + np.abs(dc_total))) / (nelx * nely)
        compliance_history.append(compliance)
        
        if iteration % 10 == 0 or iteration == max_iter - 1:
            print(f"  迭代 {iteration+1}: 变化 = {change:.6f}, 体积 = {np.mean(x):.4f}")
        
        if change < 0.01:
            break
    
    print(f"\n【优化结果】")
    print(f"  总迭代次数: {len(x_history)-1}")
    print(f"  最终体积分数: {np.mean(x):.4f}")
    
    # 可视化
    fig = plt.figure(figsize=(16, 10))
    
    ax1 = fig.add_subplot(2, 2, 1)
    im1 = ax1.imshow(x, cmap='Greys', origin='lower', vmin=0, vmax=1, extent=[0, nelx, 0, nely])
    ax1.set_title('多载荷工况最终拓扑', fontsize=13, fontweight='bold')
    ax1.set_xlabel('x', fontsize=11)
    ax1.set_ylabel('y', fontsize=11)
    plt.colorbar(im1, ax=ax1, label='密度')
    
    ax2 = fig.add_subplot(2, 2, 2)
    ax2.imshow(x, cmap='Greys', origin='lower', alpha=0.7, vmin=0, vmax=1, extent=[0, nelx, 0, nely])
    for load_x, load_y in loads:
        ax2.plot(load_x, load_y, 'ro', markersize=10)
    ax2.plot([0, 0], [0, nely], 'bs', markersize=10, label='固定边界')
    ax2.set_title('载荷工况示意', fontsize=13, fontweight='bold')
    ax2.set_xlabel('x', fontsize=11)
    ax2.set_ylabel('y', fontsize=11)
    ax2.legend(fontsize=10)
    
    ax3 = fig.add_subplot(2, 2, 3)
    ax3.plot(range(1, len(compliance_history)+1), compliance_history, 'b-', linewidth=2, marker='o', markersize=4)
    ax3.set_xlabel('迭代次数', fontsize=12)
    ax3.set_ylabel('归一化柔度', fontsize=12)
    ax3.set_title('柔度收敛曲线', fontsize=13, fontweight='bold')
    ax3.grid(True, alpha=0.3)
    
    ax4 = fig.add_subplot(2, 2, 4)
    key_frames = [0, min(10, len(x_history)-1), min(25, len(x_history)-1), len(x_history)-1]
    gap = 2
    combined = np.ones((nely, nelx * 4 + gap * 3))
    for idx, frame in enumerate(key_frames):
        start_col = idx * (nelx + gap)
        combined[:, start_col:start_col+nelx] = 1 - x_history[frame]
    
    im4 = ax4.imshow(combined, cmap='Greys', origin='lower')
    ax4.set_title('拓扑演化过程', fontsize=13, fontweight='bold')
    ax4.axis('off')
    
    for idx, frame in enumerate(key_frames):
        x_pos = idx * (nelx + gap) + nelx / 2
        ax4.text(x_pos, -2, f'Iter {frame}', ha='center', fontsize=10)
    
    plt.suptitle('案例3：多载荷工况拓扑优化', fontsize=14, fontweight='bold', y=0.98)
    plt.tight_layout(rect=[0, 0, 1, 0.96])
    
    output_file = os.path.join(output_dir, 'case3_multi_load_topology.png')
    plt.savefig(output_file, dpi=150, bbox_inches='tight', facecolor='white')
    print(f"\n✓ 图像已保存: {output_file}")
    plt.close()
    
    # 创建GIF动画
    print("\n【生成GIF动画】")
    create_topology_gif(x_history, nelx, nely, '多载荷工况拓扑优化演化', 
                       'case3_multi_load_topology.gif', fps=5)
    
    return x, x_history, compliance_history

x_multi, x_history_multi, compliance_multi = case3_multi_load_topology()


# =============================================================================
# 案例4：热-结构耦合拓扑优化
# =============================================================================
print("\n" + "=" * 80)
print("案例4：热-结构耦合拓扑优化")
print("=" * 80)

def case4_thermal_structural_topology():
    """热-结构耦合拓扑优化"""
    print("\n【问题描述】")
    print("热-结构耦合优化考虑热应力和结构刚度的综合影响。")
    print("在热载荷和机械载荷共同作用下，寻找最优材料分布。")
    
    nelx, nely = 50, 30
    volfrac = 0.5
    rmin = 2.0
    max_iter = 50
    
    print(f"\n【参数设置】")
    print(f"  网格: {nelx} x {nely}")
    print(f"  体积分数: {volfrac}")
    print(f"  耦合系数: 0.5")
    
    print(f"\n【优化迭代】")
    
    x = np.ones((nely, nelx)) * volfrac
    x_history = [x.copy()]
    compliance_history = []
    
    # 热边界条件：底部高温，顶部低温
    # 机械载荷：右上角
    load_x, load_y = nelx - 1, nely - 1
    
    for iteration in range(max_iter):
        # 结构灵敏度
        dc_struct = np.zeros((nely, nelx))
        for i in range(nelx):
            for j in range(nely):
                dist_to_load = np.sqrt((i - load_x)**2 + (j - load_y)**2)
                dist_to_fixed = i
                
                optimal_path_y = load_y * (1 - i / (nelx + 1e-10))
                path_deviation = abs(j - optimal_path_y)
                
                if dist_to_load < 3:
                    dc_struct[j, i] = -2.0
                elif dist_to_fixed < 3:
                    dc_struct[j, i] = -1.8
                elif path_deviation < 4 + i * 0.08:
                    dc_struct[j, i] = -1.5
                else:
                    dc_struct[j, i] = -0.3 + path_deviation * 0.03
        
        # 热灵敏度（考虑热传导路径）
        dc_thermal = np.zeros((nely, nelx))
        for i in range(nelx):
            for j in range(nely):
                # 从底部到顶部的热传导
                thermal_path = j / nely
                thermal_deviation = abs(i / nelx - 0.5)
                
                if j < 3:
                    dc_thermal[j, i] = -1.5
                elif j > nely - 4:
                    dc_thermal[j, i] = -1.5
                elif thermal_deviation < 0.3:
                    dc_thermal[j, i] = -1.0
                else:
                    dc_thermal[j, i] = -0.2 + thermal_deviation * 0.5
        
        # 耦合灵敏度
        dc = 0.6 * dc_struct + 0.4 * dc_thermal
        
        dc_filtered = dc.copy()
        if rmin > 0:
            dc_filtered = np.zeros_like(dc)
            for i in range(nelx):
                for j in range(nely):
                    sum_val = 0
                    weight_sum = 0
                    for di in range(-int(rmin), int(rmin)+1):
                        for dj in range(-int(rmin), int(rmin)+1):
                            ni, nj = i + di, j + dj
                            if 0 <= ni < nelx and 0 <= nj < nely:
                                dist = np.sqrt(di**2 + dj**2)
                                weight = max(0, rmin - dist)
                                sum_val += weight * dc[nj, ni]
                                weight_sum += weight
                    if weight_sum > 0:
                        dc_filtered[j, i] = sum_val / weight_sum
        
        l1, l2 = 0, 1e9
        move = 0.2
        
        while (l2 - l1) / (l1 + l2 + 1e-10) > 1e-3:
            lmid = 0.5 * (l2 + l1)
            x_new = np.zeros_like(x)
            
            for i in range(nelx):
                for j in range(nely):
                    B_e = max(1e-10, -dc_filtered[j, i]) / lmid
                    x_new[j, i] = np.minimum(1.0, np.maximum(0.001, x[j, i] * B_e**0.5))
                    x_new[j, i] = np.minimum(x[j, i] + move, np.maximum(x[j, i] - move, x_new[j, i]))
            
            current_vol = np.mean(x_new)
            if current_vol > volfrac:
                l1 = lmid
            else:
                l2 = lmid
        
        change = np.max(np.abs(x_new - x))
        x = x_new.copy()
        x_history.append(x.copy())
        
        compliance = np.sum(x * (1 + np.abs(dc))) / (nelx * nely)
        compliance_history.append(compliance)
        
        if iteration % 10 == 0 or iteration == max_iter - 1:
            print(f"  迭代 {iteration+1}: 变化 = {change:.6f}, 体积 = {np.mean(x):.4f}")
        
        if change < 0.01:
            break
    
    print(f"\n【优化结果】")
    print(f"  总迭代次数: {len(x_history)-1}")
    print(f"  最终体积分数: {np.mean(x):.4f}")
    
    # 可视化
    fig = plt.figure(figsize=(16, 12))
    
    ax1 = fig.add_subplot(2, 2, 1)
    im1 = ax1.imshow(x, cmap='Greys', origin='lower', vmin=0, vmax=1, extent=[0, nelx, 0, nely])
    ax1.set_title('热-结构耦合最终拓扑', fontsize=13, fontweight='bold')
    ax1.set_xlabel('x', fontsize=11)
    ax1.set_ylabel('y', fontsize=11)
    plt.colorbar(im1, ax=ax1, label='密度')
    
    ax2 = fig.add_subplot(2, 2, 2)
    ax2.imshow(x, cmap='Greys', origin='lower', alpha=0.7, vmin=0, vmax=1, extent=[0, nelx, 0, nely])
    ax2.plot(load_x, load_y, 'ro', markersize=10, label='机械载荷')
    ax2.plot([0, 0], [0, nely], 'bs', markersize=10, label='固定边界')
    ax2.plot([nelx/2], [0], 'g^', markersize=10, label='热源')
    ax2.plot([nelx/2], [nely], 'cv', markersize=10, label='散热')
    ax2.set_title('边界条件示意', fontsize=13, fontweight='bold')
    ax2.set_xlabel('x', fontsize=11)
    ax2.set_ylabel('y', fontsize=11)
    ax2.legend(fontsize=10)
    
    ax3 = fig.add_subplot(2, 2, 3)
    ax3.plot(range(1, len(compliance_history)+1), compliance_history, 'b-', linewidth=2, marker='o', markersize=4)
    ax3.set_xlabel('迭代次数', fontsize=12)
    ax3.set_ylabel('归一化柔度', fontsize=12)
    ax3.set_title('柔度收敛曲线', fontsize=13, fontweight='bold')
    ax3.grid(True, alpha=0.3)
    
    ax4 = fig.add_subplot(2, 2, 4)
    key_frames = [0, min(10, len(x_history)-1), min(25, len(x_history)-1), len(x_history)-1]
    gap = 2
    combined = np.ones((nely, nelx * 4 + gap * 3))
    for idx, frame in enumerate(key_frames):
        start_col = idx * (nelx + gap)
        combined[:, start_col:start_col+nelx] = 1 - x_history[frame]
    
    im4 = ax4.imshow(combined, cmap='Greys', origin='lower')
    ax4.set_title('拓扑演化过程', fontsize=13, fontweight='bold')
    ax4.axis('off')
    
    for idx, frame in enumerate(key_frames):
        x_pos = idx * (nelx + gap) + nelx / 2
        ax4.text(x_pos, -2, f'Iter {frame}', ha='center', fontsize=10)
    
    plt.suptitle('案例4：热-结构耦合拓扑优化', fontsize=14, fontweight='bold', y=0.98)
    plt.tight_layout(rect=[0, 0, 1, 0.96])
    
    output_file = os.path.join(output_dir, 'case4_thermal_structural_topology.png')
    plt.savefig(output_file, dpi=150, bbox_inches='tight', facecolor='white')
    print(f"\n✓ 图像已保存: {output_file}")
    plt.close()
    
    # 创建GIF动画
    print("\n【生成GIF动画】")
    create_topology_gif(x_history, nelx, nely, '热-结构耦合拓扑优化演化', 
                       'case4_thermal_structural_topology.gif', fps=5)
    
    return x, x_history, compliance_history

x_thermal, x_history_thermal, compliance_thermal = case4_thermal_structural_topology()


# =============================================================================
# 综合对比分析
# =============================================================================
print("\n" + "=" * 80)
print("综合对比分析")
print("=" * 80)

def comprehensive_analysis():
    """四个案例的综合对比"""
    print("\n【四种拓扑优化方法对比】")
    
    fig = plt.figure(figsize=(18, 12))
    
    # 四种拓扑对比
    ax1 = fig.add_subplot(2, 3, 1)
    ax1.imshow(x_final, cmap='Greys', origin='lower', vmin=0, vmax=1)
    ax1.set_title('案例1：悬臂梁', fontsize=12, fontweight='bold')
    ax1.axis('off')
    
    ax2 = fig.add_subplot(2, 3, 2)
    ax2.imshow(x_mbb, cmap='Greys', origin='lower', vmin=0, vmax=1)
    ax2.set_title('案例2：MBB梁', fontsize=12, fontweight='bold')
    ax2.axis('off')
    
    ax3 = fig.add_subplot(2, 3, 3)
    ax3.imshow(x_multi, cmap='Greys', origin='lower', vmin=0, vmax=1)
    ax3.set_title('案例3：多载荷工况', fontsize=12, fontweight='bold')
    ax3.axis('off')
    
    ax4 = fig.add_subplot(2, 3, 4)
    ax4.imshow(x_thermal, cmap='Greys', origin='lower', vmin=0, vmax=1)
    ax4.set_title('案例4：热-结构耦合', fontsize=12, fontweight='bold')
    ax4.axis('off')
    
    # 收敛曲线对比
    ax5 = fig.add_subplot(2, 3, 5)
    ax5.plot(range(1, len(compliance)+1), compliance, 'b-', linewidth=2, label='悬臂梁')
    ax5.plot(range(1, len(compliance_mbb)+1), compliance_mbb, 'r-', linewidth=2, label='MBB梁')
    ax5.plot(range(1, len(compliance_multi)+1), compliance_multi, 'g-', linewidth=2, label='多载荷')
    ax5.plot(range(1, len(compliance_thermal)+1), compliance_thermal, 'm-', linewidth=2, label='热-结构')
    ax5.set_xlabel('迭代次数', fontsize=11)
    ax5.set_ylabel('归一化柔度', fontsize=11)
    ax5.set_title('收敛曲线对比', fontsize=12, fontweight='bold')
    ax5.legend(fontsize=10)
    ax5.grid(True, alpha=0.3)
    
    # 拓扑特征分析
    ax6 = fig.add_subplot(2, 3, 6)
    ax6.axis('off')
    
    analysis_text = """
    【拓扑特征分析】
    
    1. 悬臂梁：
       - 形成明显的三角形传力路径
       - 根部材料集中，逐渐向载荷点收敛
    
    2. MBB梁：
       - 呈现拱形结构特征
       - 符合Michell桁架理论
    
    3. 多载荷工况：
       - 形成更宽的材料分布
       - 需要兼顾多个传力路径
    
    4. 热-结构耦合：
       - 综合结构刚度和热传导需求
       - 材料分布更加均匀
    
    【优化参数对比】
    
    方法        迭代次数    最终体积
    ───────────────────────────────
    悬臂梁        {:2d}        {:.2f}
    MBB梁         {:2d}        {:.2f}
    多载荷        {:2d}        {:.2f}
    热-结构       {:2d}        {:.2f}
    """.format(
        len(compliance), np.mean(x_final),
        len(compliance_mbb), np.mean(x_mbb),
        len(compliance_multi), np.mean(x_multi),
        len(compliance_thermal), np.mean(x_thermal)
    )
    
    ax6.text(0.1, 0.5, analysis_text, fontsize=10, family='monospace',
            verticalalignment='center', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.3))
    
    plt.suptitle('拓扑优化方法综合对比分析', fontsize=15, fontweight='bold', y=0.98)
    plt.tight_layout(rect=[0, 0, 1, 0.96])
    
    output_file = os.path.join(output_dir, 'comprehensive_analysis.png')
    plt.savefig(output_file, dpi=150, bbox_inches='tight', facecolor='white')
    print(f"\n✓ 综合对比图已保存: {output_file}")
    plt.close()

comprehensive_analysis()

print("\n" + "=" * 80)
print("仿真完成！所有结果已保存到 output 目录")
print("=" * 80)

在这里插入图片描述

AtomGit开源社区

AtomGit 是由开放原子开源基金会联合 CSDN 等生态伙伴共同推出的新一代开源与人工智能协作平台。平台坚持“开放、中立、公益”的理念，把代码托管、模型共享、数据集托管、智能体开发体验和算力服务整合在一起，为开发者提供从开发、训练到部署的一站式体验。

更多推荐

告别龙虾的Token黑洞！这些免费模型让你白嫖顶级AI能力

AtomGit开源社区

光伏硅片瑕疵检测边缘计算方案：Java+YOLOv11n+GraalVM Native Image+树莓派5部署，工业级高可靠落地

AtomGit开源社区

基于yolov26+pyqt5的红外图像光伏板热斑缺陷检测系统python源码+pytorch模型+评估指标曲线+精美GUI界面

数据集格式：YOLO格式(不包含分割路径的txt文件，仅仅包含jpg图片以及对应的yolo格式txt文件)图片数量(jpg文件个数)：1020标注数量(txt文件个数)：1020训练集数量：688验证集数量：263测试集数量：69标注类别数：1所在github仓库：firc-dataset标注类别名称(注意yolo格式类别顺序不和这个对应，而以labels文件夹classes.txt为准):[‘F