主题051：磁浮列车系统

场景描述

用一句话简述要解决的工程问题：

如何设计和分析磁浮列车的悬浮、导向、推进和制动系统，实现无接触、低摩擦、超高速的轨道交通？

---

一、引言

1.1 磁浮列车的发展历史

**磁浮列车（Maglev Train）**是一种利用磁力实现悬浮和推进的现代轨道交通工具。其发展历程可追溯至20世纪初：

（1）早期探索（1900-1960年代）

• 1909年：美国科学家Robert Goddard首次提出磁悬浮概念
• 1934年：德国工程师Hermann Kemper获得磁悬浮专利
• 1960年代：美国布鲁克海文国家实验室开始磁悬浮研究

（2）技术突破（1970-1990年代）

• 1979年：德国Transrapid技术验证车TR02达到80 km/h
• 1984年：英国伯明翰机场建成首条商业磁浮线
• 1990年代：日本MLX01超导磁浮列车达到500 km/h以上

（3）商业化运营（2000年代至今）

• 2004年：中国上海磁浮示范线开通，最高时速430 km/h
• 2015年：日本L0系超导磁浮列车达到603 km/h世界纪录
• 2020年代：中国时速600公里高速磁浮交通系统研发成功
  
  
  
  
  

1.2 磁浮列车的技术分类

磁浮列车根据悬浮原理可分为两大类：

（1）电磁悬浮（EMS - Electromagnetic Suspension）

• 原理：利用电磁铁吸力实现悬浮
• 特点：
  • 悬浮气隙小（8-12 mm）
  • 需要主动控制系统
  • 低速时即可悬浮
  • 代表：德国Transrapid、中国中低速磁浮

（2）电动悬浮（EDS - Electrodynamic Suspension）

• 原理：利用超导磁体与导电轨道的涡流排斥力实现悬浮
• 特点：
  • 悬浮气隙大（100-150 mm）
  • 自稳定，无需主动控制
  • 需要一定速度才能起飞（约100 km/h）
  • 代表：日本JR磁浮、美国Inductrack

1.3 磁浮列车的技术优势

（1）超高速运行

• 轮轨列车极限速度约350 km/h
• 磁浮列车可稳定运行500-600 km/h
• 上海磁浮示范线商业运营速度430 km/h

（2）低噪音低振动

• 无轮轨接触，消除轮轨噪声
• 空气动力学噪声成为主要噪声源
• 适合城市环境和敏感区域

（3）高爬坡能力

• 爬坡能力可达10%（轮轨约4%）
• 可适应复杂地形
• 减少隧道和桥梁建设

（4）低维护成本

• 无机械磨损部件
• 维护周期长
• 全寿命周期成本低

（5）高安全性

• 无脱轨风险
• 抱轨运行设计
• 自动控制系统

1.4 本教程的学习目标

通过本教程的学习，读者将能够：

1. 理解EMS和EDS两种悬浮系统的物理原理，掌握其数学模型和设计方法
2. 学会分析直线电机的推力特性，理解牵引和制动的工作原理
3. 掌握磁浮列车动力学建模方法，能够进行运行工况仿真
4. 了解涡流制动的原理和特性，能够设计制动系统
5. 具备使用Python进行磁浮系统仿真的能力，为实际工程设计提供支持

---

二、数学/物理模型

2.1 电磁悬浮（EMS）理论

2.1.1 电磁吸力公式

电磁铁对铁磁材料的吸力由麦克斯韦应力给出：

$$F=\frac{B^{2}A}{2{\mu }_0}=\frac{{\mu }_0{N}^2{I}^{2}A}{4g^2}$$

其中：

• $F$：电磁吸力 (N)
• $B$：气隙磁感应强度 (T)
• $A$：磁极面积 (m²)
• $N$：线圈匝数
• $I$：线圈电流 (A)
• $g$：气隙长度 (m)
• ${\mu }_0=4\pi \times 10^{-7}$ H/m：真空磁导率

关键特性：

• 吸力与电流平方成正比
• 吸力与气隙平方成反比
• 系统具有负刚度，需要主动控制

2.1.2 系统刚度分析

电磁悬浮系统的刚度定义为：

$$k=-\frac{\partial F}{\partial g}=\frac{{\mu }_0{N}^2{I}^{2}A}{2g^3}$$

负刚度问题：

• 当气隙增大时，吸力减小，导致气隙进一步增大
• 系统本身不稳定，需要主动控制
• 控制器通过调节电流维持稳定悬浮

2.1.3 主动控制系统

EMS系统的控制方程：

$$m\ddot{z}=mg-F(g,I)$$

线性化后的状态空间方程：

$$\dot{x}=Ax+Bu$$

其中状态变量 $x=[\Delta z,\Delta \dot{z}{]}^T$，控制输入 $u=\Delta I$。

PID控制器设计：

$$I(t)=I_0+K_{p}e(t)+K_i\int e(t)dt+K_d\frac{de(t)}{dt}$$

其中 $e(t)=z_{ref}-z(t)$ 是气隙误差。

2.2 电动悬浮（EDS）理论

2.2.1 涡流悬浮力

当超导磁体在导电轨道上运动时，感应涡流产生排斥力：

$$F_{lift}=\frac{B_0^{2}ab}{{\mu }_{0}h}\cdot f(v)$$

其中：

• $B_0$：超导磁体磁场强度 (T)
• $a,b$：磁体尺寸 (m)
• $h$：悬浮高度 (m)
• $f(v)$：速度相关函数

速度特性：

$$f(v)=\frac{v}{v+v_c}$$

其中 $v_c$ 是特征速度，与轨道电导率和磁体参数有关。

关键特性：

• 低速时悬浮力很小
• 速度超过100 km/h后悬浮力显著增加
• 高速时趋于饱和

2.2.2 导向力原理

EDS系统具有天然的导向稳定性：

$$F_{guidance}=-k_y\cdot y$$

当列车横向偏移时：

• 靠近侧轨道涡流增强，排斥力增大
• 远离侧轨道涡流减弱，排斥力减小
• 产生恢复力使列车回到中心位置

2.2.3 磁阻力（Drag Force）

涡流同时产生与运动方向相反的阻力：

$$F_{drag}=\frac{P_{eddy}}{v}$$

其中 $P_{eddy}$ 是涡流损耗功率。

磁阻特性：

• 低速时磁阻较大
• 高速时磁阻相对减小
• 是EDS系统的主要能量损耗

2.3 直线电机理论

2.3.1 直线同步电机（LSM）

工作原理：

• 车载超导磁体作为励磁
• 地面线圈作为电枢
• 通过控制电枢电流频率产生行波磁场

推力公式：

$$F=\frac{3\pi }{\tau }\cdot \frac{E_f{V}_s}{{\omega }_s{X}_s}\sin \delta$$

其中：

• $\tau$：极距 (m)
• $E_f$：励磁电动势 (V)
• $V_s$：定子电压 (V)
• ${\omega }_s$：同步角频率 (rad/s)
• $X_s$：同步电抗 (Ω)
• $\delta$：功角

2.3.2 直线感应电机（LIM）

工作原理：

• 车载初级绕组产生行波磁场
• 轨道作为次级（短路绕组或导电板）
• 感应电流产生推力

推力-滑差特性：

$$F=F_{max}\cdot \frac{s}{s^2+s_{max}^2}$$

其中：

• $s=(v_s-v)/v_s$：滑差率
• $v_s$：同步速度 (m/s)
• $v$：列车速度 (m/s)
• $s_{max}$：最大推力时的滑差率

2.3.3 效率分析

直线电机效率：

$$\eta =\frac{P_{out}}{P_{in}}=\frac{F\cdot v}{F\cdot v+P_{loss}}$$

损耗组成：

• 铜损：绕组电阻损耗
• 铁损：磁滞损耗和涡流损耗
• 机械损耗：空气阻力等

2.4 涡流制动理论

2.4.1 制动力产生机理

涡流制动利用电磁感应原理：

$$F_{brake}=\frac{B^{2}Lw\sigma dv}{2}$$

其中：

• $B$：制动磁场强度 (T)
• $L,w$：磁铁长度和宽度 (m)
• $\sigma$：轨道电导率 (S/m)
• $d$：轨道厚度 (m)
• $v$：列车速度 (m/s)

速度特性：

• 低速时制动力与速度成正比
• 高速时趋于饱和
• 适合高速制动，低速时需配合摩擦制动

2.4.2 制动功率

涡流制动功率：

$$P_{brake}=F_{brake}\cdot v$$

能量转换：

• 列车动能转化为轨道中的涡流热能
• 无机械磨损
• 制动能量不能回收（与再生制动不同）

2.4.3 制动距离计算

匀减速制动距离：

$$s=\frac{v_0^2}{2a}=\frac{m{v}_0^2}{2F_{brake}}$$

其中 $v_0$ 是初始速度，$a$ 是减速度。

---

三、环境准备

3.1 依赖库安装

本教程需要以下Python库：

pip install numpy matplotlib scipy

各库的作用：

• NumPy：科学计算基础库，用于数组操作和数值计算
• Matplotlib：数据可视化库，用于绘制静态图表
• SciPy：科学计算库，提供信号处理等高级功能

3.2 验证安装

运行以下Python代码验证安装：

import numpy as np
import matplotlib.pyplot as plt

print(f"NumPy版本: {np.__version__}")
print("所有依赖库安装成功！")

3.3 推荐开发环境

• Python版本：3.8及以上
• IDE推荐：VS Code、PyCharm、Jupyter Notebook
• 操作系统：Windows、Linux、macOS均可

---

四、完整代码实现

以下是完整的磁浮列车系统仿真代码，包含5个案例：

"""
磁浮列车系统仿真
=================
本代码模拟磁浮列车的核心磁系统，包括：
1. 电磁悬浮(EMS)系统磁力计算
2. 电动悬浮(EDS)系统涡流力分析
3. 直线同步电机推力计算
4. 导向系统磁场分布
5. 涡流制动特性分析
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle, FancyBboxPatch, Circle, FancyArrowPatch
from matplotlib.collections import PatchCollection
import matplotlib.animation as animation
from PIL import Image
import os

# 强制使用Agg后端，避免GUI窗口
plt.switch_backend('Agg')

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

# 创建输出目录
output_dir = os.path.dirname(os.path.abspath(__file__))

# ==================== 物理常数 ====================
MU_0 = 4 * np.pi * 1e-7  # 真空磁导率 (H/m)
G = 9.81                 # 重力加速度 (m/s²)

# ==================== 案例1: 电磁悬浮(EMS)系统 ====================
def ems_suspension_analysis():
    """
    电磁悬浮系统分析
    模拟常导磁吸型磁浮列车的悬浮力特性
    """
    print("=" * 60)
    print("案例1: 电磁悬浮(EMS)系统分析")
    print("=" * 60)
    
    # 电磁铁参数
    N = 300  # 匝数
    A = 0.1 * 0.5  # 极面积 (m²)，宽0.1m，长0.5m
    I_max = 50  # 最大电流 (A)
    
    # 气隙范围
    gaps = np.linspace(0.005, 0.03, 100)  # 5mm 到 30mm
    
    # 计算不同电流下的悬浮力
    currents = [10, 20, 30, 40, 50]  # A
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 12))
    
    # 图1: 悬浮力 vs 气隙
    ax1 = axes[0, 0]
    colors = plt.cm.viridis(np.linspace(0, 1, len(currents)))
    
    for i, I in enumerate(currents):
        # 电磁吸力公式: F = (μ₀ * N² * I² * A) / (2 * gap²)
        F = (MU_0 * N**2 * I**2 * A) / (2 * gaps**2)
        ax1.plot(gaps * 1000, F / 1000, color=colors[i], linewidth=2, 
                label=f'I = {I} A')
    
    ax1.set_xlabel('气隙 (mm)', fontsize=11)
    ax1.set_ylabel('悬浮力 (kN)', fontsize=11)
    ax1.set_title('电磁悬浮力 vs 气隙', fontsize=12, fontweight='bold')
    ax1.legend(loc='upper right')
    ax1.grid(True, alpha=0.3)
    ax1.set_yscale('log')
    
    # 图2: 悬浮力 vs 电流
    ax2 = axes[0, 1]
    gap_values = [5, 10, 15, 20, 25]  # mm
    I_range = np.linspace(5, 50, 100)
    
    colors2 = plt.cm.plasma(np.linspace(0, 1, len(gap_values)))
    for i, gap_mm in enumerate(gap_values):
        gap = gap_mm / 1000
        F = (MU_0 * N**2 * I_range**2 * A) / (2 * gap**2)
        ax2.plot(I_range, F / 1000, color=colors2[i], linewidth=2,
                label=f'gap = {gap_mm} mm')
    
    ax2.set_xlabel('电流 (A)', fontsize=11)
    ax2.set_ylabel('悬浮力 (kN)', fontsize=11)
    ax2.set_title('电磁悬浮力 vs 电流', fontsize=12, fontweight='bold')
    ax2.legend(loc='upper left')
    ax2.grid(True, alpha=0.3)
    
    # 图3: 悬浮系统刚度分析
    ax3 = axes[1, 0]
    I_operating = 30  # 工作电流
    gap_nominal = 0.010  # 额定气隙 10mm
    
    # 在额定工作点附近计算刚度
    gaps_fine = np.linspace(0.008, 0.012, 200)
    F_fine = (MU_0 * N**2 * I_operating**2 * A) / (2 * gaps_fine**2)
    
    # 计算刚度 (负刚度，需要主动控制)
    stiffness = -np.gradient(F_fine, gaps_fine) / 1000  # kN/m
    
    ax3.plot(gaps_fine * 1000, stiffness, 'b-', linewidth=2)
    ax3.axvline(x=gap_nominal * 1000, color='r', linestyle='--', 
               label=f'额定气隙 = {gap_nominal*1000} mm')
    ax3.set_xlabel('气隙 (mm)', fontsize=11)
    ax3.set_ylabel('负刚度 (kN/m)', fontsize=11)
    ax3.set_title('悬浮系统刚度特性', fontsize=12, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # 图4: EMS系统示意图
    ax4 = axes[1, 1]
    ax4.set_xlim(0, 10)
    ax4.set_ylim(0, 10)
    ax4.set_aspect('equal')
    ax4.axis('off')
    ax4.set_title('EMS磁浮系统结构示意图', fontsize=12, fontweight='bold')
    
    # 轨道
    track = Rectangle((1, 1), 8, 1, facecolor='gray', edgecolor='black', linewidth=2)
    ax4.add_patch(track)
    ax4.text(5, 1.5, '轨道 (F轨)', ha='center', fontsize=10)
    
    # 电磁铁
    magnet = Rectangle((2, 2.5), 6, 1.5, facecolor='lightblue', 
                       edgecolor='blue', linewidth=2)
    ax4.add_patch(magnet)
    ax4.text(5, 3.25, '电磁铁', ha='center', fontsize=10)
    
    # 气隙标注
    ax4.annotate('', xy=(1.5, 2.5), xytext=(1.5, 2),
                arrowprops=dict(arrowstyle='<->', color='red', lw=2))
    ax4.text(1.2, 2.25, '气隙', fontsize=9, color='red')
    
    # 悬浮力箭头
    ax4.annotate('', xy=(5, 4.2), xytext=(5, 5.5),
                arrowprops=dict(arrowstyle='->', color='green', lw=3))
    ax4.text(5.5, 5, '悬浮力', fontsize=10, color='green')
    
    # 重力箭头
    ax4.annotate('', xy=(5, 2.3), xytext=(5, 1),
                arrowprops=dict(arrowstyle='->', color='orange', lw=3))
    ax4.text(5.5, 1.5, '重力', fontsize=10, color='orange')
    
    plt.tight_layout()
    plt.savefig(os.path.join(output_dir, 'fig1_ems_suspension.png'), dpi=150, bbox_inches='tight')
    plt.close()
    
    print(f"✓ EMS悬浮分析完成")
    print(f"  额定气隙: {gap_nominal*1000} mm")
    print(f"  额定电流: {I_operating} A")
    print(f"  悬浮力: {(MU_0 * N**2 * I_operating**2 * A) / (2 * gap_nominal**2) / 1000:.2f} kN")
    
    return gaps, currents

---

五、代码深度解析

5.1 EMS悬浮力计算

# 电磁铁参数
N = 300  # 匝数
A = 0.1 * 0.5  # 极面积 (m²)，宽0.1m，长0.5m
I_max = 50  # 最大电流 (A)

# 气隙范围
gaps = np.linspace(0.005, 0.03, 100)  # 5mm 到 30mm

参数说明：

• 匝数 $N=300$：典型EMS电磁铁匝数
• 极面积 $A=0.05$ m²：单极面积
• 气隙范围：5-30 mm，覆盖工作范围

# 电磁吸力公式: F = (μ₀ * N² * I² * A) / (2 * gap²)
F = (MU_0 * N**2 * I**2 * A) / (2 * gaps**2)

公式推导：

根据麦克斯韦应力张量，电磁吸力为：

$$F=\frac{B^{2}A}{2{\mu }_0}$$

气隙磁感应强度：

$$B=\frac{{\mu }_{0}NI}{2g}$$

代入得：

$$F=\frac{{\mu }_0{N}^2{I}^{2}A}{4g^2}$$

工程意义：

• 悬浮力与电流平方成正比，需要精确控制
• 悬浮力与气隙平方成反比，气隙越小力越大
• 系统具有负刚度，必须采用主动控制

5.2 EDS涡流悬浮力计算

# 超导磁体参数
B0 = 5.0  # 磁体中心磁场 (T)
a = 0.15  # 磁体半宽度 (m)
b = 0.3   # 磁体半长度 (m)

# 导电轨道参数 (铝或铜)
sigma = 3.5e7  # 电导率 (S/m)，铝
d = 0.02  # 轨道厚度 (m)

EDS力模型：

# 简化EDS力模型
# 低速时力较小，高速时趋于饱和
v0 = 20  # 特征速度 (m/s)
F_max = (B0**2 * a * b) / (MU_0 * gap) * 1000  # 最大悬浮力 (N)
F_lift = F_max * (velocities / (velocities + v0))

物理机制：

1. 磁体运动：超导磁体在导电轨道上运动
2. 涡流感应：变化的磁场在轨道中感应涡流
3. 洛伦兹力：涡流与磁场相互作用产生排斥力
4. 速度效应：速度越高，涡流越大，但相位滞后也增大

速度特性：

• 低速时（$v\ll v_c$）：$F\propto v$
• 高速时（$v\gg v_c$）：$F\to F_{max}$
• 特征速度 $v_c$ 与轨道电导率、磁体参数有关

5.3 直线电机推力计算

# 速度范围
velocities = np.linspace(0, 150, 200)  # m/s
v_sync = 100  # 同步速度 (m/s)

# 滑差率
s = (v_sync - velocities) / v_sync

推力-滑差特性：

# 简化推力模型
F_max = 50  # 最大推力 (kN)
F_thrust = F_max * s / (s**2 + 0.1)  # 考虑转差特性

工作区域：

• $s>0$（$v<v_s$）：电动机模式，产生正向推力
• $s<0$（$v>v_s$）：发电机模式，产生制动力（再生制动）
• $s=0$（$v=v_s$）：同步运行，推力为零

控制策略：

• 恒转矩区：保持滑差率恒定，输出最大推力
• 恒功率区：电压恒定，推力随速度降低
• 弱磁区：高速时减弱磁场，扩展速度范围

5.4 涡流制动分析

# 涡流制动力模型
# 低速时与速度成正比，高速时趋于饱和
v_critical = 30  # 临界速度 (m/s)
F_brake_max = 80000  # 最大制动力 (N)

F_brake = F_brake_max * (velocities / (velocities + v_critical))

制动特性：

1. 低速区（$v\ll v_c$）：制动力与速度成正比
2. 高速区（$v\gg v_c$）：制动力趋于饱和
3. 能量转换：动能 $\to$ 涡流热能

与摩擦制动对比：

特性	涡流制动	摩擦制动
速度特性	高速有效	全速域有效
磨损	无	有
热负荷	在轨道上	在制动盘上
能量回收	否	否（再生制动可回收）
适用场景	高速制动	低速制动

---

"""
磁浮列车系统仿真
=================
本代码模拟磁浮列车的核心磁系统，包括：
1. 电磁悬浮(EMS)系统磁力计算
2. 电动悬浮(EDS)系统涡流力分析
3. 直线同步电机推力计算
4. 导向系统磁场分布
5. 涡流制动特性分析
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle, FancyBboxPatch, Circle, FancyArrowPatch
from matplotlib.collections import PatchCollection
import matplotlib.animation as animation
from PIL import Image
import os

# 强制使用Agg后端，避免GUI窗口
plt.switch_backend('Agg')

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

# 创建输出目录
output_dir = os.path.dirname(os.path.abspath(__file__))

# ==================== 物理常数 ====================
MU_0 = 4 * np.pi * 1e-7  # 真空磁导率 (H/m)
G = 9.81                 # 重力加速度 (m/s²)

# ==================== 案例1: 电磁悬浮(EMS)系统 ====================
def ems_suspension_analysis():
    """
    电磁悬浮系统分析
    模拟常导磁吸型磁浮列车的悬浮力特性
    """
    print("=" * 60)
    print("案例1: 电磁悬浮(EMS)系统分析")
    print("=" * 60)
    
    # 电磁铁参数
    N = 300  # 匝数
    A = 0.1 * 0.5  # 极面积 (m²)，宽0.1m，长0.5m
    I_max = 50  # 最大电流 (A)
    
    # 气隙范围
    gaps = np.linspace(0.005, 0.03, 100)  # 5mm 到 30mm
    
    # 计算不同电流下的悬浮力
    currents = [10, 20, 30, 40, 50]  # A
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 12))
    
    # 图1: 悬浮力 vs 气隙
    ax1 = axes[0, 0]
    colors = plt.cm.viridis(np.linspace(0, 1, len(currents)))
    
    for i, I in enumerate(currents):
        # 电磁吸力公式: F = (μ₀ * N² * I² * A) / (2 * gap²)
        F = (MU_0 * N**2 * I**2 * A) / (2 * gaps**2)
        ax1.plot(gaps * 1000, F / 1000, color=colors[i], linewidth=2, 
                label=f'I = {I} A')
    
    ax1.set_xlabel('气隙 (mm)', fontsize=11)
    ax1.set_ylabel('悬浮力 (kN)', fontsize=11)
    ax1.set_title('电磁悬浮力 vs 气隙', fontsize=12, fontweight='bold')
    ax1.legend(loc='upper right')
    ax1.grid(True, alpha=0.3)
    ax1.set_yscale('log')
    
    # 图2: 悬浮力 vs 电流
    ax2 = axes[0, 1]
    gap_values = [5, 10, 15, 20, 25]  # mm
    I_range = np.linspace(5, 50, 100)
    
    colors2 = plt.cm.plasma(np.linspace(0, 1, len(gap_values)))
    for i, gap_mm in enumerate(gap_values):
        gap = gap_mm / 1000
        F = (MU_0 * N**2 * I_range**2 * A) / (2 * gap**2)
        ax2.plot(I_range, F / 1000, color=colors2[i], linewidth=2,
                label=f'gap = {gap_mm} mm')
    
    ax2.set_xlabel('电流 (A)', fontsize=11)
    ax2.set_ylabel('悬浮力 (kN)', fontsize=11)
    ax2.set_title('电磁悬浮力 vs 电流', fontsize=12, fontweight='bold')
    ax2.legend(loc='upper left')
    ax2.grid(True, alpha=0.3)
    
    # 图3: 悬浮系统刚度分析
    ax3 = axes[1, 0]
    I_operating = 30  # 工作电流
    gap_nominal = 0.010  # 额定气隙 10mm
    
    # 在额定工作点附近计算刚度
    gaps_fine = np.linspace(0.008, 0.012, 200)
    F_fine = (MU_0 * N**2 * I_operating**2 * A) / (2 * gaps_fine**2)
    
    # 计算刚度 (负刚度，需要主动控制)
    stiffness = -np.gradient(F_fine, gaps_fine) / 1000  # kN/m
    
    ax3.plot(gaps_fine * 1000, stiffness, 'b-', linewidth=2)
    ax3.axvline(x=gap_nominal * 1000, color='r', linestyle='--', 
               label=f'额定气隙 = {gap_nominal*1000} mm')
    ax3.set_xlabel('气隙 (mm)', fontsize=11)
    ax3.set_ylabel('负刚度 (kN/m)', fontsize=11)
    ax3.set_title('悬浮系统刚度特性', fontsize=12, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # 图4: EMS系统示意图
    ax4 = axes[1, 1]
    ax4.set_xlim(0, 10)
    ax4.set_ylim(0, 10)
    ax4.set_aspect('equal')
    ax4.axis('off')
    ax4.set_title('EMS磁浮系统结构示意图', fontsize=12, fontweight='bold')
    
    # 轨道
    track = Rectangle((1, 1), 8, 1, facecolor='gray', edgecolor='black', linewidth=2)
    ax4.add_patch(track)
    ax4.text(5, 1.5, '轨道 (F轨)', ha='center', fontsize=10)
    
    # 电磁铁
    magnet = Rectangle((2, 2.5), 6, 1.5, facecolor='lightblue', 
                       edgecolor='blue', linewidth=2)
    ax4.add_patch(magnet)
    ax4.text(5, 3.25, '电磁铁', ha='center', fontsize=10)
    
    # 气隙标注
    ax4.annotate('', xy=(1.5, 2.5), xytext=(1.5, 2),
                arrowprops=dict(arrowstyle='<->', color='red', lw=2))
    ax4.text(1.2, 2.25, '气隙', fontsize=9, color='red')
    
    # 悬浮力箭头
    ax4.annotate('', xy=(5, 4.2), xytext=(5, 5.5),
                arrowprops=dict(arrowstyle='->', color='green', lw=3))
    ax4.text(5.5, 5, '悬浮力', fontsize=10, color='green')
    
    # 重力箭头
    ax4.annotate('', xy=(5, 2.3), xytext=(5, 1),
                arrowprops=dict(arrowstyle='->', color='orange', lw=3))
    ax4.text(5.5, 1.5, '重力', fontsize=10, color='orange')
    
    plt.tight_layout()
    plt.savefig(os.path.join(output_dir, 'fig1_ems_suspension.png'), dpi=150, bbox_inches='tight')
    plt.close()
    
    print(f"✓ EMS悬浮分析完成")
    print(f"  额定气隙: {gap_nominal*1000} mm")
    print(f"  额定电流: {I_operating} A")
    print(f"  悬浮力: {(MU_0 * N**2 * I_operating**2 * A) / (2 * gap_nominal**2) / 1000:.2f} kN")
    
    return gaps, currents

# ==================== 案例2: 电动悬浮(EDS)系统 ====================
def eds_suspension_analysis():
    """
    电动悬浮系统分析
    模拟超导磁浮列车的涡流悬浮力特性
    """
    print("\n" + "=" * 60)
    print("案例2: 电动悬浮(EDS)系统分析")
    print("=" * 60)
    
    # 超导磁体参数
    B0 = 5.0  # 磁体中心磁场 (T)
    a = 0.15  # 磁体半宽度 (m)
    b = 0.3   # 磁体半长度 (m)
    
    # 导电轨道参数 (铝或铜)
    sigma = 3.5e7  # 电导率 (S/m)，铝
    d = 0.02  # 轨道厚度 (m)
    
    # 速度范围
    velocities = np.linspace(0, 150, 200)  # 0 到 150 m/s (540 km/h)
    
    # 气隙范围
    gaps = np.linspace(0.05, 0.3, 100)  # 5cm 到 30cm
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 12))
    
    # 图1: 悬浮力 vs 速度
    ax1 = axes[0, 0]
    gap_values = [0.05, 0.10, 0.15, 0.20]  # m
    colors = plt.cm.viridis(np.linspace(0, 1, len(gap_values)))
    
    for i, gap in enumerate(gap_values):
        # 简化EDS力模型
        # 低速时力较小，高速时趋于饱和
        v0 = 20  # 特征速度 (m/s)
        F_max = (B0**2 * a * b) / (MU_0 * gap) * 1000  # 最大悬浮力 (N)
        F_lift = F_max * (velocities / (velocities + v0))
        ax1.plot(velocities * 3.6, F_lift / 1000, color=colors[i], linewidth=2,
                label=f'gap = {gap*100:.0f} cm')
    
    ax1.set_xlabel('速度 (km/h)', fontsize=11)
    ax1.set_ylabel('悬浮力 (kN)', fontsize=11)
    ax1.set_title('EDS悬浮力 vs 速度', fontsize=12, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    ax1.axvline(x=100, color='red', linestyle='--', alpha=0.5, label='起飞速度')
    
    # 图2: 悬浮力 vs 气隙
    ax2 = axes[0, 1]
    v_values = [50, 100, 150]  # m/s
    colors2 = plt.cm.plasma(np.linspace(0, 1, len(v_values)))
    
    for i, v in enumerate(v_values):
        F_lift = (B0**2 * a * b) / (MU_0 * gaps) * 1000 * (v / (v + 20))
        ax2.plot(gaps * 100, F_lift / 1000, color=colors2[i], linewidth=2,
                label=f'v = {v*3.6:.0f} km/h')
    
    ax2.set_xlabel('气隙 (cm)', fontsize=11)
    ax2.set_ylabel('悬浮力 (kN)', fontsize=11)
    ax2.set_title('EDS悬浮力 vs 气隙', fontsize=12, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # 图3: 导向力分析
    ax3 = axes[1, 0]
    # 横向偏移
    offsets = np.linspace(-0.1, 0.1, 100)  # -10cm 到 10cm
    v_cruise = 120  # 巡航速度 (m/s)
    gap_nominal = 0.15  # 额定气隙 (m)
    
    # 导向力与偏移成正比（恢复力）
    k_guidance = 50000  # 导向刚度 (N/m)
    F_guidance = -k_guidance * offsets
    
    ax3.plot(offsets * 100, F_guidance / 1000, 'b-', linewidth=2)
    ax3.axvline(x=0, color='red', linestyle='--', alpha=0.5)
    ax3.axhline(y=0, color='red', linestyle='--', alpha=0.5)
    ax3.set_xlabel('横向偏移 (cm)', fontsize=11)
    ax3.set_ylabel('导向力 (kN)', fontsize=11)
    ax3.set_title('EDS导向力特性', fontsize=12, fontweight='bold')
    ax3.grid(True, alpha=0.3)
    ax3.fill_between(offsets * 100, 0, F_guidance / 1000, 
                     where=(F_guidance > 0), alpha=0.3, color='green', label='恢复力')
    ax3.fill_between(offsets * 100, 0, F_guidance / 1000, 
                     where=(F_guidance < 0), alpha=0.3, color='red', label='恢复力')
    ax3.legend()
    
    # 图4: EDS系统示意图
    ax4 = axes[1, 1]
    ax4.set_xlim(0, 10)
    ax4.set_ylim(0, 10)
    ax4.set_aspect('equal')
    ax4.axis('off')
    ax4.set_title('EDS磁浮系统结构示意图', fontsize=12, fontweight='bold')
    
    # 轨道
    track_l = Rectangle((1, 4), 3.5, 0.5, facecolor='gray', edgecolor='black', linewidth=2)
    track_r = Rectangle((5.5, 4), 3.5, 0.5, facecolor='gray', edgecolor='black', linewidth=2)
    ax4.add_patch(track_l)
    ax4.add_patch(track_r)
    ax4.text(2.75, 4.25, '轨道', ha='center', fontsize=9)
    ax4.text(7.25, 4.25, '轨道', ha='center', fontsize=9)
    
    # 超导磁体
    magnet_l = Rectangle((2, 5), 1.5, 1, facecolor='lightblue', 
                         edgecolor='blue', linewidth=2)
    magnet_r = Rectangle((6.5, 5), 1.5, 1, facecolor='lightblue', 
                         edgecolor='blue', linewidth=2)
    ax4.add_patch(magnet_l)
    ax4.add_patch(magnet_r)
    ax4.text(2.75, 5.5, '超导\n磁体', ha='center', fontsize=8)
    ax4.text(7.25, 5.5, '超导\n磁体', ha='center', fontsize=8)
    
    # 车体
    vehicle = Rectangle((1.5, 6.2), 7, 2.5, facecolor='lightyellow', 
                        edgecolor='orange', linewidth=2)
    ax4.add_patch(vehicle)
    ax4.text(5, 7.5, '车体', ha='center', fontsize=10)
    
    # 气隙标注
    ax4.annotate('', xy=(4.5, 4.5), xytext=(4.5, 5),
                arrowprops=dict(arrowstyle='<->', color='red', lw=2))
    ax4.text(4.6, 4.75, '气隙', fontsize=9, color='red')
    
    # 悬浮力箭头
    ax4.annotate('', xy=(2.75, 6.2), xytext=(2.75, 7.5),
                arrowprops=dict(arrowstyle='->', color='green', lw=3))
    ax4.text(3.2, 7, '悬浮力', fontsize=9, color='green')
    
    plt.tight_layout()
    plt.savefig(os.path.join(output_dir, 'fig2_eds_suspension.png'), dpi=150, bbox_inches='tight')
    plt.close()
    
    print(f"✓ EDS悬浮分析完成")
    print(f"  超导磁场: {B0} T")
    print(f"  起飞速度: ~100 km/h")
    print(f"  巡航速度: {v_cruise*3.6:.0f} km/h")
    
    return velocities, gaps

# ==================== 案例3: 直线同步电机(LIM)推力分析 ====================
def linear_motor_analysis():
    """
    直线同步电机分析
    计算磁浮列车推进系统的推力特性
    """
    print("\n" + "=" * 60)
    print("案例3: 直线同步电机(LIM)推力分析")
    print("=" * 60)
    
    # 电机参数
    tau = 0.5  # 极距 (m)
    L = 10     # 电机长度 (m)
    I_s = 500  # 电枢电流 (A)
    N_p = 20   # 极对数
    
    # 速度范围
    velocities = np.linspace(0, 150, 200)  # m/s
    v_sync = 100  # 同步速度 (m/s)
    
    # 滑差率
    s = (v_sync - velocities) / v_sync
    s[velocities > v_sync] = (velocities[velocities > v_sync] - v_sync) / velocities[velocities > v_sync]
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 12))
    
    # 图1: 推力 vs 速度
    ax1 = axes[0, 0]
    
    # 简化推力模型
    F_max = 50  # 最大推力 (kN)
    F_thrust = F_max * s / (s**2 + 0.1)  # 考虑转差特性
    F_thrust[velocities > v_sync] = -F_thrust[velocities > v_sync]  # 再生制动区
    
    ax1.plot(velocities * 3.6, F_thrust, 'b-', linewidth=2)
    ax1.axvline(x=v_sync * 3.6, color='red', linestyle='--', 
               label=f'同步速度 = {v_sync*3.6:.0f} km/h')
    ax1.axhline(y=0, color='black', linestyle='-', alpha=0.3)
    ax1.fill_between(velocities * 3.6, 0, F_thrust, 
                     where=(F_thrust > 0), alpha=0.3, color='green', label='牵引区')
    ax1.fill_between(velocities * 3.6, 0, F_thrust, 
                     where=(F_thrust < 0), alpha=0.3, color='red', label='制动区')
    ax1.set_xlabel('速度 (km/h)', fontsize=11)
    ax1.set_ylabel('推力 (kN)', fontsize=11)
    ax1.set_title('直线电机推力 vs 速度', fontsize=12, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 图2: 推力 vs 滑差率
    ax2 = axes[0, 1]
    s_range = np.linspace(-1, 1, 200)
    F_thrust_s = F_max * s_range / (s_range**2 + 0.1)
    
    ax2.plot(s_range, F_thrust_s, 'b-', linewidth=2)
    ax2.axvline(x=0, color='red', linestyle='--', alpha=0.5)
    ax2.axhline(y=0, color='black', linestyle='-', alpha=0.3)
    ax2.set_xlabel('滑差率', fontsize=11)
    ax2.set_ylabel('推力 (kN)', fontsize=11)
    ax2.set_title('直线电机推力 vs 滑差率', fontsize=12, fontweight='bold')
    ax2.grid(True, alpha=0.3)
    
    # 图3: 效率分析
    ax3 = axes[1, 0]
    # 简化效率模型
    eta = np.zeros_like(velocities)
    eta[velocities <= v_sync] = 0.85 * (velocities[velocities <= v_sync] / v_sync)
    eta[velocities > v_sync] = 0.85 * (1 - 0.3 * (velocities[velocities > v_sync] - v_sync) / v_sync)
    eta = np.clip(eta, 0, 0.95)
    
    ax3.plot(velocities * 3.6, eta * 100, 'g-', linewidth=2)
    ax3.axvline(x=v_sync * 3.6, color='red', linestyle='--', alpha=0.5)
    ax3.set_xlabel('速度 (km/h)', fontsize=11)
    ax3.set_ylabel('效率 (%)', fontsize=11)
    ax3.set_title('直线电机效率 vs 速度', fontsize=12, fontweight='bold')
    ax3.grid(True, alpha=0.3)
    ax3.set_ylim(0, 100)
    
    # 图4: 功率分析
    ax4 = axes[1, 1]
    P_mech = F_thrust * velocities / 1000  # 机械功率 (MW)
    P_elec = np.where(P_mech > 0, P_mech / np.maximum(eta, 0.1), P_mech * np.maximum(eta, 0.1))
    
    ax4.plot(velocities * 3.6, P_mech, 'b-', linewidth=2, label='机械功率')
    ax4.plot(velocities * 3.6, P_elec, 'r--', linewidth=2, label='电功率')
    ax4.axvline(x=v_sync * 3.6, color='red', linestyle=':', alpha=0.5)
    ax4.axhline(y=0, color='black', linestyle='-', alpha=0.3)
    ax4.set_xlabel('速度 (km/h)', fontsize=11)
    ax4.set_ylabel('功率 (MW)', fontsize=11)
    ax4.set_title('直线电机功率 vs 速度', fontsize=12, fontweight='bold')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig(os.path.join(output_dir, 'fig3_linear_motor.png'), dpi=150, bbox_inches='tight')
    plt.close()
    
    print(f"✓ 直线电机分析完成")
    print(f"  同步速度: {v_sync*3.6:.0f} km/h")
    print(f"  最大推力: {F_max:.1f} kN")
    print(f"  峰值功率: {np.max(np.abs(P_mech)):.2f} MW")
    
    return velocities, F_thrust

# ==================== 案例4: 磁浮列车动力学仿真 ====================
def maglev_dynamics_simulation():
    """
    磁浮列车动力学仿真
    模拟列车的加速、巡航和制动过程
    """
    print("\n" + "=" * 60)
    print("案例4: 磁浮列车动力学仿真")
    print("=" * 60)
    
    # 列车参数
    m = 50000  # 质量 (kg)
    v_max = 150  # 最大速度 (m/s) = 540 km/h
    F_max = 50000  # 最大推力 (N)
    
    # 仿真参数
    dt = 0.1  # 时间步长 (s)
    t_total = 300  # 总仿真时间 (s)
    t = np.arange(0, t_total, dt)
    n_steps = len(t)
    
    # 初始化状态
    v = np.zeros(n_steps)  # 速度
    x = np.zeros(n_steps)  # 位置
    a = np.zeros(n_steps)  # 加速度
    F = np.zeros(n_steps)  # 推力
    
    # 运行工况
    # 0-60s: 加速
    # 60-180s: 巡航
    # 180-240s: 减速
    # 240-300s: 停车
    
    for i in range(1, n_steps):
        current_t = t[i]
        
        # 确定当前工况
        if current_t < 60:
            # 加速阶段
            target_F = F_max
        elif current_t < 180:
            # 巡航阶段
            target_F = 0
        elif current_t < 240:
            # 减速阶段
            target_F = -F_max * 0.8
        else:
            # 停车
            target_F = 0
            if v[i-1] < 1:
                v[i-1] = 0
        
        # 考虑空气阻力 (简化模型)
        F_drag = 0.5 * 1.225 * 10 * v[i-1]**2  # 简化的空气阻力
        
        # 计算净推力
        F[i] = target_F - F_drag
        
        # 计算加速度
        a[i] = F[i] / m
        
        # 更新速度和位置
        v[i] = v[i-1] + a[i] * dt
        v[i] = max(0, v[i])  # 速度不能为负
        
        x[i] = x[i-1] + v[i] * dt
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 图1: 速度-时间曲线
    ax1 = axes[0, 0]
    ax1.plot(t, v * 3.6, 'b-', linewidth=2)
    ax1.axhline(y=v_max * 3.6, color='red', linestyle='--', alpha=0.5, label='最大速度')
    ax1.axvline(x=60, color='green', linestyle=':', alpha=0.5)
    ax1.axvline(x=180, color='green', linestyle=':', alpha=0.5)
    ax1.axvline(x=240, color='orange', linestyle=':', alpha=0.5)
    ax1.set_xlabel('时间 (s)', fontsize=11)
    ax1.set_ylabel('速度 (km/h)', fontsize=11)
    ax1.set_title('列车速度-时间曲线', fontsize=12, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    ax1.text(30, 100, '加速', fontsize=10, ha='center')
    ax1.text(120, 100, '巡航', fontsize=10, ha='center')
    ax1.text(210, 100, '减速', fontsize=10, ha='center')
    ax1.text(270, 50, '停车', fontsize=10, ha='center')
    
    # 图2: 位置-时间曲线
    ax2 = axes[0, 1]
    ax2.plot(t, x / 1000, 'g-', linewidth=2)
    ax2.set_xlabel('时间 (s)', fontsize=11)
    ax2.set_ylabel('位置 (km)', fontsize=11)
    ax2.set_title('列车位置-时间曲线', fontsize=12, fontweight='bold')
    ax2.grid(True, alpha=0.3)
    
    # 图3: 加速度-时间曲线
    ax3 = axes[1, 0]
    ax3.plot(t, a / G, 'r-', linewidth=2)
    ax3.axhline(y=0, color='black', linestyle='-', alpha=0.3)
    ax3.set_xlabel('时间 (s)', fontsize=11)
    ax3.set_ylabel('加速度 (g)', fontsize=11)
    ax3.set_title('列车加速度-时间曲线', fontsize=12, fontweight='bold')
    ax3.grid(True, alpha=0.3)
    
    # 图4: 推力-速度曲线
    ax4 = axes[1, 1]
    ax4.plot(v * 3.6, F / 1000, 'purple', linewidth=2)
    ax4.axhline(y=0, color='black', linestyle='-', alpha=0.3)
    ax4.set_xlabel('速度 (km/h)', fontsize=11)
    ax4.set_ylabel('推力 (kN)', fontsize=11)
    ax4.set_title('列车推力-速度曲线', fontsize=12, fontweight='bold')
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig(os.path.join(output_dir, 'fig4_dynamics.png'), dpi=150, bbox_inches='tight')
    plt.close()
    
    print(f"✓ 动力学仿真完成")
    print(f"  最大速度: {np.max(v)*3.6:.1f} km/h")
    print(f"  总行程: {x[-1]/1000:.2f} km")
    print(f"  最大加速度: {np.max(a)/G:.3f} g")
    
    return t, v, x, a

# ==================== 案例5: 涡流制动分析 ====================
def eddy_current_braking():
    """
    涡流制动分析
    计算磁浮列车涡流制动系统的制动力特性
    """
    print("\n" + "=" * 60)
    print("案例5: 涡流制动分析")
    print("=" * 60)
    
    # 制动磁铁参数
    B_brake = 1.5  # 制动磁场 (T)
    L_magnet = 2.0  # 磁铁长度 (m)
    w_magnet = 0.3  # 磁铁宽度 (m)
    
    # 轨道参数
    sigma = 3.5e7  # 铝轨道电导率 (S/m)
    d_track = 0.02  # 轨道厚度 (m)
    
    # 速度范围
    velocities = np.linspace(0, 150, 200)  # m/s
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 12))
    
    # 图1: 制动力 vs 速度
    ax1 = axes[0, 0]
    
    # 涡流制动力模型
    # 低速时与速度成正比，高速时趋于饱和
    v_critical = 30  # 临界速度 (m/s)
    F_brake_max = 80000  # 最大制动力 (N)
    
    F_brake = F_brake_max * (velocities / (velocities + v_critical))
    
    ax1.plot(velocities * 3.6, F_brake / 1000, 'b-', linewidth=2)
    ax1.axhline(y=F_brake_max / 1000, color='red', linestyle='--', 
               label=f'最大制动力 = {F_brake_max/1000:.1f} kN')
    ax1.set_xlabel('速度 (km/h)', fontsize=11)
    ax1.set_ylabel('制动力 (kN)', fontsize=11)
    ax1.set_title('涡流制动力 vs 速度', fontsize=12, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 图2: 制动功率 vs 速度
    ax2 = axes[0, 1]
    P_brake = F_brake * velocities / 1000  # 制动功率 (kW)
    
    ax2.plot(velocities * 3.6, P_brake / 1000, 'r-', linewidth=2)
    ax2.set_xlabel('速度 (km/h)', fontsize=11)
    ax2.set_ylabel('制动功率 (MW)', fontsize=11)
    ax2.set_title('涡流制动功率 vs 速度', fontsize=12, fontweight='bold')
    ax2.grid(True, alpha=0.3)
    
    # 图3: 制动距离分析
    ax3 = axes[1, 0]
    m = 50000  # 列车质量 (kg)
    v_initial = np.linspace(50, 150, 100)  # 初始速度 (m/s)
    
    # 简化制动距离计算
    # 假设平均减速度
    a_brake = F_brake_max / m  # 最大减速度
    braking_distance = v_initial**2 / (2 * a_brake) / 1000  # km
    
    ax3.plot(v_initial * 3.6, braking_distance, 'g-', linewidth=2)
    ax3.set_xlabel('初始速度 (km/h)', fontsize=11)
    ax3.set_ylabel('制动距离 (km)', fontsize=11)
    ax3.set_title('涡流制动距离 vs 初始速度', fontsize=12, fontweight='bold')
    ax3.grid(True, alpha=0.3)
    
    # 图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')
    
    # 轨道
    track = Rectangle((1, 4), 8, 0.5, facecolor='gray', edgecolor='black', linewidth=2)
    ax4.add_patch(track)
    ax4.text(5, 4.25, '导电轨道 (铝)', ha='center', fontsize=10)
    
    # 制动磁铁
    magnet = Rectangle((2, 5), 6, 1, facecolor='lightcoral', 
                       edgecolor='red', linewidth=2)
    ax4.add_patch(magnet)
    ax4.text(5, 5.5, '制动磁铁', ha='center', fontsize=10)
    
    # 涡流示意
    for i in range(5):
        circle = Circle((2.5 + i*1.2, 4.25), 0.15, facecolor='yellow', 
                       edgecolor='orange', linewidth=1)
        ax4.add_patch(circle)
    ax4.text(5, 3.5, '感应涡流', ha='center', fontsize=9, color='orange')
    
    # 运动方向
    ax4.annotate('', xy=(8.5, 5.5), xytext=(1.5, 5.5),
                arrowprops=dict(arrowstyle='->', color='blue', lw=3))
    ax4.text(5, 6.5, '运动方向', ha='center', fontsize=10, color='blue')
    
    # 制动力箭头
    ax4.annotate('', xy=(1.5, 5.5), xytext=(3, 5.5),
                arrowprops=dict(arrowstyle='->', color='red', lw=3))
    ax4.text(1, 5.5, '制动力', fontsize=10, color='red')
    
    plt.tight_layout()
    plt.savefig(os.path.join(output_dir, 'fig5_braking.png'), dpi=150, bbox_inches='tight')
    plt.close()
    
    print(f"✓ 涡流制动分析完成")
    print(f"  最大制动力: {F_brake_max/1000:.1f} kN")
    print(f"  最大制动功率: {np.max(P_brake)/1000:.2f} MW")
    print(f"  从500km/h制动距离: ~{braking_distance[-1]:.1f} km")
    
    return velocities, F_brake

# ==================== 主程序 ====================
if __name__ == "__main__":
    print("\n" + "=" * 70)
    print("磁浮列车系统仿真")
    print("=" * 70)
    
    # 运行所有案例
    ems_suspension_analysis()
    eds_suspension_analysis()
    linear_motor_analysis()
    maglev_dynamics_simulation()
    eddy_current_braking()
    
    print("\n" + "=" * 70)
    print("所有仿真案例已完成！")
    print("=" * 70)
    print("\n生成的文件:")
    print("  - fig1_ems_suspension.png: EMS电磁悬浮分析")
    print("  - fig2_eds_suspension.png: EDS电动悬浮分析")
    print("  - fig3_linear_motor.png: 直线电机推力分析")
    print("  - fig4_dynamics.png: 列车动力学仿真")
    print("  - fig5_braking.png: 涡流制动分析")

---