第二十五篇：电-磁-热耦合分析

摘要

电-磁-热耦合是多物理场耦合仿真中的高级研究领域，涉及电磁场、温度场和结构变形之间的复杂相互作用。本主题系统阐述电-磁-热耦合的基本物理机制，包括电磁感应加热、涡流效应、磁致热效应等核心概念。详细推导电磁场控制方程（Maxwell方程组）、热传导方程和结构力学方程的耦合形式，建立三维多物理场耦合数学模型。介绍感应加热、变压器热分析、电磁热治疗等典型工程应用。通过六个工程案例的Python仿真实现，展示电-磁-热耦合问题的数值求解方法，包括时谐电磁场分析、瞬态涡流计算、多物理场顺序耦合等技术。最后探讨电-磁-热耦合在工业加热、电力设备、生物医学等前沿领域的应用前景。

关键词

电-磁-热耦合，感应加热，涡流损耗，Maxwell方程组，电磁热治疗，变压器热分析，多物理场耦合，有限元方法

---

1. 引言

1.1 电-磁-热耦合现象概述

电-磁-热耦合是指电磁场、温度场和结构场三者之间的相互作用现象，是现代工程应用中最为复杂的多物理场耦合问题之一。当导体处于交变电磁场中时，会同时产生多种物理效应：电磁感应产生涡流，涡流产生焦耳热；温度变化改变材料电磁特性；热应力引起结构变形，进而影响电磁场分布。这种多物理场的强耦合机制在感应加热、电力设备、电磁冶金等领域具有核心重要性。

电-磁-热耦合问题的复杂性主要体现在以下几个方面：

多物理场相互作用：电磁场通过涡流损耗和磁滞损耗产生热源，温度场通过材料特性（电导率、磁导率）影响电磁场分布，热应力导致结构变形改变几何边界条件。三个物理场形成闭环耦合系统。

非线性特性：铁磁材料的磁导率随磁场强度和温度变化呈现强非线性；导体电导率随温度变化；热辐射边界条件具有非线性特征。这些非线性因素使得问题求解极具挑战性。

多时间尺度差异：电磁场传播时间尺度为纳秒级，热传导时间尺度为秒级至分钟级，结构响应时间尺度介于两者之间。这种时间尺度差异要求采用特殊的时间积分策略。

材料各向异性：电工钢、永磁材料等具有磁各向异性；复合材料可能同时具有电、热、力学各向异性。各向异性材料的本构关系建模增加了问题的复杂度。

1.2 工程应用背景

电-磁-热耦合分析在众多高科技工程领域具有关键应用价值：

感应加热与熔炼：感应加热利用交变磁场在导体中感应涡流产生热量，具有加热速度快、效率高、可控性好等优点，广泛应用于金属热处理、熔炼、焊接等工艺。电-磁-热耦合分析可以优化感应器设计，预测温度分布，控制加热质量。

电力变压器：变压器在运行中因铁芯磁滞损耗、涡流损耗和绕组铜损而产生热量，温升限制是变压器设计的核心约束。电-磁-热耦合分析可以评估热点温度，优化冷却系统，预测绝缘寿命。

电磁热治疗：在肿瘤治疗中，利用交变磁场使植入的磁性纳米颗粒产热，或通过感应加热直接加热病变组织。电-磁-热耦合分析可以优化治疗参数，确保肿瘤区域达到治疗温度同时保护正常组织。

电机与发电机：旋转电机中的电磁损耗导致温升，影响绝缘性能和永磁体特性。电-磁-热耦合分析可以预测温度分布，评估冷却效果，防止永磁体退磁。

磁悬浮系统：磁悬浮列车和轴承中的电磁场与温度场耦合影响系统稳定性和承载能力。电-磁-热耦合分析可以优化磁路设计，控制温升，确保系统可靠运行。

无线充电系统：无线充电过程中，发射和接收线圈的电磁场耦合产生热量，影响充电效率和安全性。电-磁-热耦合分析可以优化线圈设计，控制温升。

1.3 本主题内容安排

本主题将系统介绍电-磁-热耦合分析的理论基础和数值方法，具体安排如下：

第2节阐述电-磁-热耦合的物理机制，包括电磁感应、涡流效应、磁滞损耗等；第3节建立电-磁-热耦合的数学模型，推导Maxwell方程组、热传导方程和结构力学方程的耦合形式；第4节介绍数值求解方法，包括有限元离散、时谐分析、瞬态求解等；第5节通过六个工程案例展示电-磁-热耦合分析的Python实现；第6节讨论结果分析与工程应用；第7节总结全文并展望未来发展方向。

---

2. 电-磁-热耦合物理机制

2.1 电磁场基本理论

2.1.1 Maxwell方程组

电磁场的基本规律由Maxwell方程组描述：

法拉第电磁感应定律：
$$\nabla \times \mathbf{E}=-\frac{\partial \mathbf{B}}{\partial t}$$

安培环路定律（含位移电流）：
$$\nabla \times \mathbf{H}=\mathbf{J}+\frac{\partial \mathbf{D}}{\partial t}$$

高斯电通量定律：
$$\nabla \cdot \mathbf{D}={\rho }_e$$

高斯磁通量定律：
$$\nabla \cdot \mathbf{B}=0$$

其中，$\mathbf{E}$为电场强度（V/m），$\mathbf{H}$为磁场强度（A/m），$\mathbf{B}$为磁感应强度（T），$\mathbf{D}$为电位移矢量（C/m²），$\mathbf{J}$为电流密度（A/m²），${\rho }_e$为电荷密度（C/m³）。

2.1.2 本构关系

电磁场量之间的关系由本构方程描述：

$$\mathbf{D}=\epsilon \mathbf{E}={\epsilon }_0{\epsilon }_r\mathbf{E}$$

$$\mathbf{B}=\mu \mathbf{H}={\mu }_0{\mu }_r\mathbf{H}$$

$$\mathbf{J}=\sigma \mathbf{E}+{\mathbf{J}}_s$$

其中，$\epsilon$为介电常数，$\mu$为磁导率，$\sigma$为电导率，${\mathbf{J}}_s$为源电流密度。

对于铁磁材料，磁导率不是常数，而是磁场强度和温度的函数：

$$\mu =\mu (H,T)$$

这种非线性关系是电-磁-热耦合的重要来源。

2.1.3 电磁势表示

引入电磁势可以简化Maxwell方程组的求解：

磁矢势$\mathbf{A}$：
$$\mathbf{B}=\nabla \times \mathbf{A}$$

电标势$\varphi$：
$$\mathbf{E}=-\nabla \varphi -\frac{\partial \mathbf{A}}{\partial t}$$

在库仑规范$\nabla \cdot \mathbf{A}=0$下，可以得到关于$\mathbf{A}$和$\varphi$的波动方程。

2.2 电磁损耗机制

2.2.1 涡流损耗

当导体处于交变磁场中时，根据法拉第定律会感应出涡流。涡流在导体电阻上产生焦耳热，称为涡流损耗。

对于时谐场（角频率$\omega$），涡流损耗功率密度为：

$$p_{eddy}=\frac{|\mathbf{J}{|}^2}{\sigma }=\sigma |\mathbf{E}{|}^2$$

对于厚度为$d$的薄板，在均匀正弦磁场中的涡流损耗为：

$$P_{eddy}=\frac{{\pi }^2{f}^2{B}_{max}^2{d}^2}{6\rho }V$$

其中，$f$为频率，$B_{max}$为最大磁感应强度，$\rho$为电阻率，$V$为体积。

为减小涡流损耗，电工钢通常采用薄片叠压结构，片间有绝缘层。

2.2.2 磁滞损耗

铁磁材料在交变磁场中磁化时，磁畴壁移动和磁矩转动需要克服阻力，产生能量损耗，称为磁滞损耗。

磁滞损耗功率与磁滞回线面积成正比：

$$p_{hyst}=f\oint HdB$$

其中，$f$为磁场频率。对于正弦激励，常用Steinmetz经验公式：

$$p_{hyst}=k_{h}f{B}_{max}^n$$

其中，$k_h$为磁滞损耗系数，$n$为Steinmetz指数（通常1.6-2.2）。

磁滞损耗与材料特性密切相关，取向硅钢的磁滞损耗远低于无取向硅钢。

2.2.3 剩余损耗

除涡流损耗和磁滞损耗外，还存在由磁畴壁运动引起的剩余损耗：

$$p_{res}=k_r(f{B}_{max}{)}^{3/2}$$

总铁损为三者之和：

$$p_{iron}=p_{hyst}+p_{eddy}+p_{res}$$

2.3 感应加热原理

2.3.1 趋肤效应

交变电流在导体中分布不均匀，电流密度随深度呈指数衰减，称为趋肤效应。趋肤深度定义为电流密度衰减到表面值$1/e$时的深度：

$$\delta =\sqrt{\frac{2}{\omega \mu \sigma }}=\sqrt{\frac{\rho }{\pi f\mu }}$$

典型材料的趋肤深度：

• 铜（20°C，50Hz）：约9.3 mm
• 铝（20°C，50Hz）：约11.7 mm
• 钢（20°C，50Hz）：约1.0 mm
• 钢（20°C，10kHz）：约0.07 mm

高频感应加热利用趋肤效应实现表面加热或深层加热控制。

2.3.2 感应加热功率

感应加热功率密度可以表示为：

$$p_{ind}=\frac{1}{2}\sigma |\mathbf{E}{|}^2=\frac{1}{2}{\omega }^2\sigma |\mathbf{A}{|}^2$$

总加热功率与频率、磁场强度、材料电导率等因素相关。对于圆柱形工件，单位长度加热功率为：

$$P=\frac{\pi \sigma {\omega }^2{r}^2{B}_0^2}{2}\cdot F(k)$$

其中，$r$为工件半径，$B_0$为表面磁场强度，$F(k)$为与工件尺寸和趋肤深度相关的修正函数。

2.3.3 磁-热耦合特性

铁磁材料在居里温度以下具有强磁性和低磁导率，加热效率高；超过居里温度后，材料变为顺磁性，磁导率急剧下降，加热效率显著降低。这种特性可用于感应加热的温度自限制。

2.4 热-结构耦合

2.4.1 热应力产生

温度分布不均匀导致热膨胀差异，产生热应力。热应变与温度变化的关系为：

$${\epsilon }_{thermal}=\alpha (T-T_0)\mathbf{I}$$

其中，$\alpha$为热膨胀系数，$T_0$为参考温度。

总应变为机械应变和热应变之和：

$$\mathbf{\epsilon }={\mathbf{\epsilon }}_{mech}+{\mathbf{\epsilon }}_{thermal}$$

2.4.2 热变形影响

热变形改变几何形状，进而影响电磁场分布。对于小变形，可以采用欧拉描述；对于大变形，需要采用拉格朗日描述或任意拉格朗日-欧拉（ALE）方法。

---

3. 电-磁-热耦合数学模型

3.1 电磁场控制方程

3.1.1 时谐电磁场

对于正弦时变电磁场，采用相量表示：

$$\mathbf{E}(\mathbf{r},t)=\text{Re}[\tilde{\mathbf{E}}(\mathbf{r})e^{j\omega t}]$$

时谐Maxwell方程组为：

$$\nabla \times \tilde{\mathbf{E}}=-j\omega \tilde{\mathbf{B}}$$

$$\nabla \times \tilde{\mathbf{H}}=\tilde{\mathbf{J}}+j\omega \tilde{\mathbf{D}}$$

对于低频问题（位移电流可忽略），得到扩散方程：

$$\nabla \times \left(\frac{1}{\mu }\nabla \times \tilde{\mathbf{A}}\right)+j\omega \sigma \tilde{\mathbf{A}}={\tilde{\mathbf{J}}}_s$$

3.1.2 磁矢势方程

在库仑规范下，磁矢势满足：

$$\nabla \times \left(\frac{1}{\mu }\nabla \times \mathbf{A}\right)+\sigma \frac{\partial \mathbf{A}}{\partial t}={\mathbf{J}}_s$$

对于二维问题（$\mathbf{A}=A_z{\mathbf{e}}_z$），简化为标量方程：

$$-\frac{\partial }{\partial x}\left(\frac{1}{\mu }\frac{\partial A_z}{\partial x}\right)-\frac{\partial }{\partial y}\left(\frac{1}{\mu }\frac{\partial A_z}{\partial y}\right)+\sigma \frac{\partial A_z}{\partial t}=J_{s,z}$$

3.2 温度场方程

3.2.1 热传导方程

温度场满足能量守恒方程：

$$\rho c_p\frac{\partial T}{\partial t}=\nabla \cdot (k\nabla T)+Q_{gen}$$

其中，$Q_{gen}$为体积热源，包括电磁损耗：

$$Q_{gen}=p_{eddy}+p_{hyst}+p_{res}$$

3.2.2 边界条件

温度场边界条件包括：

指定温度：$T=T_0$

对流换热：$-k\nabla T\cdot \mathbf{n}=h(T-T_{\infty })$

辐射换热：$-k\nabla T\cdot \mathbf{n}=\epsilon {\sigma }_{SB}(T^4-T_{\infty }^4)$

热流边界：$-k\nabla T\cdot \mathbf{n}=q_n$

3.3 耦合方程组

综合电磁场和温度场方程，得到电-磁-热耦合控制方程组：

$$\{\begin{array}{ll}\nabla \times \left(\frac{1}{\mu (T)}\nabla \times \mathbf{A}\right)+\sigma (T)\frac{\partial \mathbf{A}}{\partial t}={\mathbf{J}}_s& \text{(电磁场方程)}\\ \rho c_p\frac{\partial T}{\partial t}=\nabla \cdot (k\nabla T)+\frac{|\mathbf{J}{|}^2}{\sigma (T)}+p_{hyst}(T)& \text{(温度场方程)}\end{array}$$

这是一个强耦合非线性方程组，需要迭代求解。

---

4. 数值求解方法

4.1 有限元离散

4.1.1 电磁场弱形式

磁矢势方程的弱形式为：

$${\int }_{\Omega }\frac{1}{\mu }(\nabla \times \mathbf{W})\cdot (\nabla \times \mathbf{A})d\Omega +{\int }_{\Omega }\sigma \mathbf{W}\cdot \frac{\partial \mathbf{A}}{\partial t}d\Omega ={\int }_{\Omega }\mathbf{W}\cdot {\mathbf{J}}_{s}d\Omega$$

其中，$\mathbf{W}$为测试函数。

4.1.2 温度场弱形式

温度场方程的弱形式为：

$${\int }_{\Omega }w\rho c_p\frac{\partial T}{\partial t}d\Omega +{\int }_{\Omega }k\nabla w\cdot \nabla Td\Omega ={\int }_{\Omega }w{Q}_{gen}d\Omega +{\int }_{\Gamma }w{q}_{n}d\Gamma$$

4.2 时谐分析

对于时谐问题，采用相量表示后，电磁场方程变为复数形式：

$${\mathbf{K}}_A\tilde{\mathbf{A}}+j\omega {\mathbf{M}}_{\sigma }\tilde{\mathbf{A}}={\mathbf{f}}_J$$

其中，${\mathbf{K}}_A$为刚度矩阵，${\mathbf{M}}_{\sigma }$为电导率质量矩阵。

涡流损耗计算：

$$P_{eddy}=\frac{1}{2}{\int }_{\Omega }\sigma |\tilde{\mathbf{E}}{|}^{2}d\Omega =\frac{{\omega }^2}{2}{\int }_{\Omega }\sigma |\tilde{\mathbf{A}}{|}^{2}d\Omega$$

4.3 耦合求解策略

4.3.1 顺序耦合

顺序耦合策略交替求解电磁场和温度场：

1. 假设初始温度分布$T^{(0)}$
2. 求解电磁场方程得到${\mathbf{A}}^{(n)}$
3. 计算电磁损耗$Q_{gen}^{(n)}$
4. 求解温度场方程得到$T^{(n)}$
5. 更新材料特性$\mu (T^{(n)})$、$\sigma (T^{(n)})$
6. 检查收敛性，若不收敛返回步骤2

4.3.2 时间步进

对于瞬态问题，采用时间离散：

电磁场：隐式欧拉或Crank-Nicolson

温度场：由于热时间尺度远大于电磁时间尺度，可以采用不同时间步长（多时间尺度方法）

---

5. Python仿真实现

5.1 案例1：一维感应加热分析

问题描述：分析半无限大导体在表面交变磁场作用下的感应加热。

控制方程：一维扩散方程

$$\frac{d^2{\tilde{H}}_y}{d{x}^2}-j\omega \mu \sigma {\tilde{H}}_y=0$$

# -*- coding: utf-8 -*-
"""
主题025：电-磁-热耦合分析 - Python仿真代码
包含6个案例：
1. 一维感应加热分析
2. 圆柱工件感应加热
3. 变压器铁芯损耗分析
4. 电磁热治疗分析
5. 涡流无损检测热成像
6. 磁悬浮系统热分析
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import jv, kv
from scipy.integrate import simpson
import warnings
warnings.filterwarnings('ignore')
import os

output_dir = r'd:\文档\500仿真领域\工程仿真\多物理场耦合仿真\主题025'
os.makedirs(output_dir, exist_ok=True)

plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("="*60)
print("主题025：电-磁-热耦合分析")
print("="*60)

# ============================================================
# 案例1：一维感应加热分析
# ============================================================
print("\n" + "="*60)
print("案例1：一维感应加热分析")
print("="*60)

# 材料参数（钢）
mu = 1000 * 4 * np.pi * 1e-7  # 磁导率 (H/m)
sigma = 1e6  # 电导率 (S/m)
rho = 7850  # 密度 (kg/m³)
cp = 500  # 比热容 (J/kg·K)
k = 50  # 热导率 (W/m·K)

# 激励参数
f = 1000  # 频率 (Hz)
H_surface = 1000  # 表面磁场强度 (A/m)
omega = 2 * np.pi * f

# 趋肤深度
delta = np.sqrt(2 / (omega * mu * sigma))
print(f"趋肤深度: {delta*1000:.3f} mm")

# 解析解：磁场分布
x = np.linspace(0, 10*delta, 500)  # 深度方向
H_y = H_surface * np.exp(-x/delta) * np.cos(x/delta)

# 电流密度
J_z = H_surface / delta * np.exp(-x/delta) * np.cos(x/delta + np.pi/4)

# 功率密度（时间平均）
p_eddy = (H_surface**2 / (2 * sigma * delta**2)) * np.exp(-2*x/delta)

# 温度分布（稳态近似）
# 假设表面温度固定，求解热传导方程
T_surface = 300  # K
T_ambient = 300  # K
h_conv = 10  # 对流换热系数

# 简化：假设所有热量通过对流散失
# 单位面积总功率
P_total = simpson(p_eddy, x)
print(f"单位面积总功率: {P_total:.2f} W/m²")

# 温度升高估计
dT = P_total / h_conv
print(f"估计温升: {dT:.2f} K")

# 可视化
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# 磁场分布
ax1 = axes[0, 0]
ax1.plot(x/delta, H_y/H_surface, 'b-', linewidth=2)
ax1.axhline(y=0, color='k', linewidth=0.5)
ax1.set_xlabel('x/δ', fontsize=11)
ax1.set_ylabel('Hy/H0', fontsize=11)
ax1.set_title('Magnetic Field Distribution', fontsize=12, fontweight='bold')
ax1.grid(True, alpha=0.3)

# 电流密度分布
ax2 = axes[0, 1]
ax2.plot(x/delta, np.abs(J_z)/(H_surface/delta), 'r-', linewidth=2)
ax2.set_xlabel('x/δ', fontsize=11)
ax2.set_ylabel('|Jz|/J0', fontsize=11)
ax2.set_title('Current Density Distribution', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)

# 功率密度分布
ax3 = axes[1, 0]
ax3.semilogy(x/delta, p_eddy/np.max(p_eddy), 'g-', linewidth=2)
ax3.set_xlabel('x/δ', fontsize=11)
ax3.set_ylabel('Power Density (normalized)', fontsize=11)
ax3.set_title('Power Density Distribution', fontsize=12, fontweight='bold')
ax3.grid(True, alpha=0.3)

# 趋肤深度与频率关系
ax4 = axes[1, 1]
freq_range = np.logspace(1, 5, 100)
delta_range = np.sqrt(2 / (2*np.pi*freq_range * mu * sigma))
ax4.loglog(freq_range, delta_range*1000, 'm-', linewidth=2)
ax4.set_xlabel('Frequency (Hz)', fontsize=11)
ax4.set_ylabel('Skin Depth (mm)', fontsize=11)
ax4.set_title('Skin Depth vs Frequency', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/case1_induction_heating_1d.png', dpi=150, bbox_inches='tight')
plt.close()
print("\n✓ 案例1完成，图片已保存")

5.2 案例2：圆柱工件感应加热

问题描述：分析圆柱形钢件在感应线圈中的加热过程。

print("\n" + "="*60)
print("案例2：圆柱工件感应加热")
print("="*60)

# 几何参数
R_rod = 0.01  # 工件半径 (m)
R_coil = 0.015  # 线圈内半径 (m)
N_turns = 10  # 线圈匝数
I_coil = 100  # 线圈电流 (A)

# 材料参数
mu_steel = 1000 * 4 * np.pi * 1e-7
sigma_steel = 1e6

# 频率
f = 10000  # Hz
omega = 2 * np.pi * f

# 趋肤深度
delta = np.sqrt(2 / (omega * mu_steel * sigma_steel))
print(f"趋肤深度: {delta*1000:.3f} mm")
print(f"工件半径/趋肤深度: {R_rod/delta:.2f}")

# 径向网格
nr = 100
r = np.linspace(0, R_rod, nr)
dr = r[1] - r[0]

# 贝塞尔函数解
# H_z(r) = H_0 * J0(k*r) / J0(k*R)
# 其中 k = (1-j)/delta
k = (1 - 1j) / delta

from scipy.special import jv
H_z = np.zeros(nr, dtype=complex)
H_0 = N_turns * I_coil  # 表面磁场强度

for i in range(nr):
    H_z[i] = H_0 * jv(0, k * r[i]) / jv(0, k * R_rod)

# 电流密度
J_phi = np.zeros(nr, dtype=complex)
J_phi[1:] = -np.diff(H_z) / dr

# 功率密度
p_density = 0.5 * sigma_steel * np.abs(J_phi)**2

# 总功率
P_total = 2 * np.pi * np.sum(p_density * r) * dr
print(f"单位长度总功率: {P_total:.2f} W/m")

# 简化温度分析（假设绝热）
rho_steel = 7850
cp_steel = 500
time = 10  # s
dT = P_total * time / (np.pi * R_rod**2 * rho_steel * cp_steel)
print(f"10秒后平均温升: {dT:.2f} K")

# 可视化
fig, axes = plt.subplots(1, 3, figsize=(15, 4))

ax1 = axes[0]
ax1.plot(r*1000, np.abs(H_z)/np.abs(H_0), 'b-', linewidth=2)
ax1.set_xlabel('Radius (mm)', fontsize=11)
ax1.set_ylabel('|Hz|/H0', fontsize=11)
ax1.set_title('Magnetic Field', fontsize=12, fontweight='bold')
ax1.grid(True, alpha=0.3)

ax2 = axes[1]
ax2.plot(r*1000, np.abs(J_phi)/1e6, 'r-', linewidth=2)
ax2.set_xlabel('Radius (mm)', fontsize=11)
ax2.set_ylabel('|Jφ| (MA/m²)', fontsize=11)
ax2.set_title('Current Density', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)

ax3 = axes[2]
ax3.plot(r*1000, p_density/1e6, 'g-', linewidth=2)
ax3.set_xlabel('Radius (mm)', fontsize=11)
ax3.set_ylabel('Power Density (MW/m³)', fontsize=11)
ax3.set_title('Power Density', fontsize=12, fontweight='bold')
ax3.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/case2_cylinder_induction.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例2完成，图片已保存")

5.3 案例3：变压器铁芯损耗分析

问题描述：分析变压器铁芯在不同磁通密度和频率下的损耗特性。

print("\n" + "="*60)
print("案例3：变压器铁芯损耗分析")
print("="*60)

# Steinmetz方程参数（取向硅钢）
k_h = 0.02  # 磁滞损耗系数
n = 2.0  # Steinmetz指数
k_e = 1e-4  # 涡流损耗系数
k_r = 5e-5  # 剩余损耗系数

# 参数范围
B_range = np.linspace(0.1, 1.8, 50)  # 磁通密度 (T)
f_range = np.array([50, 100, 200, 400, 800])  # 频率 (Hz)

# 计算损耗
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

ax1 = axes[0]
colors = plt.cm.viridis(np.linspace(0, 1, len(f_range)))

for i, f in enumerate(f_range):
    p_hyst = k_h * f * B_range**n
    p_eddy = k_e * (f * B_range)**2
    p_res = k_r * (f * B_range)**1.5
    p_total = p_hyst + p_eddy + p_res
    
    ax1.plot(B_range, p_total, color=colors[i], linewidth=2, label=f'f={f}Hz')

ax1.set_xlabel('Flux Density B (T)', fontsize=11)
ax1.set_ylabel('Core Loss (W/kg)', fontsize=11)
ax1.set_title('Core Loss vs Flux Density', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 损耗分量分析
ax2 = axes[1]
f_fixed = 50  # Hz
B_fixed = 1.5  # T

p_hyst = k_h * f_fixed * B_fixed**n
p_eddy = k_e * (f_fixed * B_fixed)**2
p_res = k_r * (f_fixed * B_fixed)**1.5

losses = [p_hyst, p_eddy, p_res]
labels = ['Hysteresis', 'Eddy Current', 'Residual']
colors_pie = ['#ff9999', '#66b3ff', '#99ff99']

ax2.pie(losses, labels=labels, colors=colors_pie, autopct='%1.1f%%', startangle=90)
ax2.set_title(f'Loss Components (f={f_fixed}Hz, B={B_fixed}T)', fontsize=12, fontweight='bold')

plt.tight_layout()
plt.savefig(f'{output_dir}/case3_transformer_loss.png', dpi=150, bbox_inches='tight')
plt.close()

print(f"\n变压器损耗分析 (f={f_fixed}Hz, B={B_fixed}T):")
print(f"  磁滞损耗: {p_hyst:.4f} W/kg")
print(f"  涡流损耗: {p_eddy:.4f} W/kg")
print(f"  剩余损耗: {p_res:.4f} W/kg")
print(f"  总损耗: {p_hyst+p_eddy+p_res:.4f} W/kg")
print("✓ 案例3完成，图片已保存")

5.4 案例4：电磁热治疗分析

问题描述：分析磁纳米颗粒在交变磁场中的产热及组织温度分布。

print("\n" + "="*60)
print("案例4：电磁热治疗分析")
print("="*60)

# 组织参数
k_tissue = 0.5  # 热导率 (W/m·K)
rho_tissue = 1000  # 密度 (kg/m³)
cp_tissue = 3600  # 比热容 (J/kg·K)
perfusion = 0.001  # 血液灌注率 (1/s)
T_blood = 37  # 血液温度 (°C)

# 纳米颗粒参数
M_s = 400e3  # 饱和磁化强度 (A/m)
H_ac = 10000  # 交变磁场强度 (A/m)
f_magnetic = 100e3  # 磁场频率 (Hz)
concentration = 0.01  # 体积浓度

# 比吸收率 (SAR) 计算
# SAR = (π * μ0 * f * H * M_s * concentration) / ρ
mu0 = 4 * np.pi * 1e-7
SAR = (np.pi * mu0 * f_magnetic * H_ac * M_s * concentration) / rho_tissue
print(f"比吸收率 SAR: {SAR:.2f} W/kg")

# 球形肿瘤模型
R_tumor = 0.005  # 肿瘤半径 (m)
nr = 100
r = np.linspace(0, 3*R_tumor, nr)

# 稳态温度分布（Pennes生物热方程简化）
# d²T/dr² + (2/r)dT/dr + (SAR - perfusion*(T-T_blood))/k = 0

# 解析解（均匀热源球体）
T_center = T_blood + SAR / (3 * k_tissue) * R_tumor**2
T_surface = T_blood + SAR / (3 * k_tissue) * R_tumor**2 / (1 + k_tissue/(h_conv*R_tumor))

# 数值求解
T_steady = np.zeros(nr)
h_conv = 10  # 表面换热系数

for i in range(nr):
    if r[i] <= R_tumor:
        # 肿瘤内部
        T_steady[i] = T_blood + SAR / (6 * k_tissue) * (3*R_tumor**2 - r[i]**2)
    else:
        # 肿瘤外部
        T_steady[i] = T_blood + (SAR * R_tumor**3) / (3 * k_tissue * r[i])

# 考虑血液灌注的修正
T_steady = T_steady - (SAR * perfusion * rho_tissue * cp_tissue / k_tissue) * (r**2)

print(f"\n热治疗温度分析:")
print(f"  肿瘤中心温度: {T_center:.2f} °C")
print(f"  肿瘤表面温度: {T_surface:.2f} °C")
print(f"  治疗温度范围: 42-45 °C (热疗)")
print(f"  消融温度: >50 °C")

# 可视化
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

ax1 = axes[0]
ax1.plot(r*1000, T_steady, 'r-', linewidth=2)
ax1.axvline(x=R_tumor*1000, color='k', linestyle='--', label='Tumor boundary')
ax1.axhline(y=42, color='orange', linestyle='--', alpha=0.7, label='Hyperthermia')
ax1.axhline(y=50, color='red', linestyle='--', alpha=0.7, label='Ablation')
ax1.set_xlabel('Radius (mm)', fontsize=11)
ax1.set_ylabel('Temperature (°C)', fontsize=11)
ax1.set_title('Temperature Distribution', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 参数敏感性分析
ax2 = axes[1]
H_range = np.linspace(5000, 20000, 50)
SAR_range = (np.pi * mu0 * f_magnetic * H_range * M_s * concentration) / rho_tissue
T_max_range = T_blood + SAR_range / (3 * k_tissue) * R_tumor**2

ax2.plot(H_range/1000, T_max_range, 'b-', linewidth=2)
ax2.axhline(y=42, color='orange', linestyle='--', alpha=0.7, label='Hyperthermia threshold')
ax2.axhline(y=50, color='red', linestyle='--', alpha=0.7, label='Ablation threshold')
ax2.set_xlabel('Field Strength H (kA/m)', fontsize=11)
ax2.set_ylabel('Max Temperature (°C)', fontsize=11)
ax2.set_title('Temperature vs Field Strength', fontsize=12, fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/case4_hyperthermia.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例4完成，图片已保存")

5.5 案例5：涡流无损检测热成像

问题描述：利用涡流加热和热成像检测材料缺陷。

print("\n" + "="*60)
print("案例5：涡流无损检测热成像")
print("="*60)

# 平板参数
L_plate = 0.1  # 板长 (m)
W_plate = 0.1  # 板宽 (m)
t_plate = 0.005  # 板厚 (m)

# 材料参数
sigma_al = 3.5e7  # 铝电导率
mu_al = 4 * np.pi * 1e-7
k_al = 237

# 缺陷参数
defect_x, defect_y = 0.05, 0.05  # 缺陷位置
defect_size = 0.01  # 缺陷尺寸
defect_depth = 0.002  # 缺陷深度

# 涡流激励
f_eddy = 1000  # Hz
I_coil = 5  # A
N_coil = 50

# 简化模型：计算表面温度分布
nx, ny = 50, 50
x = np.linspace(0, L_plate, nx)
y = np.linspace(0, W_plate, ny)
X, Y = np.meshgrid(x, y)

# 无缺陷区域功率密度
omega = 2 * np.pi * f_eddy
delta_al = np.sqrt(2 / (omega * mu_al * sigma_al))
H_surface = N_coil * I_coil
p_normal = 0.5 * sigma_al * (omega * mu_al * H_surface * delta_al)**2

# 缺陷区域（电导率降低，功率密度变化）
defect_mask = ((X - defect_x)**2 + (Y - defect_y)**2 < (defect_size/2)**2)
sigma_defect = sigma_al * 0.1  # 缺陷处电导率降低
p_defect = 0.5 * sigma_defect * (omega * mu_al * H_surface * delta_al)**2

# 功率密度分布
P_density = np.ones((ny, nx)) * p_normal
P_density[defect_mask] = p_defect

# 稳态温度分布（简化）
T_surface = 300 + P_density / 100  # 简化的温升计算

print(f"\n涡流检测分析:")
print(f"  铝趋肤深度: {delta_al*1000:.3f} mm")
print(f"  正常区域功率密度: {p_normal/1e3:.2f} kW/m³")
print(f"  缺陷区域功率密度: {p_defect/1e3:.2f} kW/m³")
print(f"  温度对比度: {np.max(T_surface) - np.min(T_surface):.2f} K")

# 可视化
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

ax1 = axes[0]
im1 = ax1.contourf(X*1000, Y*1000, P_density/1e6, levels=20, cmap='hot')
ax1.set_xlabel('x (mm)', fontsize=11)
ax1.set_ylabel('y (mm)', fontsize=11)
ax1.set_title('Eddy Current Power Density', fontsize=12, fontweight='bold')
ax1.set_aspect('equal')
cbar1 = plt.colorbar(im1, ax=ax1)
cbar1.set_label('Power Density (MW/m³)', fontsize=10)

ax2 = axes[1]
im2 = ax2.contourf(X*1000, Y*1000, T_surface, levels=20, cmap='jet')
ax2.set_xlabel('x (mm)', fontsize=11)
ax2.set_ylabel('y (mm)', fontsize=11)
ax2.set_title('Surface Temperature', fontsize=12, fontweight='bold')
ax2.set_aspect('equal')
cbar2 = plt.colorbar(im2, ax=ax2)
cbar2.set_label('Temperature (K)', fontsize=10)

plt.tight_layout()
plt.savefig(f'{output_dir}/case5_thermography.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例5完成，图片已保存")

5.6 案例6：磁悬浮系统热分析

问题描述：分析磁悬浮轴承中的电磁损耗和温升。

print("\n" + "="*60)
print("案例6：磁悬浮系统热分析")
print("="*60)

# 磁悬浮轴承参数
B_bias = 0.5  # 偏置磁通密度 (T)
B_control = 0.2  # 控制磁通密度 (T)
f_switch = 10000  # 开关频率 (Hz)

# 材料参数
mu_core = 5000 * 4 * np.pi * 1e-7  # 铁芯磁导率
sigma_core = 1e6  # 铁芯电导率
rho_core = 7800
cp_core = 500

# 几何参数
V_core = 0.001  # 铁芯体积 (m³)
A_surface = 0.05  # 散热表面积 (m²)

# 损耗计算
# 磁滞损耗
k_h = 0.05
p_hyst = k_h * f_switch * (B_bias + B_control)**2

# 涡流损耗
k_e = 1e-4
p_eddy = k_e * (f_switch * B_control)**2

# 铜损（简化）
I_coil = 2  # A
R_coil = 5  # Ohm
p_copper = I_coil**2 * R_coil / V_core

p_total = p_hyst + p_eddy + p_copper

print(f"\n磁悬浮轴承损耗分析:")
print(f"  磁滞损耗: {p_hyst:.2f} W/m³")
print(f"  涡流损耗: {p_eddy:.2f} W/m³")
print(f"  铜损: {p_copper:.2f} W/m³")
print(f"  总损耗密度: {p_total:.2f} W/m³")
print(f"  总损耗功率: {p_total * V_core:.2f} W")

# 稳态温度
h_ambient = 10  # 环境换热系数
T_ambient = 300  # K
T_steady = T_ambient + p_total * V_core / (h_ambient * A_surface)
print(f"  稳态温度: {T_steady:.2f} K ({T_steady-273.15:.2f} °C)")

# 瞬态温度分析
time = np.linspace(0, 3600, 1000)  # 1小时
tau = rho_core * cp_core * V_core / (h_ambient * A_surface)  # 时间常数
T_transient = T_ambient + (T_steady - T_ambient) * (1 - np.exp(-time/tau))

# 不同控制策略比较
strategies = {
    'Low Bias': {'B_bias': 0.3, 'B_control': 0.2},
    'Medium Bias': {'B_bias': 0.5, 'B_control': 0.2},
    'High Bias': {'B_bias': 0.7, 'B_control': 0.2},
    'PWM Control': {'B_bias': 0.5, 'B_control': 0.3}
}

results = {}
for name, params in strategies.items():
    p_h = k_h * f_switch * (params['B_bias'] + params['B_control'])**2
    p_e = k_e * (f_switch * params['B_control'])**2
    p_c = I_coil**2 * R_coil / V_core
    p_t = p_h + p_e + p_c
    T_s = T_ambient + p_t * V_core / (h_ambient * A_surface)
    results[name] = {'power': p_t, 'temp': T_s}

# 可视化
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# 瞬态温度
ax1 = axes[0, 0]
ax1.plot(time/60, T_transient-273.15, 'b-', linewidth=2)
ax1.axhline(y=T_steady-273.15, color='r', linestyle='--', label='Steady state')
ax1.set_xlabel('Time (min)', fontsize=11)
ax1.set_ylabel('Temperature (°C)', fontsize=11)
ax1.set_title('Transient Temperature Rise', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 损耗分量
ax2 = axes[0, 1]
loss_types = ['Hysteresis', 'Eddy Current', 'Copper']
loss_values = [p_hyst*V_core, p_eddy*V_core, p_copper*V_core]
colors = ['#ff9999', '#66b3ff', '#99ff99']
ax2.bar(loss_types, loss_values, color=colors, alpha=0.7)
ax2.set_ylabel('Power (W)', fontsize=11)
ax2.set_title('Loss Breakdown', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3, axis='y')

# 不同策略比较
ax3 = axes[1, 0]
names = list(results.keys())
temps = [results[n]['temp']-273.15 for n in names]
powers = [results[n]['power']*V_core for n in names]

x_pos = np.arange(len(names))
bars = ax3.bar(x_pos, temps, alpha=0.7)
ax3.set_xticks(x_pos)
ax3.set_xticklabels(names, rotation=15, ha='right')
ax3.set_ylabel('Temperature (°C)', fontsize=11)
ax3.set_title('Temperature for Different Strategies', fontsize=12, fontweight='bold')
ax3.grid(True, alpha=0.3, axis='y')

ax4 = axes[1, 1]
ax4.scatter(powers, temps, s=100, alpha=0.7)
for i, name in enumerate(names):
    ax4.annotate(name, (powers[i], temps[i]), fontsize=8, ha='center')
ax4.set_xlabel('Total Power (W)', fontsize=11)
ax4.set_ylabel('Temperature (°C)', fontsize=11)
ax4.set_title('Power vs Temperature', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/case6_maglev_thermal.png', dpi=150, bbox_inches='tight')
plt.close()
print("\n✓ 案例6完成，图片已保存")

print("\n" + "="*60)
print("所有案例仿真完成！")
print("="*60)
print("\n生成的图片文件:")
print("  - case1_induction_heating_1d.png")
print("  - case2_cylinder_induction.png")
print("  - case3_transformer_loss.png")
print("  - case4_hyperthermia.png")
print("  - case5_thermography.png")
print("  - case6_maglev_thermal.png")
print("="*60)