- Email: [email protected]

Contents lists available at ScienceDirect

Cryogenics journal homepage: www.elsevier.com/locate/cryogenics

Research paper

Multiplicities and thermal runaway of current leads for superconducting magnets Rizos N. Krikkis Institute of Thermal Research, 2 Kanigos Str, PO Box 106 77, Athens, Greece

a r t i c l e

i n f o

Article history: Received 31 October 2016 Received in revised form 29 December 2016 Accepted 23 January 2017 Available online 9 February 2017 Keywords: Current lead Temperature blow-up thermal runaway Superconducting magnet Numerical bifurcation analysis Multiplicity

a b s t r a c t The multiple solutions of conduction and vapor cooled copper leads modeling current delivery to a superconducting magnet have been numerically calculated. Both ideal convection and convection with a finite heat transfer coefficient for an imposed coolant mass flow rate have been considered. Because of the nonlinearities introduced by the temperature dependent material properties, two solutions exist, one stable and one unstable regardless of the cooling method. The limit points separating the stable form the unstable steady states form the blow-up threshold beyond which, any further increase in the operating current results in a thermal runway. An interesting finding is that the multiplicity persists even when the cold end temperature is raised above the liquid nitrogen temperature. The effect of various parameters such as the residual resistivity ratio, the overcurrent and the variable conductor cross section on the bifurcation structure and their stabilization effect on the blow-up threshold is also evaluated. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction Superconducting magnets maintained at cryogenic temperatures are powered from current leads at room temperature, Wilson [1], Iwasa [2]. The thermal connection between the room and cryogenic environments through the current leads introduces a significant heat leak to the cryostat and a substantial theoretical and experimental work has been carried out in the effort to minimize this heat leak to the cryogenic liquid and optimize the lead design with respect to material and geometry selection [3–8]. The basic configurations examined consist of vapor and conduction cooled (or ‘dry’) leads. An interesting feature of the mechanism of heat generation produced by an applied current and its balance by conduction is the multiple steady states. Jones et al. [9] employed temperature and purity dependent material properties to determine burn-out limits for copper current leads cooled by helium vapor. During numerical overcurrent simulations two steady states were calculated, in certain cases depending on the initial conditions and on the rate the current was increased and reached its final value. Aharonian et al. [10] and Buchs and Hyman [11] demonstrated that for a fixed current, geometry and vapor flow rate two solutions exist. The stability of the two solutions was assessed numerically i.e. by introducing a small perturbation on the temperature profile and calculating the response in the time domain. A solution resulting in thermal runaway was considered as unstable. More recently E-mail address: [email protected] http://dx.doi.org/10.1016/j.cryogenics.2017.01.009 0011-2275/Ó 2017 Elsevier Ltd. All rights reserved.

Hanzelka [12] reported that for certain combinations of copper lead geometry and applied current operating in vacuum with radiative exchange, as many as three solutions were found. As it is obvious from the literature cited in the preceding paragraphs most of the work has been focused on the design and optimization of current leads and only a few studies considered the multiple solutions and the stability aspects of the problem. As a matter of fact Hull [13] argued that the multiplicity may be an artifact of the numerical methods employed. Today, the field of numerical bifurcation analysis has become mature and bifurcation mechanisms are widely accepted as decisive phenomena for explaining and understanding stability and structural change [14–16]. Within this framework the aim of the present study is to numerically explore the multiplicity and blow-up (thermal runaway) features of copper current leads delivering current from a relatively warm environment to superconducting magnets at cryogenic temperatures. Multiple solutions exist not only for the vapor cooled leads but for conduction cooled leads as well. The solution structure is analyzed with sufficient bifurcation diagrams describing the effects of the residual resistivity ratio, the conductor geometry, the cold end temperature and the overcurrent on the multiplicity regions and the blow-up threshold. 2. Analysis Consider a copper conductor of cylindrical geometry with a variable cross section AðXÞ, length L, thermal conductivity K, elec^ and specific heat C, as it is schematically trical resistivity q

9

R.N. Krikkis / Cryogenics 83 (2017) 8–16

Nomenclature A c C F G h H I k K L _ m p P q Q RRR t T u v V x X

conductor cross sectional area [m2] ðC=C ref Þ reduced specific heat capacity [–] conductor specific heat capacity [J/(kg K)] flow number, Eq. (9) [–] generation number, Eq. (8) [–] ðH=Href Þ reduced heat transfer coefficient [–] heat transfer coefficient [W/(m2 K)] current through lead [A] ðK=K ref Þ reduced thermal conductivity [–] conductor thermal conductivity [W/(m K)] conductor length, Fig. 1 [m] coolant mass flow rate [kg/s] profile tapper ratio, Eq. (3) [–] wetted perimeter [m] heat load, Eq. (14) [W] reduced heat load, Eq. (15) [–] residual resistivity ratio [–] time [sec] temperature [K] conduction-convection parameter (CCP) [–] ðV=V ref Þ dimensionless voltage [–] voltage difference across the lead [V] ðX=LÞ dimensionless distance along conductor, Fig. 1 [–] distance along conductor, Fig. 1 [m]

depicted in Fig. 1. The warm end is maintained at ambient temperature, say T H ¼ 300 K and the cold end at liquid helium temperature T L ¼ 4:2 K. A helium gas stream of constant mass flow rate _ is used to cool the conductor. Assuming that the conductor is m thermally thin so that transverse temperature gradients may be neglected the energy balance for the lead and the cooling gas take the form:

X, x

TH , AH

y z

dimensionless transverse coordinate, Fig. 1 [–] ðI=Iref Þ2 current overload factor [–]

Greek symbols a thermal diffusivity [m2/s] b ðccp AÞg =ðcC ref AH Þ time scaling factor [–] c density [kg/m3] H ðT=DT ref Þ dimensionless temperature [–] k eigenvalue [–] ^ =q ^ ref Þ reduced conductor electrical resistivity [–] q ðq q^ conductor electrical resistivity [X m] s ðat=L2 Þ dimensionless time [–] Subscripts g gas coolant H warm end ðx ¼ 1Þ of lead L cold end ðx ¼ 0Þ of lead LP reference to limit points ref reference value s reference to steady state Superscript ð0Þ derivative with respect to x

cCðTÞA

@T @ @T q^ ðTÞ ¼ KðTÞA HPðT T g Þ þ I2 ; @t @X @X A

½cðT g Þcp ðT g ÞAg

@T g @T _ p ðT g Þ g HPðT T g Þ; ¼ mc @t @X

ð1Þ

ð2Þ

where T is the conductor temperature, c its density, cp is the coolant specific heat capacity, T g its temperature, H is the convective heat transfer coefficient and I is the applied current. The design of the conductor cross section A has been a subject of optimization as well, either with simple geometries, Eckert et al. [17], Jiahui et al. [18] or as a variational problem, Okolotin and Bol’shakov [19]. For the purposes of the present study and in order to reduce the number of the parameters involved in the analysis a simple linear profile has been considered

yðXÞ ¼ p þ ð1 pÞðX=LÞ;

p P 1:

ð3Þ

In terms of the above profile the cross sectional area and the wetted perimeter may be expressed as

I

AðXÞ ¼ AH y2 ðXÞ;

PðXÞ ¼ PH yðXÞ:

ð4Þ

Introducing dimensionless variables

L

x ¼ X=L; H ¼ T=DT ref ; Hg ¼ T g =DT ref ;

s ¼ at=L2 ; y2 ¼ A=AH

^ =q ^ ref ; c ¼ C=C ref ; z ¼ ðI=Iref Þ2 h ¼ H=Href ; k ¼ K=K ref ; q ¼ q

m

;

ð5Þ the partial differential equations describing the temperature distribution of the conductor and the cooling gas take the form:

TL , AL

liquid

y

cy2

b

@H @ zGq 2 @H u2 yhðH Hg Þ þ 2 ; ¼ ky @ s @x @x y

@ Hg @ Hg ¼ u2 hðH Hg Þ F ; @s @x

ð6Þ

ð7Þ

where b ¼ ðccp AÞg =ðcC ref AH Þ is a time scaling factor. The current Fig. 1. Physical model and coordinate system.

overload factor z is introduced since the leads may not always

10

R.N. Krikkis / Cryogenics 83 (2017) 8–16

operate under steady state conditions, Chang and Kim [20] or the nominal current may be exceeded under fault conditions Xu et al. [21]. Three important dimensionless numbers appear in the above relationships, namely the generation number G the flow number F and the conduction-convection parameter u. The generation number is defined as:

G¼

L AH

2

^ I q K DT 2

!

:

ð8Þ

ref

Physically it provides a measure of the ratio of the Joule heating to the heat dissipated by conduction. The flow number is defined as:

_ pL _ p DT ref mc mc F¼ ¼ K ref AH ðK DTÞref AH =L

ð9Þ

and it measures the ratio of the convective to conductive cooling, whereas the conduction-convection parameter

u2 ¼

Href L2 ; K ref ðA=PÞH

ð10Þ

is extensively used in extended surfaces and conjugate heat transfer problems [22–25] and strongly affects the solution. As it will become evident in the next sections all three numbers have a profound effect on the bifurcation structure. Under steady state conditions the partial differential equations Eqs. (6) and (7) reduce to a system of ordinary differential equations for the conductor and gas temperatures HðxÞ and Hg ðxÞ: 0

ðky H0 Þ ¼ u2 yhðH Hg Þ 2

zGq ; y2

F H0g u2 hðH Hg Þ ¼ 0; 0

ð11Þ ð12Þ

0 g

where H and H denote differentiation with respect to x. The boundary conditions at the cold and warm ends are:

Hð0Þ ¼ HL ;

Hð1Þ ¼ HH ;

Hg ð0Þ ¼ HgL :

ð13Þ

Of particular interest is the heat leak transmitted to the cryogenic liquid

dT AH 2 q ¼ KA ¼ p ðK DTÞref ðkH0 ÞL : dX X¼0 L

ð14Þ

Substituting the ratio AH =L from the definition of the generation number in Eq. (8), the dimensionless heat load reads:

Q¼

q 1=2

^ K DTÞref p2 ðI2 q

ðkH0 Þ ¼ pﬃﬃﬃﬃ L : G

ð15Þ

0

ðk#0 Þ þ kH H0 #0 þ ðkHH H0s þ kH H00s þ zGqH ckÞ# u2 ð# wÞ ¼ 0; 2

ð18Þ Fw0 þ bkw u2 ð# wÞ ¼ 0;

ð19Þ

where c ¼ cðHs Þ, k ¼ kðHs Þ, kH ¼ ð@[email protected] HÞHs ðxÞ , kHH ¼ ð@ 2 [email protected] H2 ÞHs ðxÞ and primes denote differentiation with respect to x. The corresponding boundary conditions are:

#ð0Þ ¼ #ð1Þ ¼ wð0Þ ¼ 0:

ð20Þ

If all eigenvalues are negative then the steady state solution under consideration is stable (and denoted with a continuous line on the bifurcation diagrams). If on the other hand at least one eigenvalue is positive, the steady state solution is unstable (denoted with a dashed line). For the numerical solution of Eqs.(18), (19) the steady state solutions Hs ðxÞ and Hgs ðxÞ must be available, so Eqs. (11), (12) are attached to Eqs. (18), (19) forming an extended boundary value problem. The eigenfunctions were normalized using the condition #0 ð0Þ ¼ 1, Seydel [14]. 2.2. Material properties The material properties of copper are considered as temperature dependent whereas the thermal conductivity and the electrical resistivity depend additionally on the RRR. Semi-empirical correlations that fit the experimental measurements with reasonable accuracy have been employed for the thermal conductivity and the electrical resistivity as proposed by Hust and Lankford [27]. For the copper specific heat the correlation provided by Jones ^ on et al. [9] was used. The nonlinear dependence of the K and q temperature and RRR is depicted in Fig. 2 and in Fig. 3 respectively. In this respect the use of the Wiedemann-Franz law is avoided. 2.3. Numerical methods The two-point boundary value problem described by Eqs. (11), (12) has been solved numerically. In order to ensure an accurate numerical solution, two different methods have been employed, resulting in identical results under the strict tolerances imposed. The first one utilizes a multi-shooting Runge-Kutta formula pair of order 8(7), Hairer et al. [28] and the second one a spline collocation method described by Ascher et al. [29]. Continuation along the various branches has been carried out along the lines suggested by Seydel [14]. For the computation of the singular points an extended problem is being formed from the partial derivatives of Eqs. (11),

As it will be shown in Section 3, Eqs. (11), (12) have two solutions and an assessment of the stability is essential since only stable solutions are physically realizable. The stability of a certain steady states Hs ðxÞ, Hgs ðxÞ to small perturbations, #ðxÞ, wðxÞ i.e.:

Hðx; sÞ ¼ Hs ðxÞ þ #ðxÞ expðksÞ;

ð16Þ

Hg ðx; sÞ ¼ Hgs ðxÞ þ wðxÞ expðksÞ;

ð17Þ

will be determined by the eigenvalues k of the corresponding linearized problem with respect to the steady state, Michelsen and Villadsen [26]. Substituting Eqs. (16) and (17) into Eqs. (6) and (7) and considering a uniform cross section ðy ¼ 1Þ and a constant heat transfer coefficient ðh ¼ 1Þ for simplicity, the eigenvalue problem describing the stability of the steady states for the lead and the gas read:

thermal conductivity [W/(mK)]

2.1. Stability

10

00 30

4

00 10 10

0 30

3

0 10

10

R= RR

2

1

30

10

100

temperature [K] Fig. 2. Copper thermal conductivity with RRR as a parameter.

11

R.N. Krikkis / Cryogenics 83 (2017) 8–16

300 p = 1, z = 1

-8

10

stable unstable limit points

blow-up threshold

10-9

(dΘ/dx)L

electrical resistivity [Ωm]

200

RRR = 20 40

40

60

10-10

80

30

110

10

100

120 RRR = 160

20

temperature [K]

thermal runaway

100 90 80 70 60 40 50 60

10

15

20

25

30

35

40

45

50

55

60

generation number G

Fig. 3. Copper electrical resistivity with RRR as a parameter.

Fig. 4. Blow-up threshold on the ðG; H0L Þ plane with RRR as parameter.

(12) with respect to the parameters, according to Witmer et al. [30]. The solution of the nonlinear boundary value problem arising in the stability section, Eqs.(18), (19), is carried out iteratively. For the ideal convection case, where only Eq. (18) is involved, the difficulty of obtaining starting values for the eigenvalues is tackled by neglecting the #0 term in the left hand side of Eq. (18) and assuming further, a constant value for the thermal conductivity. This simplification reduces the eigenproblem of Eq. (18) to a Sturm-Liouville one, where eigenvalues and eigenfunctions may be computed easily to any specified index, order and accuracy, Pruess and Fulton [31]. The eigenvalues calculated in this way are used as starting values for the complete problem in Eq. (18) and the iterative shooting algorithm converges within a few iterations. For the case of a finite heat transfer coefficient where both Eqs.(18), (19) have to be solved the eigenvalues calculated for the ideal convection are being used as starting values. Blow-up solutions of Eqs. (6), (7) are accompanied with very abrupt changes and steep gradients and must be accurately reproduced in numerical simulations. Thus a moving mesh algorithm based on the equidistribution principle has been adopted, Budd et al. [32] and Huang et al. [33]. It is worth mentioning that although a constant heat transfer coefficient (i.e. h ¼ 1) will be used for comparing the theoretical model with experimental data in Section 3.4, the methods and the code developed apply to a general heat transfer coefficient of the form: h ¼ hðx; h; hg Þ. 3. Results and discussion 3.1. Conduction cooled leads This is the simplest case where u ¼ F ¼ 0 and Eq. (11) is reduced to: 0

ðky H0 Þ þ 2

zGq ¼ 0: y2

ð21Þ

As it can be seen from the above relationship, the conductor temperature HðxÞ will be a function of the generation number G and of the overcurrent z. It will additionally depend on the material properties, RRR and on the cross section geometry through the taper ratio p. Because of the nonlinearities involved in Eq. (21), the corresponding bifurcation diagram H0L versus G consists of multivalued curves, which are depicted in Fig. 4 for a uniform cross section conductor and a range of values of the RRR. More specifically for a given value of the generation number two solutions exist. One stable, indicated with a continuous curve and one unstable, marked with a dashed line. The source of the multiplicity is the nonlinear dependence of electrical resistivity on temperature. For a copper conduc-

tor with a typical purity of RRR ¼ 100 the thermal conductivity variations are of one order of magnitude as it can be seen from Fig. 2. In contrast the variations in the electrical resistivity are of two orders of magnitude (Fig. 3) and for this reason they are the main contribution to the multiplicity structure observed, Gurevich and Mints [34]. Sample stability results are summarized in Table 1 where several eigenvalues have been calculated for various RRR’s and a positive eigenvalue is associated with the instability of the upper branch of the solutions. Each curve exhibits a singular point and in particular a limit point GLP , defined by the relationship:

dG=dH0L ¼ 0;

ð22Þ

which separates the stable from the unstable solution branches. As G approaches the limit point GLP the eigenvalues of the stable and the unstable solutions tend to zero. It is important to mention that beyond the limit point that is for G > GLP no solution exists. This is of paramount practical importance because if the operating current during a transient or a fault exceeds the limit point, Eqs. (6), (7) can exhibit unbounded growth in finite time, that is thermal runaway or temperature blow-up, as it will be shown later (see Fig. 9). The phenomenon has been described qualitatively by Maddock et al. [35], Gurevich and Mints [36] and rigorously by Bandle and Brunner [37]. For instance, such an event is believed to be the reason for the electric faults in the 13 kA circuits of the Large Hadron Collider back in 2008, [38,39]. Therefore the line connecting the limit points is also the threshold for thermal runaway. A similar structure is observed for the relationship between the heat load Q and the generation number G shown in Fig. 5 together with the points where Q is being minimized (optimized). The calculated warm end temperature gradient is presented in Fig. 6 where the optimum points are shown as well. For this problem it is well known that for the optimization the condition H0H ¼ 0 must hold and this is an indirect verification of the validity of the bifurcation analysis, since during continuation all solutions should be calculated for the parameters range under consideration. The temperature distribution corresponding to stable and unstable states are shown in Fig. 7 for various generation numbers and a constant cross section conductor. For this particular case Gopt ﬃ 24. Below Gopt the stable temperature profiles are monotonic whereas above Gopt there exist a peak closer to the warm end, i.e. HH is exceeded locally. For the unstable solutions all temperature profiles are non-monotonic. Fig. 8 shows the effect of the cross section profile on the overcurrent margin. The maximum overcurrent prior to thermal runaway for each curve is determined by the corresponding limit point zLP . A variable cross section is accompanied with a substantial improvement of the blow-up margin since the

12

R.N. Krikkis / Cryogenics 83 (2017) 8–16

Table 1 Eigenvalues ki for the steady states with RRR ¼ 40, z ¼ 1 and constant cross section. G ¼ 10

i

1 2 3 4 5

G ¼ 24

G ¼ 35

Stable ss

Unstable ss

Stable ss

Unstable ss

Stable ss

Unstable ss

3.77 19.70 46.94 83.40 129.06

+15.08 50.74 84.99 127.40 177.96

2.54 16.87 42.37 77.70 121.82

+6.44 55.69 91.871 136.62 189.89

1.19 13.98 36.88 69.71 111.91

+1.83 31.42 60.82 99.28 146.68

100 80 60

30 limit points minimum points

stable unstable

10

24 35

40

20 p = 1, z = 1

35

20

temperature Θ(x)

heat load Q

24

10 9 8 7

RRR = 40

120 160

60 80

6

10 8 6

G = 10

4 2

stable unstable

1 0.8 0.6

blow-up threshold

5

RRR = 40, p = 1.0, z = 1.0

0.4

4 10

0.0

15

20

25

30

35

40

45

50

0.1

0.2

0.3

55

0.4

0.5

0.6

0.7

0.8

0.9

1.0

distance along conductor x

generation number G

Fig. 7. Temperature variation along the conductor length. Fig. 5. Blow-up threshold on the ðG; Q Þ plane with RRR as parameter.

30 40

RRR = 40

p = 1, z = 1

20

15

0

-40

(dΘ/dx)H

40

60

80

RRR = 160 120

24

-60 -80

blow-up threshold

-100 -120

G = 10 15 20 24 blow-up threshold

4

p=1

2

-160 limit points minimum points

-180 -200 10

10 9 8 7 6 5

3

stable unstable

-140

blow-up threshold

20

heat load Q

-20

stable unstable limit points

G = 10

20

p=2

1

2

3

4

5

6 7 8 9 10

20

2

15

20

25

30

35

40

45

50

55

generation number G Fig. 6. Blow-up threshold on the ðG; H0H Þ plane with RRR as parameter.

increased area can sustain a higher current. The blow-up phenomenon for conduction cooled leads is presented in Fig. 9. Starting form a stable steady state Hs where G ¼ 35, p ¼ 1 and z ¼ 1 the current is increased in a step-like transient (see insert in Fig. 9). From Fig. 8 the blow-up threshold may be estimated as zLP ﬃ 1:07. The temperature gradually increases around the point x ﬃ 0:7 (blow-up point) until roughly s ¼ 0:54. At s ¼ 0:55 the temperature has dramatically increased around the blow-up point whereas at s ¼ 0:55017 the blow-up is clearly manifested, since the heat dissipated by conduction can no longer match the heat generation rate produced by the applied current. Another salient feature of Eq. (21) shown in Fig. 10, is that the multiplicity is persisting even when the cold end temperature is being raised above the liquid nitrogen temperature. In the temper-

current overload factor z = (I/Iref)

Fig. 8. Blow-up threshold on the ðz; Q Þ plane with p as parameter.

ature range 100–300 K the thermal conductivity is almost constant but the electrical resistivity increases by almost an order of magnitude. Thus the nonlinear dependence of the electrical resistivity on temperature is again the source of the multiplicity and the blowup behavior. This is particularly important when the metal is coupled to a High Temperature Superconducting (HTS) material to form a composite lead either in the full superconducting mode or in the current sharing one [40–46]. HTS leads are usually embedded in a metal matrix so part or the whole transport current may be carried by the metal during a superconducting-normal transition. In such a case a blow-up phenomenon similar to the one described for metal leads may take place. Indeed the hot spot curves calculated by Wesche and Fuchs [41] in their Fig. 6, simulating a complete loss of coolant for composite HTS leads (Bi-2212

13

R.N. Krikkis / Cryogenics 83 (2017) 8–16

30

180

140 120 100

0.55017

5

20

4 2

zLP

1

0.2 0.3 0.4 0.5 0.6 0.7 0.8 dimensionless time τ

80

0.54 0.53

60

10 9 8 7 6 5

blow-up threshold

0.0

4

0.2

3

0.4

0.52

40

0.6 Θs

τ = 0.50

20

2

0.8 1.0

RRR = 40, G = 35, p = 1

0 0.0

0.1

stable unstable

p = 1, z = 1, RRR = 40

0.55

3

heat load Q

overcurrent z(τ)

160

conductor temperature Θ(x)

limit points minimum points

6

0.2

0.3

0.4

0.5

0.6

0.7

0.8

F = 1.2

0.9

1.0

1 10

distance along conductor x

30

40

50

limit points minimum points

10 9 8 7

stable unstable

limit points minimum points

heat load Q

heat load Q

5

100 60

4

blow-up threshold

3

40

100

2

77

blow-up threshold

6

60

40

30

20 TL = 4.2K

5

stable unstable p = 1, z = 1 F = 1, RRR = 40

6

RRR = 40, p = 1, z = 1

77

70 80 90

Fig. 12. Blow-up threshold on the ðG; Q Þ plane with F as parameter.

30

10 9 8 7

60

generation number G

Fig. 9. Thermal runaway for conduction cooled leads during current step-like change.

20

20

20 TL = 4.2K

4 8

10

12

14

16

18

20

22

24

26

28

30

32

34

36

38

1 10

generation number G

15

20

25

30

35

40

45

50

55

60

65

70

generation number G

Fig. 10. The effect of the cold end temperature on the solution structure. Fig. 13. Blow-up threshold on the ðG; Q Þ plane with T L as parameter. 30

limit points minimum points

20

180

stable unstable

3.0

3 2

1 0.9 0.8 0.7 0.6

10

140

conductor temperature Θ (x)

heat load Q

10 9 8 7 6 5 4

p = 1.0 1.5 2.0 2.5

120 100 80 60

overcurrent z(τ )

160 RRR = 40

2.5 2.0

zLP

1.5

0.88407

1.0

0.88

0.3 0.4 0.5 0.6 0.7 0.8 0.9 dimensionless time τ

30

40

50 60 70 8090 100

200

τ = 0.75

40 Θs

300

generation number G Fig. 11. Effect of cross section profile on the heat load Q.

0.86

0.85

RRR = 40, G = 30 F = 0.4, p = 1

20 20

0.87

0 0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

distance along conductor x Fig. 14. Thermal runaway for vapor cooled leads during current step-like change.

bulk material and Bi-2223/Ag tapes), are very similar to the blowup profiles for the conduction cooled copper leads in Fig. 9 and the vapor cooled copper leads in Fig. 14. Similar behavior is encountered in superconducting devices where during several experi-

ments it was observed that the quenching to a normal state of current carrying high temperature superconducting wires, tapes or films was followed by the sample local destruction due to overheating, Vysotsky et al. [47], Romanovskii and Watanabe [48].

R.N. Krikkis / Cryogenics 83 (2017) 8–16

Maeda and Yanagisawa [49] underscored the inherent difficulties in the protection of HTS magnets from abrupt thermal runaways. Therefore thermal runaway is a common problem for the components of the superconducting magnet including the coil and both the HTS and the metallic parts of the current leads. Constructing a diagram qualitatively similar to Fig. 8 for example, where the curves for each component are displayed together against a common bifurcation measure, will show which component is more susceptible to blow-up, i.e. which limit point ‘‘comes first”. Apart from lowering the heat leak to the cryogenic liquid a variable cross section introduces a stabilizing effect since the blow-up threshold is substantially shifted to higher generation number values enabling higher operating currents as presented in Fig. 11.

325 300

lead temperature T (x) [K]

14

275

Measurements from Fang et al. (2013) RRR = 40, p = 1.0 F = 0.084, z = 1.0

250 225

G = 5.76 (I = 300 A)

200

G =10.24 (I = 400 A)

175 150 theoretical, u = 2.96 experimental

125

3.2. Vapor cooled leads, ideal convection

100 0.0

When the lead is cooled by a gas stream and the convective cooling is considered as ideal, that is Hg ! H, Eq. (11) takes the following form: 0

ðky H0 Þ þ 2

zGq ¼ F H0 : y2

ð23Þ 0

Hence, the introduction of the linear term F H is not expected to change qualitatively the multiplicity and stability structure. Indeed as it is shown in Fig. 12, the presence of the convective flow reduces the heat load, enables higher operating currents and extends the blow-up threshold. However, the multiplicity is still present and the structure is quite similar to the ‘dry’ case where F ¼ 0. The multiplicity is also observed when the cold end temperature is raised as depicted in Fig. 13 because the heat generation is the dominant factor controlling the temperature profile. Blow-up behavior is quite similar to the conduction cooled case since the heat generated by the current is significantly higher than the heat dissipated by conduction and convection combined, as shown in Fig. 14. 3.3. Vapor cooled leads, finite heat transfer coefficient This is the most general case since the variations of the gas temperature along the conductor are taken into consideration and the bifurcation analysis is carried out for the system of Eqs. (11) and (12) and an additional parameter, namely the conductionconvection parameter is introduced. The effect of the finite heat transfer coefficient is shown in compact form in Fig. 15 where the limit points have been projected of the ðu; GÞ plane for flow numbers in the range 0:1 6 F 6 0:8. As expected the limit points are bounded from one side by the common conduction limit as

60

0.8

p = 1, z = 1, RRR = 40

generation number G

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

Fig. 16. Theoretical calculations compared with measurements from Fang et al. [50].

u ! 0 and by the ideal convection limit as u ! 1 which is different for each value of the flow number, as shown in Fig. 12. The overshooting described by Chang et al. [8] is also observed here and it is more pronounced for low flow numbers. Since the limit points define the blow-up threshold, operating above the conduction limit line, a stoppage of flow fault will cause a thermal runaway unless protection devices are activated. 3.4. Comparison with experiment Fig. 16 shows a comparison between the theoretical calculations and a set of experimental data conducted by Fang et al. [50] for cylindrical copper leads cooled by nitrogen vapor. For this configuration a constant heat transfer coefficient was assumed H ¼ 38:6 W=ðm2 KÞ, the flow number is F ¼ 0:0836 corresponding _ ¼ 1:0 105 kg/s and the conductionto a flow rate of m convection parameter is u ¼ 2:96. The measured temperatures at the lead ends were used as boundary conditions. For this value of u no appreciable difference between the two convection models is expected according to Fig. 15 and the agreement between theory and experiment is reasonable. Apart from the reasons described by the authors in [50] for the deviations observed, it is believed that by adopting a position dependent Nusselt number (Graetz’s solution) instead of the asymptotic one with Nu ¼ 4:36 for laminar flow and a constant heat flux in a tube, the accuracy might have been even better. Since the measurement of the temperature distribution across the lead is difficult, a convenient indicator of the average temperature is the voltage difference between the warm and hot ends:

Z V ¼I 0

52

0.2

distance along conductor x

56 0.6

0.1

L

^ =AÞ dX: ðq

ð24Þ

Employing the reduced variables in Eq. (5) and introducing a refer48

^ K DTÞref , the dimensionless voltage ence voltage as V 2ref ¼ ðq during a current transient zðsÞ may be expressed as:

0.4 0.3

44

v ðsÞ ¼ ½zðsÞG1=2

0.2 40

F = 0.1 conduction limit

36 0.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

4.5

5.0

conduction-convection parameter u Fig. 15. Projection of the limit points on the ðu; GÞ plane. T L ¼ 4:2 K, T H ¼ 300 K.

Z

0

1

y2 qðx; sÞ dx

v ðsÞ ð25Þ

A qualitative comparison between the blow-up voltages for brass leads from Lee et al. [51] and copper leads, is depicted in Fig. 17. The comparison is restricted to the time interval close to the beginning of the runaway process, because the transient waveforms of the overcurrent used in the experiment and in the present simulation are certainly different. The initial conditions and the lead con-

R.N. Krikkis / Cryogenics 83 (2017) 8–16

time t [sec]

280 300 320 340 360 380 400 420 440 460 480 500 900 5.0 theoretical, present study measurements, Lee et al. (2001)

4.6

890 880

4.4 4.2

870

4.0

860

3.8

850

V [mV]

dimensionless voltage v (τ)

4.8

3.6 840 3.4 830

3.2 3.0 0.10

0.15

0.20

0.25

0.30

0.35

0.40

0.45

820 0.50

dimensionless time τ Fig. 17. Qualitative comparison of runaway voltages for copper and brass leads. Measurements adapted from Lee et al. [51].

figuration for the numerical computation is the same as in Fig. 14. Despite the differences in the material properties and the cooling method, the runaway features appear remarkably similar, at least during the early stages. This may be explained by the fact that the runaway phenomenon is dominated by the Joule heating and the electrical resistivity’s nonlinear dependence on the temperature rather than the thermal conductivity, the specific heat capacity or the cooling method. 4. Conclusions A numerical bifurcation analysis has been carried out for vapor and conduction cooled copper current leads powering a superconducting magnet at cryogenic temperature. The problem was based on one-dimensional axial conduction and two solutions have been calculated, one stable and one unstable, regardless of the cooling method. The limit points separating the stable from the unstable solutions form the temperature blow-up threshold since an applied current exceeding the limit point would result in a thermal runaway. Another worth mentioning finding is that the multiplicity is persisting even if the cold end temperature is raised to 100 K because of the dominant effect of the heat generation compared to the heat dissipation by convection and/or conduction. A possible application of the present model for the blow-up limits, to the complete superconducting magnet including the coil and composite leads with a HTS at the cold part is also discussed. A variable cross section appears to have a stabilization effect by extending the blow-up threshold allowing for higher operating currents. References [1] Wilson Martin N. Superconducting magnets. New York: Oxford University Press; 2002. [2] Iwasa Yukikazu. Case studies in superconducting magnets. Design and operational issues. 2nd ed. New York: Springer; 2009. [3] McFee R. Optimum input leads for cryogenic apparatus. Rev Sci Instrum 1959;30:98–102. [4] Buyanov Yu L, Fradkov AB, Yu Shebalin I. A review of current leads for cryogenic devices. Cryogenics 1975;15:193–200. [5] Katheder H, Schappals L. Design and test of a 10 kA gas-cooled current-lead for superconducting magnets. IEEE Trans Magn 1981;17:2071. [6] Maehata K, Ishibashi K, Wakuta Y. Design chart of gas-cooled current leads made of copper of different RRR values. Cryogenics 1994;34:935–40. [7] Lehtonen J, Korpela A, Mikkonen R, Tuisku J. Numerical optimization of conduction cooled current leads for low current applications. Cryogenics 2005;45:336–42. [8] Chang H-M, J Byun J, Jin H-B. Effect of convection heat transfer on the design of vapor-cooled current leads. Cryogenics 2006;46:324–32.

15

[9] Jones MC, Yeroshenko VM, Starostin A, Yaskin LA. Transient behaviour of helium-cooled current leads for superconducting power transmission. Cryogenics 1978;18:337–43. [10] Aharonian G, Hyman LG, Roberts L. Behaviour of power leads for superconducting magnets. Cryogenics 1981;21:145–51. [11] Buchs K, Hyman LG. Investigation of multiple solutions for gas cooled current leads and results of survey information on current leads. Cryogenics 1983;23:362–6. [12] Hanzelka P. Current leads in vacuum space of cryogenic systems. Cryogenics 1999;39:955–61. [13] Hull JR. High temperature superconducting leads for cryogenic apparatus. Cryogenics 1989;29:1116–23. [14] Seydel R. Practical bifurcation and stability analysis. 3rd ed. New York: Springer-Verlag; 2009. [15] Govaerts WJF. Numerical methods for bifurcations of dynamical equilibria. Philadelphia: SIAM; 2000. [16] Golubitsky M, Scaeffer D. Singularities and groups in bifurcation theory, vols. 1 and 2. New York: Springer; 1985. [17] Eckert D, Endig N, Lange F. Gas cooled current leads of variable cross section for helium cryostats. Cryogenics 1970;10:138–41. [18] Jiahui Z, Juntao T, Ming Q, Bin Y, Panpan C. Optimal design and experimental study on gas cooled current lead with varying cross-section. IEEE Trans Appl Supercond 2010;20:1705–9. [19] Okolotin VS, I Bol’shakov V. Gas cooled current leads with variable crosssection. Cryogenics 1972;12:227–9. [20] Chang H-M, Kim MJ. Optimization of conduction current leads with unsteady operating current. Cryogenics 2009;49:210–6. [21] Xu QJ, Zhu ZA, Liu LQ, Zong ZG, Huang SK, Yi CL, et al. Study on unsymmetrical phenomena of vapor cooled current leads. Cryogenics 2006;46: 676–82. [22] Krikkis RN, Razelos P. Optimum design of spacecraft radiators with longitudinal rectangular and triangular fins. ASME J Heat Transfer 2002;124:805–11. [23] Krikkis RN, Sotirchos SV, Razelos P. Analysis of multiplicity phenomena in longitudinal fins under multi-boiling conditions. ASME J Heat Transfer 2004;126:1–7. [24] Krikkis RN. On the multiple solutions of boiling fins with heat generation. Int J Heat Mass Transfer 2015;80:236–42. [25] Krikkis RN. Laminar conjugate forced convection over a flat plate. multiplicities and stability. Int J Therm Sci 2017;111:204–14. [26] Michelsen ML, Villadsen J. Solution of differential equation models by polynomial approximation. Englewood Cliffs, NJ: Prentice-Hall; 1978. [27] Hust JG, Lankford AB. Thermal conductivity of aluminium, copper, iron, and tungsten for temperatures from 1 K to the melting point. Boulder, CO: National Bureau of Standards) NBSIR 84-300; 1984. [28] Hairer E, Nørsett SP, Wanner G. Solving ordinary differential equations I. Nonstiff problems. New York: Springer-Verlag; 1987. [29] Ascher UM, Mattheij RMM, Russel RD. Numerical solution of boundary value problems for ordinary differential equations. 2nd ed. Philadelphia: SIAM; 1995. [30] Witmer G, Balakotaih V, Luss D. Finding singular points of two-point boundary value problems. J Comput Phys 1986;65:244–50. [31] Pruess S, Fulton CT. Mathematical software for Sturm-Liouville problems. ACM Trans Math Softw 1993;17:360–76. [32] Budd CJ, Huang W, Russell RD. Moving mesh methods for problems with blowup. SIAM J Sci Comput 1996;17:305–27. [33] Huang W, Ma J, Russell RD. A study of moving mesh PDE methods for numerical simulation of blowup in reaction diffusion equations. J Comput Phys 2008;227:6532–52. [34] Gurevich AVI, Mints RG. Localized waves in inhomogeneous media. Sov Phys Usp 1984;27:19–41. [35] Maddock BJ, James GB, Norris WT. Superconducting composites: heat transfer and steady state stabilization. Cryogenics 1969;9:261–73. [36] Gurevich AVI, Mints RG. Self-heating in normal metals and superconductors. Rev Mod Phys 1987;59:941–99. [37] Bandle C, Brunner H. Blowup in diffusion equations: a survey. J Comput Appl Math 1998;97:3–22. [38] Werweij AP. Thermal runaway of the 13 kA busbar joints in the LHC. IEEE Trans Appl Supercond 2010;20:2155–9. [39] Willering GP, Bottura L, Fessia P, Scheuerlein C, Werweij AP. Thermal runaways in LHC interconnections: experiments. IEEE Trans Appl Supercond 2011;21:1781–5. [40] Mumford FJ. Superconducting current-leads made from high T c superconductor and normal metal conductor. Cryogenics 1989;29:206. [41] Wesche R, Fuchs AM. Design of superconducting current leads. Cryogenics 1994;34:145–54. [42] Ballarino A. High temperature superconducting current leads for the Large Hadron Collider. IEEE Trans Appl Supercond 1999;9:523. [43] Lee H, Kim HM, Iwasa Y, Kim K, Arakawa P, Laughon G. Development of vaporcooled HTS-copper 6-kA current lead incorporating operation in the currentsharing mode. Cryogenics 2004;44:7–14. [44] Yang JY et al. Experimental study of the new type of HTS elements for current leads to be applied to the nuclear fusion devices. IEEE Trans Appl Supercond 2012;22:4801204. [45] Drotziger S et al. Investigation of HTS current leads under pulsed operation for JT-60SA. IEEE Trans Appl Supercond 2012;22:4801704.

16

R.N. Krikkis / Cryogenics 83 (2017) 8–16

[46] Heller R, Bauer P, Savoldi L, Zanino R, Zappatore A. Predictive 1-D thermalhydraulic analysis of the prototype HTS current leads for the ITER correction coils. Cryogenics 2016;80:325–32. [47] Vysotsky VS, Sytnikov VE, Rakhmanov AL, Ilyin Y. Analysis of stability and quench in HTS devices – new approaches. Fusion Eng Des 2006;81:2417–24. [48] Romanovskii VR, Watanabe K. Operating modes of high-T c composite superconductors and thermal runaway conditions under current charging”. Supercond Sci Technol 2006;19:541–50.

[49] Maeda H, Yanagisawa Y. Recent developments in high temperature superconducting magnet technology (review). IEEE Trans Appl Supercond 2014;24:4602412. [50] Fang J, Yu T, Li ZM, Wei B, Qiu M, Zhang HJ. Comparison and analysis of the efficiency of heat exchange of copper rod and copper wires current lead. Physica C 2013;494:364–9. [51] Lee H, Arakawa P, Efferson KR, Iwasa Y. Helium vapor – cooled brass current leads: experimental and analytical results. Cryogenics 2001;41:485–9.