多物理场耦合算法

电磁热耦合：工程设置与诊断验证

从焦耳热密度与趋肤深度估算出发，给出频域电磁与瞬态热串联的网格要求、损耗分离参数、ANSYS 物理环境配置和电磁热参数取值表。

林工程师16 分钟更新于 2026/9/22
CAE结果诊断与可信度验证工程设置与参数选择多物理场耦合算法电磁热耦合焦耳热趋肤深度工程

电磁热耦合：工程设置与诊断验证

电磁热耦合的设置顺序与多数多物理场相反：应当先算电磁、再定热网格，因为导体内的热源密度完全由电流密度与频率决定，而频率又通过趋肤深度决定电磁网格的最小尺寸。电磁热耦合的典型流程是电磁求解器算出损耗密度，热求解器据此算温度场，再把温度回代修正电导率。这类结果最容易被两类伪影污染：把趋肤效应算成体积均匀损耗，以及把温度反馈迭代的未收敛当作稳态。

趋肤深度决定网格

交变电流集中在导体表层，特征深度为

δ=2ωμσ=1πfμ0μrσ\delta=\sqrt{\frac{2}{\omega\mu\sigma}}=\sqrt{\frac{1}{\pi f\mu_0\mu_r\sigma}}

对铜（μr=1\mu_r=1、σ=5.96×107 S/m\sigma=5.96\times10^7\ \mathrm{S/m}），50 Hz 时 δ=1/(π×50×4π×10−7×5.96×107)≈9.2×10−3 m=9.2 mm\delta=\sqrt{1/(\pi\times50\times4\pi\times10^{-7}\times5.96\times10^7)}\approx9.2\times10^{-3}\ \mathrm{m}=9.2\ \mathrm{mm}；1 MHz 时 δ\delta 按 50/106\sqrt{50/10^6} 缩小，约 6.5×10−5 m=0.065 mm6.5\times10^{-5}\ \mathrm{m}=0.065\ \mathrm{mm}。

网格必须能在表层内放下至少 3 个单元，即导体表层单元尺寸 he≤δ/3h_e\le\delta/3：50 Hz 时为 3.1 mm，1 MHz 时骤降到 0.022 mm。这个反差解释了为什么高频电磁热问题必须用表层扫掠网格或阻抗边界条件，而不能用体网格均匀剖分。

趋肤深度与网格/频率的匹配

交流下电流集中在表面，特征深度为

δ=2ωμσ\delta=\sqrt{\frac{2}{\omega\mu\sigma}}

铜在 50 Hz 时：ω=2π×50=314.2 rad/s\omega=2\pi\times50=314.2\ \text{rad/s}、μ=μ0=1.2566×10−6 H/m\mu=\mu_0=1.2566\times10^{-6}\ \text{H/m}，ωμσ=314.2×1.2566×10−6×5.96×107=2.35×104\omega\mu\sigma=314.2\times1.2566\times10^{-6}\times5.96\times10^{7}=2.35\times10^{4}，故 δ=2/2.35×104=9.2×10−3 m\delta=\sqrt{2/2.35\times10^{4}}=9.2\times10^{-3}\ \text{m}，即 9.2 mm。导体半径若只有 3 mm，则 r/δ=0.33<1r/\delta=0.33<1，电流近似均匀，可以用体损耗；半径 15 mm 时 r/δ=1.63r/\delta=1.63，必须分层剖分。

判定规则：导体表层 1.5δ 范围内至少布置 3 层单元。9.2 mm 对应首层厚度约 3.1 mm、第二层 3.1 mm、第三层 3.0 mm。若网格首层厚 10 mm 而 δ 只有 9.2 mm，损耗会被系统性高估——这是"电流密度峰值算大了"的最常见来源。频率升高到 1 kHz 时 δ∝1/f\delta\propto1/\sqrt{f}，降到 9.2/20=2.06 mm9.2/\sqrt{20}=2.06\ \text{mm}，同一网格立刻失效，因此频率变化后必须重查网格。

关键参数取值表

参数 典型取值 依据
铜电导率 σ\sigma 5.96×107 S/m5.96\times10^7\ \mathrm{S/m} 20 ℃ 标准值
电阻温度系数 0.0039 K−10.0039\ \mathrm{K^{-1}} 铜，需反馈到 σ(T)\sigma(T)
表层单元尺寸 heh_e ≤δ/3\le\delta/3 分辨趋肤深度
电磁时间步 谐波分析无需步长 周期平均热源
热时间步 ≤τth/20\le\tau_{th}/20 分辨热瞬态
磁滞系数 khk_h 150 W/(m3T1.8Hz)150\ \mathrm{W/(m^3T^{1.8}Hz)} 电工钢实测拟合
涡流系数 kek_e 0.5 W/(m3T2Hz2)0.5\ \mathrm{W/(m^3T^2Hz^2)} 叠片厚度相关

ANSYS 物理环境配置

多物理场串联用物理环境（physics environment）管理，先把电磁设置写成文件，切换单元类型后再写热设置，求解时交替读入：

/PREP7
ET,1,SOLID236,0          ! 边棱元, 频域电磁
MP,MURX,1,1
MP,RSVX,1,1.68E-8        ! 铜电阻率, Ohm*m
/SOLU
ANTYPE,HARMIC
HARFRQ,50                ! 50 Hz
PHYSICS,WRITE,emag       ! 写出电磁物理环境
FINISH

/PREP7
ET,1,SOLID278             ! 切换到热单元
MP,KXX,1,400
PHYSICS,WRITE,thermal
FINISH

/SOLU
PHYSICS,READ,emag
SOLVE
PHYSICS,READ,thermal
LDREAD,HGEN,,,,2,'emag','rst'   ! 读入焦耳热作为体热源
ANTYPE,TRANS
SOLVE

LDREAD,HGEN 把电磁结果中的热生成率直接读成热分析的体载荷，是串联耦合的关键一行；若忘记这一步，热分析会得到零热源却依然收敛，属于最难察觉的设置错误。

焦耳热与损耗密度的量级

导电介质中的焦耳热密度为

q=∣J∣2σq=\frac{\lvert\mathbf{J}\rvert^2}{\sigma}

铜的电导率 σ=5.96×107 S/m\sigma=5.96\times10^7\ \mathrm{S/m}。取电流密度 J=5×106 A/m2J=5\times10^6\ \mathrm{A/m^2}，则 q=(5×106)2/5.96×107=4.19×105 W/m3q=(5\times10^6)^2/5.96\times10^7=4.19\times10^5\ \mathrm{W/m^3}，即 419 kW/m³。这个量级必须先算出来，才能判断后续要不要考虑温度对电导率的反馈——铜的电阻温度系数约 0.0039 K−10.0039\ \mathrm{K^{-1}}，温升 100 K 会让 σ\sigma 下降近 30%，进而改变热源本身。

铁磁材料的损耗还要加上磁滞与涡流两部分，用损耗分离写成

Pv=khfBp1.8+ke(fBp)2P_v=k_h f B_p^{1.8}+k_e (f B_p)^2

第一项是磁滞损耗，第二项是经典涡流损耗。取 kh=150 W/(m3⋅T1.8⋅Hz)k_h=150\ \mathrm{W/(m^3\cdot T^{1.8}\cdot Hz)}、ke=0.5 W/(m3⋅T2⋅Hz2)k_e=0.5\ \mathrm{W/(m^3\cdot T^2\cdot Hz^2)}、f=50 Hzf=50\ \mathrm{Hz}、Bp=1.5 TB_p=1.5\ \mathrm{T}，得 Pv=150×50×1.51.8+0.5×(50×1.5)2≈150×50×2.05+0.5×5625≈1.54×104+2.81×103≈1.82×104 W/m3P_v=150\times50\times1.5^{1.8}+0.5\times(50\times1.5)^2\approx150\times50\times2.05+0.5\times5625\approx1.54\times10^4+2.81\times10^3\approx1.82\times10^4\ \mathrm{W/m^3}。按密度 7650 kg/m³ 折算约 2.4 W/kg，与常用电工钢在 1.5 T、50 Hz 下的损耗量级一致。

频域电磁到瞬态热的串联

工程上常做两步串联：先用频域（时谐）电磁求出周期平均热源，再把热源映射到热网格上做瞬态传热。周期平均焦耳热为

qˉ=12∣J^∣2σ\bar{q}=\frac{1}{2}\frac{\lvert\hat{\mathbf{J}}\rvert^2}{\sigma}

其中 J^\hat{\mathbf{J}} 是复相量幅值。之所以能用周期平均值，是因为电磁时间尺度（微秒）与热时间尺度（秒）相差六个数量级以上，温度在一个电磁周期内几乎不变。判据是两个尺度之比小于 10−310^{-3}，否则必须做瞬态电磁与瞬态热的双向耦合。

热源映射后，表面散热主导的温升可用一维估算校核：

ΔTmax=qL22k+qLh\Delta T_{max}=\frac{qL^2}{2k}+\frac{qL}{h}

取 q=4.19×105 W/m3q=4.19\times10^5\ \mathrm{W/m^3}、导体半厚 L=5 mmL=5\ \mathrm{mm}、铜 k=400 W/(m⋅K)k=400\ \mathrm{W/(m\cdot K)}、自然对流 h=10 W/(m2⋅K)h=10\ \mathrm{W/(m^2\cdot K)}，得第一项仅 0.013 K，第二项 4.19×105×0.005/10=210 K4.19\times10^5\times0.005/10=210\ \mathrm{K}。可见温升几乎全部由表面散热能力决定，内部导热根本不是瓶颈——这直接告诉设计者该加散热面积而不是加铜。

故障模式与判定试验

现象 根因 判定试验
高频下损耗远低于手册值 网格未分辨趋肤深度 表层加密到 δ/3\delta/3，看损耗是否上升
温升随迭代发散 电导率温度反馈未加限幅 给 σ(T)\sigma(T) 设上下界，观察是否稳定
热分析收敛但温度不升 未执行 LDREAD,HGEN 检查热源体载荷是否为零
叠片涡流损耗偏高 用实体代替叠片各向异性 改用分层或等效电阻率重算
频率翻倍损耗不按平方增长 未计入磁滞项 分离两项，单独扫频验证斜率
损耗比手算高 2 倍以上 电流密度用了峰值而非有效值 按 q=J2/σq=J^2/\sigma 手算并与输出积分对比
高频下损耗几乎不随频率变 网格未分辨趋肤层 算 δ\delta，检查表层 1.5δ 内层数是否 ≥3
温升比集总估算低一半 散热边界用了过大 hh 或漏掉辐射 用 ΔT=P/(hA)\Delta T=P/(hA) 反算等效 hh
温度反馈迭代 20 次仍漂移 恒流/恒压激励混淆 确认激励类型，核对 q∝ρq\propto\rho 的方向
三维温度远低于一维估算 热源体积积分漏项 核对 q⋅Vq\cdot V 与 PinP_{in} 账本
瞬态算 5 分钟就报稳态 时间常数被低估 算 τ=ρcpV/(hA)\tau=\rho c_pV/(hA)，至少积分 3τ
温度场出现网格状波纹 电磁与热网格不一致且未守恒映射 检查损耗密度映射是否做体积加权归一

能量账本：第一诊断量

电磁求解器输出的焦耳损耗密度为

q=J2σ=ρJ2q = \frac{J^{2}}{\sigma} = \rho J^{2}

其中 JJ 为电流密度、σ\sigma 为电导率、ρ=1/σ\rho=1/\sigma 为电阻率。诊断第一步是把 qq 在导体体积上积分，与输入电功率对比：

Ploss=∫Vq dV=?Pin−Pstored−PradP_{loss}=\int_V q\,dV \stackrel{?}{=} P_{in}-P_{stored}-P_{rad}

取铜导体 σ=5.96×107 S/m\sigma=5.96\times10^{7}\ \text{S/m}、J=5 A/mm2=5.0×106 A/m2J=5\ \text{A/mm}^2=5.0\times10^{6}\ \text{A/m}^2，则 q=(5.0×106)2/5.96×107=4.19×105 W/m3q=(5.0\times10^{6})^{2}/5.96\times10^{7}=4.19\times10^{5}\ \text{W/m}^3，即 0.419 MW/m³。若导体体积 V=1.0×10−4 m3V=1.0\times10^{-4}\ \text{m}^3，则 Ploss=4.19×105×1.0×10−4=41.9 WP_{loss}=4.19\times10^{5}\times1.0\times10^{-4}=41.9\ \text{W}。这个数就是后续热分析的唯一热源，任何偏离 41.9 W 超过 2% 的输入都应先查电流密度定义（峰值还是有效值）与单位换算。

与解析稳态解和实测温升对照

用集总热阻估算平衡温度，作为三维热分析的独立对照。自然对流 h=10 W/(m2⋅K)h=10\ \text{W/(m}^2\cdot\text{K)}、散热面积 A=0.02 m2A=0.02\ \text{m}^2、损耗 41.9 W，则

ΔT=PlosshA=41.910×0.02=209.5 K\Delta T=\frac{P_{loss}}{hA}=\frac{41.9}{10\times0.02}=209.5\ \text{K}

该温升对多数绝缘材料不可接受，说明必须改用强迫风冷（h≈50 W/(m2⋅K)h\approx50\ \text{W/(m}^2\cdot\text{K)} 时 ΔT=41.9\Delta T=41.9 K）或增大散热面积。这个集总算例的价值在于：它给出温度量级的上下界，任何三维结果若远低于 41.9 K 而散热条件只有自然对流，几乎必定是边界条件设置错误。

热时间常数 τ=ρcpV/(hA)\tau=\rho c_p V/(hA)：取铜 ρcp=3.45×106 J/(m3⋅K)\rho c_p=3.45\times10^{6}\ \text{J/(m}^3\cdot\text{K)}、V=1.0×10−4 m3V=1.0\times10^{-4}\ \text{m}^3、hA=0.2 W/KhA=0.2\ \text{W/K}，得 τ=3.45×106×1.0×10−4/0.2=1725 s\tau=3.45\times10^{6}\times1.0\times10^{-4}/0.2=1725\ \text{s}，约 29 分钟。这意味着瞬态算例至少要算 3τ≈86 分钟才能接近稳态，只算 5 分钟就宣称"稳态温度"是常见的误判。

# 电磁-热单向/双向耦合的最小验证脚本
sigma0, alpha, T0 = 5.96e7, 3.93e-3, 20.0   # S/m, 1/K, degC
J, V, h, A = 5.0e6, 1.0e-4, 10.0, 0.02       # A/m^2, m^3, W/m^2K, m^2
T = 20.0
for it in range(30):
    sigma = sigma0 / (1.0 + alpha * (T - T0))    # 电导率随温度下降
    q = J * J / sigma                            # 焦耳损耗密度
    P = q * V                                    # 总损耗
    T_new = T0 + P / (h * A)                     # 集总热平衡
    if abs(T_new - T) < 1.0e-3 * max(T, 1.0):
        break
    T = 0.5 * T + 0.5 * T_new                    # 松弛 0.5
print("P_loss = %.2f W, T = %.1f degC, iters = %d" % (P, T, it))

电阻率的温度反馈与迭代收敛

温度改变电导率，电导率又改变损耗，形成反馈：

ρ(T)=ρ0[1+α(T−T0)]\rho(T)=\rho_{0}\left[1+\alpha\left(T-T_{0}\right)\right]

铜在 20 ℃ 时 ρ0=1.68×10−8 Ω⋅m\rho_0=1.68\times10^{-8}\ \Omega\cdot\text{m}、α=3.93×10−3 K−1\alpha=3.93\times10^{-3}\ \text{K}^{-1}。若稳态温度为 100 ℃，则 ρ=1.68×10−8×[1+3.93×10−3×80]=1.68×10−8×1.314=2.21×10−8 Ω⋅m\rho=1.68\times10^{-8}\times[1+3.93\times10^{-3}\times80]=1.68\times10^{-8}\times1.314=2.21\times10^{-8}\ \Omega\cdot\text{m}，电阻率上升 31.4%。在电流受控的前提下损耗 q∝ρq\propto\rho，因此损耗同步上升 31.4%；若电流受控而电压固定，损耗反而下降约 24%。这两种情形的结论相反，诊断时必须先确认激励是恒流还是恒压。

反馈迭代用松弛 Tk+1=(1−ω)Tk+ω H(Tk)T^{k+1}=(1-\omega)T^{k}+\omega\,\mathcal{H}(T^{k})，ω\omega 取 0.3～0.7。收敛判据用相对温度残差 ∣ΔT∣/ΔTref<10−3|\Delta T|/\Delta T_{ref}<10^{-3}，且必须同时报告迭代次数；若迭代 20 次仍在漂移 0.5 K 以上，说明温度系数与损耗模型不自洽，而不是"收敛慢"。

参考资料

  1. Steinmetz C.P., “On the law of hysteresis,” Transactions of the AIEE, 9, 1892.
  2. Bertotti G., “General properties of power losses in soft ferromagnetic materials,” IEEE Transactions on Magnetics, 24(1), 1988.
  3. Ida N., Numerical Modeling for Electromagnetic Non-Destructive Evaluation, Chapman & Hall, 1995.
  4. Bossavit A., Computational Electromagnetism, Academic Press, 1998.
  5. ANSYS Inc., Maxwell 3D User’s Guide, Release 2021R2.
  6. Bastos J.P.A., Sadowski N., Electromagnetic Modeling by Finite Element Methods, Marcel Dekker, 2003.
  7. Ida N., Bastos J.P.A., Electromagnetics and Calculation of Fields, 2nd ed., Springer, 1997.
  8. Biro O., Preis K., “On the use of the magnetic vector potential in the finite element analysis of three-dimensional eddy currents,” IEEE Transactions on Magnetics, 25(4), 1989.
  9. Incropera F.P., DeWitt D.P., Bergman T.L., Lavine A.S., Fundamentals of Heat and Mass Transfer, 7th ed., Wiley, 2011.
  10. Meeker D., Finite Element Method Magnetics: User’s Manual, 2019.
  11. IEC 60287-1-1, Electric cables — Calculation of the current rating — Part 1-1: Current rating equations and calculation of losses, 2014.
  12. Bungartz H.-J., Schäfer M. (eds.), Fluid-Structure Interaction: Modelling, Simulation, Optimisation, Springer, 2006.