In previous section, a MIITs calculation is preformed to estimate the maximum temperature after the magnet quench with a given resistance of dump resistor and
Chapter 5. Quench Protection for Superconducting Magnets
material properties. However, the PCS is a special case, in which the radiation damage could affect the maximum temperature after quench significantly. With the consideration of the effect of radiation damage on the RRR of conductor, the maximum temperature estimated from the MIITs calculation becomes higher than the allowable temperature of 200 K for PCS. To investigate the radiation effect on the behavior of magnet quench, a quench simulation is preformed.
5.2.1 Overview of Quench Simulation
Figure 5.2 shows the flow chart of quench simulation. It consists of a thermal model with a detailed coil geometry as described in chapter 3 and a electrical model of quench protection circuit as given in Sec. 5.1.1. Due to massive computation the connection of coils is ignored, the quench simulation is preformed with a single coil to give a conservative study. In order to implement the irradiation effects, the profile of neutron flux is utilized to calculated the degradation of RRR, additionally, the result of coil temperature described in chapter 3 is inputted as the initial temperature. To trigger the magnet quench, a heater is fired to create a hot spot inside coil with a given time duration. Then the evolution of quench propagation is calculated as follows:
Figure 5.2: Overview of the program flow of quench simulation. The main flow of simulation is shown on left, and the right one is the flow of calculation of joule heating.
1). The material properties in each element are updated with the profiles of magnetic field, RRR and temperature.
2). The joule heating in each element of conductor is calculated with current and resistance, in which the transition of superconductor is solved with an iterative method. The details will be described in Sec. 5.2.2 later.
3). The resistance of each conductor is summed up to achieve the coil resistance for the calculation of current decay.
4). The quench detection is validated with a criteria that the coil voltage is higher than 0.1 V during a duration of 0.1 sec.
5). The current is obtained by solving the differential equation of electrical circuit with a Runge-Kutta method, then the magnetic filed is estimated from current.
6). The thermal analysis is preformed to calculate the temperature in coil.
7). The minimum time step is calculated from the thermal diffusivity and mesh size to avoid the divergence of calculation.
8). The calculation is repeated until the current is decayed to a value closed to zero.
5.2.2 Joule Heating
In order to obtain the coil temperature after quench, the heat generation term in the heat transfer equation (3.3) is replace by
Q=
(Qh(t) +Qj, conductor
Qh(t), others (5.4)
where Qh is the heat generation provided by heater,Qj is the joule heating, respectively.
Since the current only flows in the conductor, the joule heating is not allowed in the element such as support shell and aluminum strip.
Figure 5.3: The equivalent circuit at the superconducting transition state (left) and normal conducting state (right).
The heat generation of joule heating is calculated with three cases: 1). superconducting state, in which the heat generation is equal to zero, 2). transition state, in which the
Chapter 5. Quench Protection for Superconducting Magnets
equivalent circuit consists of superconductor, copper matrix and aluminum stabilizer in parallel, and 3). normal conducting state, in which the equivalent circuit only is contained the copper matrix and aluminum stabilizer as given in Fig. 5.3. In the transition state, the correlation of current and voltage is given by
Esc = Ec( Isc
Ic(B, T))n = ρCu
ACuICu= ρAl
AAlIAl, (5.5)
I = IAl+Isc+ICu, (5.6)
where Ec is the reference electric field to determine the critical current, n is the n-value, Isc is the current in superconductor,Ic is the critical current of superconductor, ρCu is the resistivity of copper,ACu is the cross-sectional area of copper, ICu is the current in copper matrix, and IAl is the current in aluminum-stabilizer, respectively. Here, the reference Ec is determined to be 0.1 µV/cm, and n-value is set to be 40 in this simulation. The current shared with the aluminum stabilizer can be calculated by coupling the equation (5.6) and (5.6) as
ρAl
AAlIAl = ρCu
ACu(I−IAl−Isc), (5.7)
then we get
IAl = ρCu/ACu
ρCu/ACu+ρAl/AAl(I−Isc), (5.8) where AAl is the cross-sectional area of aluminum stabilizer. Combining the equation (5.5) with (5.6), we obtain
f(Isc) =Ec( Isc
Ic(B, T))n− ρCu
ACu(I−IAl−Isc) = 0. (5.9) Since the resistivity and critical current can be determined by temperature and magnetic field, the current in superconductor Isc can be obtained from the function (5.9) by root-finding. The Isc is calculated with a secant method as described in Appendix A.2.
Then, the heat generation of joule heating is expressed as
Qj =
0, Ic > 0,Esc/Ec < 0.8 (I−IAl)ICuRCu+IAl2 RAl, Ic > 0,Esc/Ec > 0.8 ICu2 RCu+IAl2 RAl, Ic= 0
. (5.10)
Since the heat generation exists below the reference critical current, the iterative calculation starts when the term of Esc/Ec= (Isc/Ic)n is higher than 0.8 in simulation.
Current Redistribution in Stabilizer
For the current flowing in an aluminum-stabilizer, the mechanism is more complex.
The behavior of current redistribution in aluminum stabilizer is occurred after the quench due to the resistivity of aluminum, and current dominated in the aluminum-stabilizer is a function of time as well as a magnetic field [97], [98]. As reported in Ref. [99], [100], [101], the current redistribution changes the quench propagation velocity significantly
in a highly stabilized superconducting cables so that the conventional equation given in Ref. [102] to calculate the quench propagation velocity cannot be used at high operation current region.
A simple function is given by A.V. Garilin et al. to express the behavior of current redistribution [103] in the quench simulation for the ATLAS central solenoid, however, the expression is not mentioned in detail. Thus the Garilin’s function to model the current distribution is modified with a calculation from the Maxwell equation.
Ignored the effect of magnetic field, the current redistribution can be obtained from the following Maxwell equations:
1
µ∇ ×B = J +∂E
∂t , (5.11)
∇ ×E = −∂B
∂t , (5.12)
where B is the magnetic field,J is the current density,E is the electrical field,µ is the magnetic permeability, and is the dielectric constant, respectively. In a conductor, the second term of equation (5.12) can be ignored. Taking the curl of equation (5.12), then we obtain
∇ × ∇ ×E =−∂(∇ ×B)
∂t (5.13)
Inputing the equation (5.12) and the Ohm’s law: J = σE, the equation can be written as
∇ × ∇ ×E+σµ∂E
∂t = 0 (5.14)
where σ is the electrical conductivity. Due to the differential relations [104], the term with the curls is given by
∇ ×(∇ ×E) = ∇(∇ ·E)−(∇ · ∇)E (5.15) Since the magnetic permeability is equal to the vacuum permeability of 4π×10−7 H/m in the aluminum stabilizer, when ∇ ·E= 0, this equation can be written as
− ∇2E+σµ0
∂E
∂t = 0 (5.16)
Applying the relation of electrical resistivity, ρ = 1/σ and the Ohm’s law, the govern equation of current diffusion can be written as
∂J
∂t = ρ
µ0∇2J (5.17)
in which the magnetic diffusivity is defined as Dm = ρ
µ0
. (5.18)
Chapter 5. Quench Protection for Superconducting Magnets
Compared with the heat transfer equation, the current diffusion is similar to the heat diffusion. Thereby, according to the Ref. [49], the exact solution of this equation for a slab of width 2a is found
∆J =
∞
X
n=1
Cncos(nπx
2a )exp(−Dmn2π2
4a2 t) (5.19)
where t is the time, ∆J is the change of current density with respect to the initial current density, and Cn is determined from the initial condition and following equation:
Cn= 1 a
Z 2a 0
Jicos(nπx
2a )dx= 4
nπJi, n = 1,3,5, ... (5.20)
Diffusion Region
dAl
lsc
wsc
l
w
Figure 5.4: Current diffusion in Al-stabilized conductor. The current diffusion region with the width of dAl starts from the region of Rutherford cablelsc×wsc until fill of the entire stabilizer l×w.
From the equation (5.19), the current changes with a characteristic time τm = 4a2µ0
π2ρ . (5.21)
So that the diffusion length dAl as shown in Fig. 5.4 can be expressed as follows by assuming the diffusion length is linear with the change of time.
dAl(t) = s
π2ρAl(t−tcs)
4µ0 (5.22)
where tcs is the time at the moment when the current sharing is occurred,t is the time, and ρAl is the electrical resistivity of aluminum, respectively. Thus the cross-sectional area of aluminum stabilizer in equation (5.8) and (5.9) is calculated as
AAl = (wsc+ 2dAl)(lsc+ 2dAl)−ACu−AN bT i, dAl ≤ l−lsc
2 and dAl ≤ w−wsc
2 (5.23) where l×w is the cross section of 15 × 4.73 mm2 for the conductor, andlsc ×wsc is the cross section of 8.05 × 2.4 mm2 for Rutherford cable.
5.2.3 Electrical Circuit Model
As shown in Fig. 1.6, the differential equation for quench protection circuit can be written as
LtotdI
dt +I(Rdump+Rcoil) +Vd = 0 (5.24) where Vdis the diode voltage. By summing up the resistance in each element of conductor, the coil resistance Rcoil is calculated as
Rcoil =
(0, Esc/Ec<0.8
RAlRCu
RAl+RCu, Esc/Ec>0.8 (5.25) in which the coil resistance not only is a temperature, RRR and magnetic field dependent parameter but also relies on the radiation damage in this case. The voltage of diode is implemented as a function of current provided from the Shockley diode equation,
Vd(I) =Vtln(I
Is + 1) (5.26)
where Is is the saturation current, and Vt is the thermal voltage, respectively. A diode with a thermal voltage of 0.153 V and a saturation current of 200 mA is used in this simulation. To obtain the current in each time step, this differential equation is solved numerically with a iterative method, the fourth-order Runge-Kutta method, which is introduced in Appendix A.1.
The effect of the magnetization is ignored, and the change of magnetic field is obtained by assuming the magnetic field changes as same as current,
Bt=Bt−1(1 + ∆I
It−1) (5.27)
where Bt−1 is the magnetic field at the previous time step, and ∆I = It−It−1 is the difference of the current between present and previous time step, respectively.
Total Inductance in Coil
During the process of magnet quench, the inductive voltage is changed by the current decay, which could cause a higher voltage drop if the current changes faster. In order to investigate the change of voltage after quench, because only the vector potential along the azimuth exists in a solenoid the profile of coil total inductance is calculated from the result of magnetic field simulation as
L= Φ I0 = 1
I0 Z
Aθdlθ (5.28)
where Φ is the magnetic flux, I0 is the operating current of 2700 A, Aθ is the vector potential along the azimuth andl is the length along the azimuth of solenoid, respectively.
The vector potential is derived from the magnetic field simulation. Then the total inductance at each position (i, j) can be obtained from
Li,j = 2πri,jAi,j
I0 (5.29)
Chapter 5. Quench Protection for Superconducting Magnets
where the ri,j is the radius. Figure 5.5 shows the calculated coil total inductance. The difference between the maximum and minimum of total inductance in each coil is within 1.3 mH, and the peak of total inductance is about 3.3 mH located in CS1 coil.
1.02 1.00 0.98 0.96 0.94 0.92 0.90 0.88 Z [m]
0.68 0.70 0.72 0.74 0.76 0.78 0.80 0.82
R [m]
1.65
1.80 1.95
2.10 2.25
1.51 1.601.69 1.78 1.87 1.96 2.05 2.14 2.23 2.32
Inductance [mH]
0.6 0.4 0.2 0.0 0.2 0.4
Z [m]
0.68 0.70 0.72 0.74 0.76 0.78 0.80 0.82
R [m]
2.25 2.50
2.50 2.75
3.00 3.25
2.055 2.190 2.325 2.460 2.595 2.730 2.865 3.000 3.135 3.270
Inductance [mH]
0.8 1.0 1.2 1.4 1.6 1.8 2.0
Z [m]
0.68 0.69 0.70 0.71 0.72 0.73 0.74 0.75
R [m]
1.80 1.95
2.10 2.25 2.40
1.72 1.80 1.88 1.96 2.04 2.12 2.20 2.28 2.36 2.44
Inductance [mH]
2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9
Z [m]
0.68 0.70 0.72 0.74 0.76 0.78
R [m]
1.40 1.80 1.60 2.00
1.25 1.35 1.45 1.55 1.65 1.75 1.85 1.95 2.05 2.15
Inductance [mH]
Figure 5.5: Estimated total inductance in CS0 (upper left), CS1 (upper right), MS1 (lower left) and MS2 (lower right).