超导失超模拟
理论与物理
失超现象的本质
老师，失超就是超导突然失效吗？它有多危险？
问得好。失超是指超导体局部温度超过临界温度 $T_c$ 而转变为常导态的现象。在超导态下电阻为零，电流可以无损耗流动，但只要有一点转变为常导态，该部分就会突然产生电阻。
产生电阻后，就会在那里产生焦耳热，对吧？
没错。而且超导磁体中流过大电流。以LHC的超导双极磁体为例，为了产生8特斯拉的磁场，流过的电流高达约11,850 A。一个线圈的电感约为100 mH，因此储存的磁场能量为
这相当于一辆汽车以160公里时速行驶时的动能。一旦发生失超，这部分能量将在磁体内部转化为热能，如果无法控制，磁体可能会熔化，或者液氦急剧汽化导致管道破裂。
诶，真的发生过这样的事故吗？
2008年9月，LHC的3-4区段发生的事故非常有名。超导连接处的一个焊点不良（电阻值仅约220 nΩ）引发了失超。失超保护系统在某个区间未能按设计工作，电弧放电击穿了氦容器。约6吨液氦瞬间汽化，导致53台双极磁体及周边基础设施受损。最终，修复和安全强化措施花费了14个月。
220纳欧姆的不良导致停机14个月…所以失超传播模拟才这么重要啊。
正是如此。失超模拟主要有三个目的：
- 热点温度预测 — 确认温度上升最高的点是否超过磁体材料的允许温度（通常300~400 K）
- 失超检测时间确定 — 量化保护系统做出反应前的允许延迟时间
- 保护电路参数优化 — 确定泄放电阻值、失超加热器的布置和发热量
控制方程
失超传播在数学上如何描述呢？
基本是在热传导方程中加入焦耳热源项。超导线的温度分布 $T(\mathbf{x}, t)$ 遵循以下方程：
各项含义如下：
- $\rho \, c_p$ — 体积比热容 [J/(m³·K)]。超导线是复合材料（NbTi + Cu + 环氧树脂），因此使用有效值
- $k(T)$ — 热导率 [W/(m·K)]。极低温下温度依赖性极强（铜在4 K时约600 W/(m·K)，30 K时约2000 W/(m·K)，300 K时约400 W/(m·K)）
- $\rho_\text{el}(T, B)$ — 电阻率 [Ω·m]。超导态（$T < T_c$）下为零，常导态下为铜基体的电阻率。在 $T_c$ 附近存在电流共享的过渡区域
- $J$ — 电流密度 [A/m²]
- $Q_\text{ext}$ — 外部热源（如失超加热器）[W/m³]
$\rho_\text{el}$ 在超导态为零，在常导态为有限值…也就是说，这项是失超的“开关”对吧？
很敏锐。准确地说，超导体还有另一个临界温度，称为电流共享温度 $T_\text{cs}$。当 $T < T_\text{cs}$ 时，所有电流流过超导细丝，发热为零。当 $T_\text{cs} < T < T_c$ 时，部分电流分流到铜基体，开始发热。当 $T > T_c$ 时，完全转变为常导态。精确建模这个过渡区域是模拟精度的关键。
其中 $T_\text{op}$ 是运行温度，$I_\text{op}$ 是运行电流，$I_c$ 是临界电流。
常导区传播速度（NZPV）
一旦失超开始，会以多快的速度扩散？
常导区传播速度（NZPV: Normal Zone Propagation Velocity）是描述失超扩散速度最重要的参数。假设一维稳态传播，可以用Wilson的经典公式近似：
其中 $C(T_c)$ 是临界温度下的体积比热容。
具体速度是多少呢？
不同材料差异巨大。这正是超导磁体设计的核心问题。
| 线材 | $T_c$ [K] | 典型NZPV | 特征 |
|---|---|---|---|
| NbTi | 9.2 | 5〜20 m/s | LHC等应用广泛。NZPV快，保护相对容易 |
| Nb₃Sn | 18.3 | 1〜10 m/s | 计划用于FCC・ITER-CS。脆性是挑战 |
| REBCO (HTS) | 约92 | 0.1〜50 mm/s | NZPV极慢 → 局部过热风险严重 |
| BSCCO-2223 | 约110 | 1〜100 mm/s | 带状导体。各向异性强 |
请注意，REBCO（高温超导体）的NZPV只有NbTi的1/100到1/1000。NZPV慢意味着产生的热量被限制在狭窄区域，导致热点温度急剧上升。这就是HTS磁体保护设计极其困难的原因。
热点温度
热点温度具体怎么计算？
估算热点温度上限时使用绝热假设。作为最坏情况，假设失超发生点的热传导为零。则热平衡为：
整理两边并积分，得到：
左边是材料固有的函数，以温度上限 $T_\text{max}$ 为参数。右边是电流的时间积分（$J^2 \cdot t$ 的累积值），取决于保护电路的响应速度。这个积分称为MIITS（Mega I²t）。
也就是说，如果保护电路能快速衰减电流，右边就会变小，从而抑制热点温度，对吧！
没错。设计标准是，保护电路设计需满足 $T_\text{max} < 300$ K（环氧浸渍线圈）或 $T_\text{max} < 150$ K（高性能磁体）。从失超检测到能量抽取开始的延迟时间 $t_\text{det}$ 越长，MIITS越大，热点温度就越高，因此检测速度是生命线。
失超保护电路
保护电路具体是什么机制？
主要有三种方式。实际磁体中会组合使用。
(1) 能量抽取（Energy Extraction）
外部设置泄放电阻 $R_d$，检测到失超后打开开关，将电流导向泄放电阻。电流呈指数衰减：
其中 $R_q(t)$ 是随时间增大的线圈常导电阻。时间常数 $\tau = L/(R_d + R_q)$ 越小，能量抽取越快。但线圈端子间电压 $V = I \times R_d$ 不能超过绝缘耐压，因此 $R_d$ 有上限约束。LHC将最大电压限制在约1 kV。
(2) 失超加热器（Quench Heater）
在线圈外表面粘贴不锈钢箔加热器，检测到失超后通以大电流脉冲（数百A）。这使线圈整体被故意转为常导态，将储存的能量分散到整个线圈。这样可以避免局部热点。
(3) CLIQ（Coupling-Loss Induced Quench）
CERN开发的新方法，通过在线圈内部诱发交流损耗来快速扩展失超。已被LHC高亮度（HL-LHC）的内侧三重四极磁体采用。比传统的失超加热器更快、更可靠。
保护的基本原理就是“快速抽取能量”或“均匀扩散到整体”，或者两者并用，对吧。
理解得很完美。而为了优化这些保护方式的参数（$R_d$ 值、加热器延迟时间、CLIQ电容），失超传播的FEM模拟是必不可少的。
“迷你大爆炸”的守护者 — CERN的超导磁体与失超保护
LHC有1,232台超导双极磁体排列，总储存能量约达10 GJ（相当于2.4吨TNT炸药）。每个磁体区段都有独立的失超保护系统，能在10毫秒内检测到异常的电压上升并启动能量抽取。2008年事故后，CERN逐一重新测量了所有连接处的超导接合电阻，并新增了6,000个安全阀和强化型失超检测系统。这项改造花费了约4,000万瑞士法郎（约60亿日元）。“一个纳欧姆的不良导致60亿日元的修理费”，这个事实说明了失超模拟的重要性。
引发失超的机制
- 导体运动（Conductor Motion）：电磁力导致的导体微小滑动产生摩擦热。是NbTi磁体最常见的失超原因，也是“训练失超”的主因。
- 环氧树脂裂纹：浸渍树脂的微小裂纹引发局部发热。原因是极低温下环氧树脂的热收缩和应力集中。
- 束流损失（对加速器而言）：带电粒子束的一部分撞击超导线圈造成局部加热。LHC使用束流损失监测器（BLM）进行监控。
- 交流损耗：变化磁场下的磁滞损耗、耦合损耗、涡流损耗。在脉冲磁体或快速循环磁体中占主导。
- 连接处电阻：超导线材接头处的电阻性发热。LHC 2008年事故即由此引起。
量纲分析与典型参数值
| 参数 | 符号 | SI单位 | NbTi @ 4.2 K 典型值 |
|---|---|---|---|
| 临界温度 | $T_c$ | K | 9.2 K |
| 临界电流密度 | $J_c$ | A/m² | $3 \times 10^9$ (5 T, 4.2 K) |
| 铜的电阻率 | $\rho_\text{Cu}$ | Ω·m | $2 \times 10^{-10}$ (RRR=100) |
| 体积比热容 | $C$ | J/(m³·K) | $\sim 800$ (4.2 K) |
| 热导率（轴向） | $k_\parallel$ | W/(m·K) | $\sim 600$ |
| 热导率（横向） | $k_\perp$ | W/(m·K) | $\sim 1$ (含绝缘层) |
| NZPV | $v_q$ | m/s | 5〜20 |
| MIITS允许值 | — | MA²·s | 10〜30 |
数值解法与实现
FEM公式化
要用FEM求解失超的控制方程，会是什么样的公式化呢？
从热传导方程的伽辽金弱形式出发。使用试验函数 $w$：
将温度场用形函数 $N_i$ 近似并离散化为 $T^h = \sum_i N_i T_i$，最终得到以下矩阵方程：
其中：
- $\mathbf{C}$ — 热容矩阵（$C_{ij} = \int_\Omega N_i \rho c_p N_j \, d\Omega$）。具有温度依赖性
- $\mathbf{K}$ — 热传导矩阵（$K_{ij} = \int_\Omega \nabla N_i \cdot k \nabla N_j \, d\Omega$）。具有温度依赖性
- $\mathbf{F}$ — 由焦耳发热+外部加热+冷却边界条件构成的载荷向量
由于 $\mathbf{C}$、$\mathbf{K}$、$\mathbf{F}$ 都依赖于温度 $T$，因此是强非线性问题。牛顿-拉弗森法迭代求解每个时间步是标准方法。
时间积分方案
时间方向的离散化怎么做？失超是快速现象，时间步长似乎很重要。
观察得很准。失超传播的FEM中，标准使用后向欧拉法（一阶精度的隐式法）或广义梯形法则（$\theta$法）：
$\theta = 1$ 是后向欧拉（无条件稳定），$\theta = 0.5$ 是克兰克-尼科尔森（二阶精度但可能振荡）。
失超初始阶段，常导前沿移动速度快，需要将时间步长细化到 $\Delta t \sim 10^{-5}$ 秒。另一方面，在能量抽取阶段，时间尺度为 $\tau \sim 0.1$ 秒量级，可以增大 $\Delta t$。实际工作中，使用自适应时间步进，根据温度变化率自动调整步长是铁律。
时间步长$10^{-5}$秒，整体衰减是$0.1$秒量级，那岂不是需要大约$10^4$步？
如果使用自适应步进，开始时步长小（$10^{-5}$ s），温度分布平滑后逐渐增大（$10^{-3}$ s）。通常收敛需要1,000〜5,000步左右。COMSOL的BDF求解器会自动完成这个操作。
材料非线性的处理
听说极低温下材料物性变化剧烈，怎么处理呢？
这是失超模拟最大的难点。以铜的比热为例，4 K时只有约0.1 J/(kg·K)，而300 K时约385 J/(kg·K)。也就是说变化了3,000倍以上。热导率也在数百倍的尺度上变化。
实际工作中使用以下数据源：
- NIST Cryogenic Material Properties Database — 铜、铝、不锈钢等标准材料在4〜300 K范围内的 $c_p(T)$, $k(T)$, $\rho_\text{el}(T)$ 数据已整理好
- Bottura 的 NbTi 临界表面模型 — $J_c(T, B)$ 的经验公式。描述了磁场和温度的双重依赖性
- 在COMSOL中使用查找表插值（分段线性或3
なった
詳しく
報告