Superconducting Quench Simulation

Good question. Quench is a phenomenon where a part of a superconductor transitions to a normal conducting state by exceeding its critical temperature $T_c$. In the superconducting state, electrical resistance is zero, so current flows without loss. However, if even one point transitions to the normal state, electrical resistance suddenly appears at that part.

Exactly. Moreover, large currents flow in superconducting magnets. Taking the LHC superconducting dipole magnet as an example, a current of about 11,850 A flows to generate an 8 Tesla magnetic field. The inductance of one coil is about 100 mH, so the stored magnetic field energy is

This is equivalent to the kinetic energy of a passenger car traveling at 160 km/h. When a quench occurs, this energy is converted into heat inside the magnet. If not controlled, the magnet could melt, or liquid helium could rapidly vaporize causing pipes to rupture.

The accident in LHC Sector 3-4 in September 2008 is famous. A poor superconducting joint (resistance value only about 220 nΩ) triggered the quench. The quench protection system did not function as designed in some sections, and an electrical arc discharge punctured the helium vessel. About 6 tons of liquid helium vaporized instantly, damaging 53 dipole magnets and surrounding infrastructure. Ultimately, repairs and safety reinforcement took 14 months.

Exactly. The objectives of quench simulation are threefold:

1. Hotspot Temperature Prediction — Verify that the point with the highest temperature rise does not exceed the magnet material's allowable temperature (typically 300–400 K)
2. Quench Detection Time Determination — Quantify the allowable delay time before the protection system reacts
3. Protection Circuit Parameter Optimization — Determine dump resistance values, quench heater placement, and heating power

The basis is the heat conduction equation with a Joule heating source term. The temperature distribution $T(\mathbf{x}, t)$ of the superconducting wire follows this equation:

The meaning of each term is as follows:

• $\rho \, c_p$ — Volumetric heat capacity [J/(m³·K)]. Superconducting wires are composites (NbTi + Cu + epoxy), so effective values are used.
• $k(T)$ — Thermal conductivity [W/(m·K)]. At cryogenic temperatures, temperature dependence is extreme (copper is about 600 W/(m·K) at 4 K, about 2000 W/(m·K) at 30 K, and about 400 W/(m·K) at 300 K).
• $\rho_\text{el}(T, B)$ — Electrical resistivity [Ω·m]. Zero in the superconducting state ($T < T_c$), equal to the copper matrix resistivity in the normal state. There is a current sharing transition region near $T_c$.
• $J$ — Current density [A/m²]
• $Q_\text{ext}$ — External heating source (e.g., quench heaters) [W/m³]

Sharp observation. To be precise, superconductors have another critical temperature called the current sharing temperature $T_\text{cs}$. For $T < T_\text{cs}$, all current flows through the superconducting filaments, generating no heat. For $T_\text{cs} < T < T_c$, part of the current diverts to the copper matrix, and heating begins. For $T > T_c$, it becomes fully normal conducting. Accurately modeling this transition region is key to simulation accuracy.

Here, $T_\text{op}$ is the operating temperature, $I_\text{op}$ is the operating current, and $I_c$ is the critical current.

The Normal Zone Propagation Velocity (NZPV) is the most important parameter representing the speed of quench spread. Assuming one-dimensional steady-state propagation, it can be approximated by Wilson's classical formula:

Here, $C(T_c)$ is the volumetric heat capacity at the critical temperature.

It varies dramatically by material. This is a core issue in superconducting magnet design.

Note that REBCO (High-Temperature Superconductor) NZPV is only 1/100 to 1/1000 that of NbTi. A slow NZPV means generated heat is confined to a narrow region, causing hotspot temperature to rise rapidly. This is why HTS magnet protection design is extremely difficult.

The adiabatic assumption is used to estimate the upper limit of hotspot temperature. As a worst-case scenario, heat conduction from the quench origin is assumed to be zero. Then the heat balance is:

We start from the Galerkin weak form of the heat conduction equation. Using a test function $w$:

Approximating the temperature field with shape functions $N_i$ as $T^h = \sum_i N_i T_i$ and discretizing, we finally get this matrix equation:

Where:

• $\mathbf{C}$ — Heat capacity matrix ($C_{ij} = \int_\Omega N_i \rho c_p N_j \, d\Omega$). Temperature dependent.
• $\mathbf{K}$ — Thermal conductivity matrix ($K_{ij} = \int_\Omega \nabla N_i \cdot k \nabla N_j \, d\Omega$). Temperature dependent.
• $\mathbf{F}$ — Load vector consisting of Joule heating + external heating + cooling boundary conditions.

Since $\mathbf{C}$, $\mathbf{K}$, and $\mathbf{F}$ all depend on temperature $T$, it becomes a strongly nonlinear problem. The standard approach is to solve each time step iteratively using the Newton-Raphson method.

Good point. In FEM for quench propagation, the Backward Euler method (first-order accurate implicit method) or the generalized trapezoidal rule ($\theta$-method) is standardly used:

$\theta = 1$ is Backward Euler (unconditionally stable), $\theta = 0.5$ is Crank-Nicolson (second-order accurate but may oscillate).

In the initial stage of a quench, the normal front moves rapidly, so time steps need to be as fine as $\Delta t \sim 10^{-5}$ seconds. On the other hand, during the energy extraction phase, the time constant is on the order of $\tau \sim 0.1$ seconds, so $\Delta t$ can be larger. In practice, using adaptive time stepping to automatically adjust step width based on the rate of temperature change is the golden rule.

With adaptive stepping, steps are fine initially ($10^{-5}$ s) and gradually become coarser ($10^{-3}$ s) as the temperature distribution smooths out. Typically, it converges within about 1,000–5,000 steps. COMSOL's BDF solver does this automatically.

This is the biggest challenge in quench simulation. Taking copper's specific heat as an example, at 4 K it's only about 0.1 J/(kg·K), but at 300 K it's about 385 J/(kg·K). That's a change of over 3,000 times. Thermal conductivity also changes by a factor of several hundred.

In practice, the following data sources are used:

• NIST Cryogenic Material Properties Database — $c_p(T)$, $k(T)$, $\rho_\text{el}(T)$ for standard materials like copper, aluminum, stainless steel are compiled for 4–300 K.
• Bottura's NbTi Critical Surface Model — Empirical formula for $J_c(T, B)$. Describes dependence on both magnetic field and temperature.
• In COMSOL, lookup table interpolation (piecewise linear or 3