Florida State University Libraries Electronic Theses, Treatises and Dissertations

The Graduate School

2004

Numerical Simulation of Quench Propagation in Superconducting Magnets by Using High Order Methods Shaolin Mao

Follow this and additional works at the FSU Digital Library. For more information, please contact [email protected]

THE FLORIDA STATE UNIVERSITY COLLEGE OF ENGINEERING

NUMERICAL SIMULATION OF QUENCH PROPAGATION IN SUPERCONDUCTING MAGNETS BY USING HIGH ORDER METHODS

By

SHAOLIN MAO

A Dissertation submitted to the Department of Mechanical Engineering in partial fulfillment of the requirements for the degree of Doctor of Philosophy

Degree Awarded: Fall Semester, 2004

The members of the Committee approve the dissertation of Shaolin Mao defended on October 18, 2004.

Cesar A. Luongo Professor Directing Dissertation

David A. Kopriva Outside Committee Member

Chiang Shih Committee Member

David A. Cartes Committee Member

Approved:

Chiang Shih, Chair, Department of Mechanical Engineering

Ching-Jen Chen, Dean, College of Engineering

The Office of Graduate Studies has verified and approved the above named committee members.

ii

To my parents with love

iii

ACKNOWLEDGEMENTS

I would like to thank the members of the Dissertation Committee at Florida State University for many useful suggestions during my writing of this work. I also would like to thank the staff at the Center for Advanced Power Systems, Florida State University, and Magnet Science and Technology Division, National High Magnetic Field Laboratory, for many valuable discussions during my research on superconducting magnets: Dr. John R. Miller, Dr. Andrey V. Gavrillin, Dr. Soren O. Prestemon (now at LBNL), Dr. Steven W. Van Sciver, and Dr. Stephen L. Woodruff. I want to thank Dr. Luca Bottura at CERN for his help and advice in using GANDALF code. I want to thank Dr. Chris K.W. Tam for his great courses. I would also like to thank Dr. David A. Kopriva for his instruction in my research work, his kind help saved me a lot time during my work. I especially thank my dissertation advisor Dr. Cesar A. Luongo for his enormous help, thoughtful support, and unbiased instruction on my writing. I also want to thank Prof. Yue Sheng Li at Zhong-Shan University, who lead me to the door of scientific computing. Finally, I would like to thank my family for their patience and support during this long-time endeavor.

iv

TABLE OF CONTENTS

List of Tables ............................................................................................................................vii List of Figures ..........................................................................................................................viii Abstract .....................................................................................................................................xii

1. INTRODUCTION ..................................................................................................................1

2. THE MECHANISMS OF QUENCH PROPAGATION ......................................................10 2.1 Physical description of quench propagation...................................................................11 2.2 Modeling of quench propagation ...................................................................................15 2.3 Boundary conditions ......................................................................................................25

3. THE ENERGY BALANCE MODEL AND NUMERICAL SOLUTION ...........................27 3.1 Thermal analysis of subcooled superfluid helium in CICC ...........................................27 3.2 Numerical algorithm ......................................................................................................33 3.3 Model results and discussions ........................................................................................35 3.4 Conclusions ....................................................................................................................51

4. NUMERICAL SIMULATION OF QUENCHING PROPAGATION AT EARLY PHASE BY HIGH ORDER SCHEMES ............................................................................52 4.1 Mathematical Model ......................................................................................................53 4.2 High order finite difference schemes .............................................................................55 4.3 Adaptive mesh technique ...............................................................................................58 4.4 Numerical stability analysis ...........................................................................................60 4.5 Numerical damping and testing......................................................................................66 4.6 Benchmarking ................................................................................................................71

v

4.7 Conclusions ....................................................................................................................79

5. DISCONTINUOUS GALERKIN SPECTRAL ELEMENT SIMULATION OF QUENCHING PROPAGATION IN CICCs ......................................................................80 5.1 Introduction....................................................................................................................81 5.2 The Discontinuous Galerkin application for Real Fluid (helium)..................................83 5.3 Implementation of discontinuous Galerkin spectral element methods ..........................87 5.4 Numerical results............................................................................................................92 5.5 Discussion ....................................................................................................................113

6. CONCLUSIONS.................................................................................................................115 6.1 Suggestions for future work .........................................................................................118

APPENDIX A. ADAPTIVE MESH TECHNIQUES .............................................................119

APPENDIX B. NOMENCLATURE ......................................................................................125

REFERENCES .......................................................................................................................127

BIOGRAPHICAL SKETCH ..................................................................................................133

vi

LIST OF TABLES

Table 3.1: Record of high-current operations, superconducting outsert ..................................36 Table 4.1: Nominal operating parameters in Arp’s case...........................................................72 Table 4.2: Nominal operating parameters in Ando’s case ........................................................75

vii

LIST OF FIGURES

Figure 2.1: Cross section views of CICC used for the three subcoils at NHMFL 45-T hybrid magnets ....................................................................................................12 Figure 2.2: Schematic of quench propagation in a CICC connected to a big helium cryostat ................................................................................................................13 Figure 2.3: A simple description of hot-spot calculation model .............................................16 Figure 2.4: A typical one-dimensional model of quench propagation ....................................18 Figure 2.5: Schematic of quasi one-dimensional CICC ..........................................................20 Figure 2.6: The cross section of three-dimensional model .....................................................24 Figure 3.1: Schematic of the quasi 1D energy balance model ................................................29 Figure 3.2: Schwmatic of the NHMFL 45-T superconducting outsert. There are 6 layers in the coil “A”, layer 1 is the highest magnetic field .............................................29 Figure 3.3: Evolution of total AC losses (hysteresis plus coupling) in layer 1/coil A during a 5 A/s current ramp ...................................................................................38 Figure 3.4: Evolution of index heating in layer 1/coil A. Ramp rate is 5 A/s, index number n=5 ............................................................................................................38 Figure 3.5: Index heating for layer 1/coil A during and following a 5 A/s current ramp (layer is adiabatic and insert is off)........................................................................39 Figure 3.6: Peak temperature in layer 1/coil A during a 5 A/s current ramp (layer is adiabatic and insert is off) ....................................................................................39 Figure 3.7: Temperature evolution of layer 1/coil A during a 5 A/s ramp ..............................40 Figure 3.8: Temperature evolution of layer 1/coil A during a 2.5 A/s ramp ...........................40 Figure 3.9: Temperature distribution in the six layers of coil “A” at t=2033s (end of a 5 A/s current ramp)..................................................................................................41 Figure 3.10: Peak temperature in layer 1/coil A for a different current ramp with the same index number n=15...............................................................................................41 Figure 3.11: Temperature, operating current density, and critical current density at the center of layer 1/coil A during a 5 A/s current ramp to 10 kA (layer is adiabatic, insert coil is off) ...................................................................................42 Figure 3.12: Index heating for layer 1/coil A during a 5 A/s current ramp with heat transfer between layers includes, insert coil is off ...............................................44 Figure 3.13: Comparison of index heating in layer 1/coil A for adiabatic and nonadiabatic condition (5 A/s current ramp, insert is off)..........................................45 Figure 3.14: Temperature evolution of layer 1/coil A during 5 A/s ramp and layer is adiabatic, insert coil is off ....................................................................................45 Figure 3.15: Temperature distribution in different layers with a 5 A/s ramp at time t=500 s (non-adiabatic condition, insert is off) ...............................................................46

viii

Figure 3.16: Magnetic field distribution along layer 1 and 6 of coil A at the end of 5 A/s current ramp, t=2000 s, rated I=10 kA.................................................................47 Figure 3.17: Peak temperature in layer 1/coil A during and after a 5 A/s current ramp for different index number (n=5, 15) .........................................................................48 Figure 3.18: Influence of index number on temperature evolution at non-adiabatic condition and insert coil is turned on ...................................................................48 Figure 3.19: Initial temperature distribution of layer 1/coil A used to calculate the recovery time ........................................................................................................49 Figure 3.20: Thermal recovery time of layer 1/coil A under adiabatic condition in the superfluid helium regime......................................................................................50 Figure 4.1: Schematic of quasi 1D CICC, three components are connect through heat transfer and friction force .....................................................................................54 Figure 4.2: Choices of special stencils in high-order finite difference schemes near the interfaces of two different mesh-size sub domains ..............................................60 Figure 4.3: Real part and imaginary part of the two roots in the complex plane of ω ∆t ( ω ∆t is a real variable) .................................................................................62 Figure 4.4: Real part and imaginary part of the two roots in the complex plane of ω ∆t ( ω ∆t is a pure imaginary variable) ...............................................................63 Figure 4.5: Stability region of the second-order and fourth-order time integration schemes in the complex pane ................................................................................63 Figure 4.6: Numerical dispersion errors as indicated by the dependence of the numerical wave number α on the α of PDE........................................................................65 Figure 4.7: Comparison of numerical solution by sixth order schemes at two time t=40, and 80.....................................................................................................................68 Figure 4.8: Comparison of numerical solution with the exact solution for a simple wave equation with a step input. The wave propagates 300 time periods ......................69 Figure 4.9: Comparison of numerical solution of discontinuity propagation when numerical damping is too large to smear the solution...........................................70 Figure 4.10: Dispersion errors for some finite difference schemes to wave equation propagation...........................................................................................................70 Figure 4.11: Evolution of helium quenching pressure in the Arp’s case ..................................73 Figure 4.12: Evolution of conductor temperature in the middle position in the Arp’s case, the conduct was recovered due to the helium.......................................................74 Figure 4.13: Normal zone evolution in a short sample of NbTi with a 26 m in length and a diameter of 5.97 mm ..........................................................................................76 Figure 4.14: Evolution of quenching pressure at the central position of the NbTi CICC, operating current I=1.5 kA...................................................................................76 Figure 4.15: Evolution of conductor temperature at the central position of the NbTi CICC, operating current I=1.5 kA........................................................................77 Figure 4.16: Evolution and distribution of induced- flow velocity of helium inside the NbTi CICC ...........................................................................................................77

ix

Figure 4.17: Evolution and distribution of conductor temperature inside the NbTi CICC ......78 Figure 4.18: Evolution and distribution of helium quenching pressure inside the NbTi CICC.....................................................................................................................78 Figure 5.1: Second-order convergence of the finite difference methods ................................81 Figure 5.2: A typical spectral convergence of Fourier collocation methods ..........................82 Figure 5.3: Interpolation errors by using second-degree polynomials to fit curve of heat capacity of supercritical helium............................................................................86 Figure 5.4: Stability region of Runge-Kutta methods of order one to four, the inner most is the first order methods, and outer most is the fourth-order methods................91 Figure 5.5a: Evolution and distribution of quenching pressure in a short CICC sample when the external heat pulse is 8× 103 W/m using DG-SEM ..............................93 Figure 5.5b: Evolution and distribution of quenching pressure in a short CICC sample when the external heat pulse is 8× 103 W/m using GANDALF code..................93 Figure 5.6a: Conductor temperature distribution in a short sample using DG-SEM................94 Figure 5.6b: Conductor temperature distribution in a short sample using GANDALF code ...94 Figure 5.7a: Distribution of induced-flow velocity of helium in a short sample using DGSEM with an external pulse of 8× 103 W/m ........................................................95 Figure 5.7b: Distribution of induced-flow velocity of helium in a short sample using GANDALF code with an external pulse of 8× 103 W/m.....................................95 Figure 5.8: Distribution and evolution of numerical solution in a short sample by GANDLAF code ..................................................................................................96 Figure 5.9: Distribution and evolution of quenching pressure in the same sample when a large external pulse of 4 ×10 4 W/m is set...........................................................97 Figure 5.10: Distribution and evolution of conductor temperature in the same sample with a large external pulse of 4 ×10 4 W/m.........................................................97 Figure 5.11: Convergence comparison by using DG-spectral element methods for different time steps and mesh size, I=1.5 kA ......................................................98 Figure 5.12: Convergence comparison by using DG-spectral element methods for different time steps and mesh size, I=1.8 kA ......................................................99 Figure 5.13: Convergence test of conductor temperature by using DG-spectral element methods for different time step and mesh size, I=1.5 kA ..................................99 Figure 5.14a: Distribution and evolution of helium density by using DG-spectral element methods, I=1.5 kA............................................................................................100 Figure 5.14b: Distribution and evolution of helium density by using GANDALF code, I=1.5 kA ...........................................................................................................100 Figure 5.15a: Distribution and evolution of quenching pressure by DG-spectral element methods, I=1.5 kA............................................................................................101 Figure 5.15b: Distribution and evolution of quenching pressure by using GAN DALF code, I=1.5 kA..................................................................................................101

x

Figure 5.16a: Distribution and evolution of induced- flow velocity of helium by DGspectral element methods, I=1.5 kA..................................................................102 Figure 5.16b: Distribution and evolution of helium induced-flow velocity by using GANDALF code, I=1.5 kA...............................................................................102 Figure 5.17a: Distribution and evolution of conductor temperature by DG-spectral element methods, I=1.5 kA ..............................................................................103 Figure 5.17b: Distribution and evolution of conductor temperature by using GANDALF code, I=1.5 kA ...................................................................................................103 Figure 5.18: Validation of the physical model with experimental data and the numerical solution for I=1.5 kA in Ando’s case .................................................................105 Figure 5.19: Validation of the physical model with experimental data and the numerical solution for I=1.8 kA in Ando’s case ................................................................105 Figure 5.20: Evolution of conductor temperature by GANDALF code and DG-spectral element methods for I=1.5 kA in Ando’s case ...................................................106 Figure 5.21: Evolution of conductor temperature by GANDALF code and DG-spectral element methods for I=1.8 kA in Ando’s case ...................................................107 Figure 5.22: Evolution of quenching pressure by GANDALF code and DG-spectral element methods for I=1.5 kA in Ando’s case ...................................................107 Figure 5.23: Evolution of quenching pressure by GANDALF code and DG-spectral element methods for I=1.8 kA in Ando’s case ...................................................108 Figure 5.24: Evolution of quenching pressure by GANDALF code and DG-spectral element methods for I=2.0 kA in Ando’s case ...................................................109 Figure 5.25: Thermal- hydraulic quench back (THQB) phenomenon is observed if the friction factor is large enough.............................................................................110 Figure 5.26: Effect of convection of helium will determine the numerical results of quenching pressure .............................................................................................111 Figure 5.27: Evolution of maximum quenching pressure for I=2.0 kA in Ando’s case .........113 Figure A.1: Time integration methods for adaptive mesh......................................................119 Figure A.2: Inverse interpolation of I(x) is used to generate an equi-distribute mesh ...........121 Figure A.3: A typical moving grid trajectories for the algorithm used with nondimensional scale. µ1 = 10, µ 2 = 1000 .................................................................122 Figure A.4: A typical moving grid trajectories for the algorithm used with nondimensional scale. µ1 = 50, µ 2 = 1000 ................................................................123

xi

ABSTRACT

Large-scale superconducting magnet systems operate in a high current environment. System protection during quench (loss of superconductivity) is always a key consideration. Thermalhydraulic analysis of superconducting magnets is used to provide criteria for the design and operation of such facilities, and research on quench propagation is a main topic in thermalhydraulic analysis. In this study, numerical simulation is applied to compressible helium flow and coupled heat transfer in cable-in-conduit superconductor (CICC). Helium flow in CICCs and quench propagation simulation are relevant engineering problems. These flows have high Reynolds number and low Mach number. The large non- linear source terms due to Joule heating and friction force introduce difficulties for most numerical techniques. The main consideration in this dissertation is to seek some numerical methods with high accuracy (resolution) and efficiency. In the first part of this study, a simple physical model, the energy balance model, was proposed for the first time to track the superfluid helium (He II) and normal helium (He I) fronts, which is very critical to simulate He II/He I transition in large-scale superconducting magnets. This new model was used to ana lyze the thermal stability (quench behavior) of the NHMFL 45-T hybrid magnets system. This model resulted in high efficiency of numerical simulation of thermal stability analysis compared to complicated 1D quench propagation model. Besides, adaptive mesh techniques were introduced to improve numerical efficiency. This model was successfully predicted quench behavior, right protection, and showed that index heating was the cause to drive quench in the NHMFL 45-T hybrid magnets system. In the second part of this study, two high-order methods were applied to simulate quench propagation in CICC magnets. The main consideration in this dissertation was to seek some numerical methods with high accuracy (resolution) and efficiency. The first high-order finite different methods, dispersion-relation-preserving (DRP) schemes were used to simulate quench

xii

propagation in CICC at the early phase in order to decrease the phase errors (dispersion errors). High-order methods were introduced for the first time to simulate quenc h propagation in superconducting magnets. This research was focusing on the early phase of quench propagation, which helps to understand the quench propagation mechanism. Detailed numerical algorithm was discussed, and numerical benchmarking results were obtained. We conclude that the research on the early phase of quench propagation is an important issue to understand quench propagation in detail. The second high-order methods, discontinuous Galerkin (DG) spectral element methods (SEM) were applied to overcome numerical difficulties encountered by most classical methods. DG-SEM is very robust to simulate quench propagation for which traditional methods were only highly efficient either at early stage or later stage of quench propagation in CICC. Detailed numerical benchmarks showed that DG-spectral element methods have advantages in dealing with large non- linear source term problems (quench propagation) by avoiding artificial damping and keeping the system numerical conservation. DG-spectral element methods showed high accuracy (high resolution in large gradient regions), they also demonstrated high numerical stability for large external disturbances. DG-spectral element methods can be parallelized to speed up the simulation while retaining high accuracy. Roe’s approximate Riemann solver for real fluid/gas (liquid helium) was processed successfully by using polynomial interpolation to fit heat capacity of helium in Riemann integral for the first time, which is a key step to apply DGspectral element methods in quench propagation simulation. From this study we conclude that high-order finite difference schemes and DG-spectral element methods obtain better results in accuracy and efficiency than lower-order methods. Finally, a Fortran code was developed based on the study of high-order methods in this dissertation.

xiii

CHAPTER ONE INTRODUCTION

Superconducting magnets have many applications in research and industry, especially when high magnetic fields are needed or economical or technical reasons limit electrical power availability. Because of the high cost associated with the design and construction of large-scale magnets, system protection (under abnormal operating conditions) is always a key consideration. Most large-scale magnets are cooled with supercritical helium, typically at a temperature near 4 K and pressure below 1 MPa, although sometimes superfluid helium is used. The most striking properties of superfluid helium are its complete lack of viscosity and its very high thermal conductivity. It is the high thermal conductivity that makes superfluid helium attractive as a coolant. Heat conduction in stagnant superfluid helium obeys the Gorter-Mellink law instead of classic Fourier’s law [1]. It is worthwhile to understand how large-scale superconducting magnets are built. A very common configuration in large-scale magnets is to use cable- in-conduit conductors (CICC), which are the preferred choice for magnets when operating conditions require a reliable and stiff design [1]. The conductors of such magnets consist of superconducting composite wires braided into cables and placed in a strong jacket, which contains helium fluid. In some of these conductors, the cable is held tightly by the jacket which is usually swaged around the conductor, and the supercritical helium circulates through the interstices of the cable. This type of conductor is called a cable- in-conduit conductor, the strands of the cable are unrestrained except where they cross each other, wire motion is possible in such conductors. In other types of internally cooled conductors, the cable is soldered together and the helium fluid circulates through a tube down its center [2]. CICCs are insulated with glass fiber braid. Vacuum impregnation with epoxy resin has been the most popular technique so far. Pure copper is usually added to superconducting composites in order to promote dynamic stability and protection from quench [3].

1

For the design of most large-scale superconducting magnets, a stability analysis is necessary. The aim of the stability analysis is to see if the magnet will work, and if it is stable. We are especially interested in the transition when the superconducting magnet goes to the normal state (resistive state). Because of abnormal operating conditions such as external mechanical energy (friction, etc.) or any external heat deposition, the superconductor goes above the stability margin very quickly, usually in a time scale of one or two seconds. This transition to thermal runaway is called “quench” propagation. A great disappointment in the early days of magnet construction was the discovery that the performance of superconducting wires, when wound into magnets, always fell far short of what might have been expected from the results on short samples of conductors. This effect came to be known as ‘coil degradation’, which is due to quenched coils. Quench is a term commonly used to describe the process which occurs when any part of a magnet goes from the superconducting to the normal resistive state. Because normal-state resistivities and current densities are both high, intense local heating will ensue, taking the quench point and surrounding region to a temperature far above critical. This is an irreversible process causing the entire stored energy of the magnet to be dissipated as heat. One can only turn off the current supply, wait for the magnet to cool down and then try again [3]. Many early researchers have put lots of effort on finding methods to avoid quench. Some strategies to stabilize magnets are: training, to add more copper to the conductor, or to reserve some void space during the winding of the coils for helium to percolate. When enough cooling is provided, the quenching magnet can be cooled back down to the superconducting state, which is called recovery. Magnets stabilized in this way are called cryostable. To the contrary, some magnets with less copper in the conductor are operated with higher current densities than cryostable magnets. This type of magnet is metastable [2]. In the 1980s, superconductor stability analysis became a more interesting field of research. Stability interests us because it is indispensable for continuity of operation. Truly stable magnets were built by immersing the magnets in a pool of boiling helium and reducing the overall current density by adding copper until the Joule power per unit conductor surface in the normal state was less than the minimum film boiling heat flux. This kind of stability became known as cryostability because it was achieved by removing the Joule heat with cooling. Cryostability is attractive because it relieved us of the burden of knowing exactly what thermal perturbations the

2

magnet will suffer. Some workers, endeavoring to maintain cryostability, have tried to improve heat transfer between the conductor and helium. These efforts have taken several different roads, such as improvement of boiling heat transfer by treatment of the conductor surface, renunciation of pool boiling in favor of internally cooling the superconductors, and use of superfluid helium as coolant [2]. The stability of CICC, and more generally any interna lly cooled conductor, is strongly influenced by the motion induced in the helium during recovery. Once an initial normal zone has been created, most of the Joule heat it produces is transferred to the helium adjacent to it. The helium, being confined laterally by the conduit, expands longitudinally because of the Joule heating. The velocity it attains as a result of such heating can be quite high, with values of tens of meters per second or more being possible. Such high velocities of flow over the surface of the conductor substantially augment the transfer of heat from the conductor to the helium and thereby strongly influence the stability of the conductor [4]. The effect of the induced flow on the stability of CICC is a complicated story involving a sometimes multi- valued stability margin. The stability margin of a CICC is defined as the largest sudden, uniform heat deposition in the strands of the cable from which recovery of the superconducting state is still possible. It is limited by the available heat capacity of the helium between the bath temperature and the current sharing threshold temperature. This is because recovery from a thermal excursion takes only tens of milliseconds, whereas the residence time of helium in the coil is typically tens of minutes. So for the purpose of deciding the issue of stability, the helium inventory can be considered limited. Stability margin has been measured by introducing a heat pulse to the conductor and observing whether the conductor subsequently recovered or que nched. In some experiments, the double sequence recovery-quench-recovery-quench was observed. In other words, sometimes the conductor recovered from small pulses, quenched for large pulses, recovered for yet larger pulse and again quenched for largest puls es. This type of multiple stability has been repeatedly observed in the laboratory [2]. The most interesting parameters of thermal hydraulic analysis of superconducting quench propagation are the normal zone length, and normal zone propagating speed evaluated through experimental or numerical approach. Some theoretical analysis, such as equal-area theorem, minimum propagating zone (MPZ), can be found in [3], [5].

3

Analytic expressions were proposed in the 1980s based on the similarity solution approach. The aim is to approximate the thermal hydraulic process during quench using some simple onedimensional models. Dresner has developed an approximate analytic expression for the maximum pressure during quench [6]. The similarity solution has been shown to be in good agreement with experiment. The helium expulsion velocity at the ends of a complete normalizing cable was also proposed by Dresner [6], again under some approximating assumptions such as friction force and constant properties of helium. Other theoretical work has been done by Shajii and Freidberg [7], giving an approximate expression for quench propagation speed and pressure distribution, neglecting the inertia term and assuming the same temperature distribution in conductor and helium. A very important facet of quench simulation is to calculate the normal zone propagation correctly. Some early results were given by Dresner [8], [9]. Because the main quench propagation mechanism in CICC is a hot helium expulsion, Dresner suggested that the velocity of the normal zone propagation equals the local velocity of expansion of the helium. After introducing a dimensionless similarity, Shajii and Freidberg developed asymptotic expressions for the normal zone propagation speed and the pressure increase in the normal zone [7]. Before expanding on the subject of the numerical simulation of quench propagation problem, it is necessary to review the mathematical models developed in the past two decades. This problem is so complicated that it is impossible to find a practical approach to completely solve it in three-dimensions, even with the most powerful computing facilities. Hence, most numerical work has been based on one-dimensional or quasi two-dimensional models by introducing some reasonable simplifications to the physical model of quench propagation. Arp gave a computational

illustration

of

one-dimensional time-dependent

quench

in

forced-flow

superconductors. He also described the physical boundary conditions through a intuitive analogue [10]. Marinucci took into account transients in the heat transfer coefficient for analysis of stability and quench characteristics. Finite difference in space and fourth order Runge-Kutta time integration were used to solve the system [11]. There are some important works at the end of the 1980s and the beginning of the 1990s. The early numerical work done by Wong is based on a one and a half dimensional model [12]. A finite volume method and staggered grid technique were used. However, the calculation is computationally expensive for complete a three-dimensional model. Luongo et al. developed a

4

1D model by assuming the same temperature between helium and metal. Semi-discretization was used to obtain a system of ODE or DAE (differential algebraic equations) and then a method of line (MOL) was applied to solve these equations (time integration) [13]. Bottura gave a complete description of the 3D model of the quench problem in CICC [14]. Based on some assumptions he reduced it to a 1D model. FEM was applied to solve the system [15], [16], while using the pressure drop equation instead of the momentum equation [17]. Shajii and Freidberg started by simplifying the general one-dimensional model to obtain a uniform compact system. A spline collocation method was used to solve the simplified governing equations [18]. Cheng analyzed thermally induced transient flow of single-phase supercritical helium in a forced- flow cooling system, applying an implicit scheme to the convection term and an explicit scheme to the diffusion Fourier heat conduction term [19]. Ito discussed a two-dimensional model of coolant flow of supercritical helium [20]. Koizumi et al. investigated transient heat transfer from strands to helium and increased the accuracy of the numerical simulation. The governing equations were discretized by the Beam-Warming scheme where implicit time integration was used to avoid the Courant-Friedrichs-Lewy (CFL) condition restriction [21]. He also improved his work further by solving the thermal boundary layer equations [21], [22]. Wang et al. have introduced an adaptive mesh method to solve the quench governing equations with good efficiency and accuracy [23], [24]. The review above encompasses most of the work on quench simulation of the past 20 years, but by no means is an exhaustive listing. A large-scale superconducting magnet is a very complex system. It is unrealistic to directly measure the normal zone propagating speed and the pressure rise inside the CICC [25], [26]. The transition from the superconducting to the resistive state in a CICC is usually fast and takes place in a few seconds, yet knowledge of this behavior is crucial for the design of a protection scheme. The fluid flow in a CICC cooled by supercritical helium is in the turbulent regime, in which the Reynolds number is up to 105 or higher. Unique properties of helium fluid, heat capacity discontinuity of normal helium and superfluid helium, strong thermal coupling effects between the fluid and the strands/conduit wall all prove to be a mathematically stiff problem. Large temperature gradients in the quench front lead to sharp variations in the thermal properties of the solid and liquid, which may increase by two to three orders of magnitude over a short distance, also posing a computational problem. The re exists a moving front during the quench propagation which is a special boundary, different from shock or other discontinuities [27], [28]. Since the

5

moving front (edge of the normal zone) is a few centimeters in length, compared to 100-1000 meters of the CICC, high computational resolution is difficult even in 1D, and nearly impossible in 2D or 3D. There are also two time scales: the quench front velocity and the speed of sound in liquid helium, again, that create a computational difficulty. All these factors make the analytic solution of the quench problem extremely difficult, if not virtually impossible, leading to the need in practice to develop numerical simulation techniques. Quench propagation in CICC magnets is a special moving front problem, a coupled threeelement system problem. If superfluid helium is used in CICC magnets, He II/He I transition front (phase change) is another moving boundary problem before quench propagation happens in CICC. Therefore, it is very important to choose numerical methods to improve the resolution of the large gradient region of solutions. A detailed investigation of numerical approaches of phase change problems can be found in Crank [29]. There are many efficient methods to solve phase change problems with moving boundaries, such as front-tracking method [29], fixed grid method [30], phase field method [31], and level set method [32]. Front-tracking method needs to calculate the phase change fronts explicitly at each time step. Interpolation is usually used near the fronts. The implementation is straightforward, however, this method may sometimes be difficult or even impossible to track a complicated moving front directly. Enthalpy method is a kind of fixed grid method, which avoids calculating the moving front explicitly, and the position appears, a posteriori, being a feature of this method. The major advantage of fixed grid methods is its versatility to three-dimensional problems. Phase field methods also avoid the need to explicitly track the moving interface. In place of introducing the enthalpy function the phase field function is supplied based on Landau-Ginzburg theory. The level set method has gained popularity in recent years, especially tracking complicated fronts, such as dendritical configurations. In the current study front-tracking methods developed for conventional solid/liquid or liquid/gas phase change transformations are applied to solve He II/He I transition moving boundary problem. Also moving mesh methods are introduced to improve the quality of the solution. In the first part of this study, we propose a compact physical model for thermal stability analysis in superconducting magnets (long time quench initiation in CICC). A goal of this dissertation is to seek an effective approach to simulate quenc h propagation. In quench propagation simulation, computational time (algorithm efficiency) is always an important

6

consideration, especially for long time simulation. It usually needs more than three weeks to use the UNIX workstation (SunBlade 2000) to simulate three hours of evolution of helium flow in CICC by using usual quasi-1D quench models. Therefore, we need to seek a compact model effectively to simulate this long time transition phenomenon (He II/He I) in CICC magnets, which may involve hours of real time thermal transients. A compact physical model will be proposed for the first time based on experimental and operational observation. This simple energy balance model is used to track the He II/He I fronts in superconducting magnets which use He II (static usually) as the coolant. As an example, we use this model to analyze the NHMFL 45-Tesla hybrid magnet system at the National High Magnetic Field Laboratory. The numerical results by this model will be compared with operational observations and experimental data. To conclude this part of my study, a new physical model will be introduced in Chapter 3, and numerical algorithms are discussed to solve this model, and application to the NHMFL 45Tesla hybrid magnets system will be discussed. In the second part of this study, we will focus on a quasi-1D quench propagation model and with an interest in the unsteady problem in a short quench period (often several seconds). Therefore, computational accuracy is another important consideration in quench propagation simulation. As a matter of fact, there is a tradeoff between computation accuracy and efficiency, or the expense of computation. As we reviewed before, there are many choices to solve this problem (quench propagation in a short time). Finite difference method can be easily implemented [33], finite element method has its merit when applied to domains of complicated geometry [34], spectral method contains exponential accuracy for smooth solution [35], finite volume method inherits the conservation of flux to governing equations by introducing staggered grids [12], [36]. Even though there are some commercial codes for quench propagation simulation, for example, the finite element GANDALF code, there is still a need to improve the computational efficiency for long time simulation and to improve computational accuracy (resolution in large gradient region of solution) and more numerical stability to large external disturbances. Another goal of this dissertation is to seek new effective solver for the compressib le subsonic Navier-Stokes equation (helium flow) and coupled heat equation (solid materials) for accurate discretization in time and space. To identify a preferred approach to solve the problem of interest (quench propagation in CICC for a short period) is not a trivial matter. Quench propagation

7

simulation is too complicated to find an absolutely advantageous solver. Different solvers were designed specially for different aspects of the application. Most numerical methods discussed above are low-order sche mes, either in spatial or temporal discritization. As suggested by the name, low-order methods need smaller time step or finer grid to keep the same computational accuracy compared with high-order methods. For wave propagation problems, low-order methods will introduce much more dispersion errors (phase errors in wave propagation) than high-order methods. In this study, two new high-order methods, dispersion-relation-preserving scheme (DRP) in chapter 4 and the discontinuous Galerkin spectral element methods in Chapter 5, will be implemented to solve the quench propagation problems. Fourth-order spatial and temporal discretization are used when implementing DRP method [37]. The DRP scheme has the advantage over traditional finite difference schemes of decreasing dispersion errors (phase errors in wave propagation). Considering the early phase of quench propagation in CICC, the diffusion is very small (friction force in the momentum equations and the coupling heat transfer in the energy equation are both small), and the quench propagation is more like the acoustic wave problem. We can keep very small dispersion errors and diffusion errors by using DRP schemes for quench propagation at the early stage. Discontinuous Galerkin (DG) spectral element methods (SEM) are also considered in this study. Spectral element method [38] can be regarded as a high-order finite element method (which is similar to p-type finite element method) that combines the geometrical flexibility of finite element methods (h-type) with the rapid convergence of spectral techniques [39], [40], [41]. DG-SEM has the advantage in processing the boundary conditions, imposition of the numerical flux across the interfaces between adjacent elements in the computational domain. DG-SEM also overcomes the difficulties encountered in central difference methods (artificial damping term needed) by applying some approximate Riemann Solver. Besides, DG-SEM makes the imposition of boundary condition easier and without degeneration of order of accuracy near the boundary of the computational domain. In the application of DG-SEM to quench propagation in CICC magnets, the most difficult part is the implementation of Riemann solver for a real fluid/gas (compressible helium flow). In this study, a new Roe’s approximate Riemann solver for liquid helium will be proposed and discussed in Chapter 5. This work is the basis of the implementation of discontinuous Galerkin spectral element methods to quench

8

propagation simulation in CICC magnets. Detailed numerical results will also be shown in this study to benchmark the new methods in quench propagation. DG-SEM demonstrates high resolution of solution, and high numerical stability to large external disturbances. In the end of this study, our computational results will be benchmarked against with an established commercial code, such as GANDALF, and with experimental data. In short, we will study quench propagation computational simulation in two aspects, one is to find a simple physical model effectively for long time evolution simulation, another is to seek high-order methods for a short time quench propagation simulation with high accuracy and efficiency at the same time.

9

CHAPTER TWO THE MECHANISMS OF QUENCH PROPAGATION

The quench problem is of interest both as a complicated “unsolved” physical problem and as a good vehicle for describing and illustrating the efficiency and versatility of our numerical techniques. Most importantly it is of a first principle understanding for the design robust magnet protection schemes. It is not possible to develop efficient, accurate, and stable numerical techniques for quench simulation without due consideration of the underlying physics and mathematics. Because most large-scale magnets are built with CICCs, we will concentrate our attention on CICC, where numerical simulation is about internal thermal-hydraulics. If a conductor is in a superconducting state, there will be no losses (ideally) in the magnet, and the helium mass flow or static supercritical helium will remove all the heat originated from non-perfect thermal shielding of the outside environment, and other heat sources, such as conduction [42]. At the design condition, the cooling is enough to keep the superconducting coil in stable operation. However, due to unforeseen disturbance, such as the increase of steady-state heat input, the slipping among components under the action of the Lorentz force, cracking of the epoxy, de-bonding of layer-to-layer insulation, or index heating [43], [44], additional heat could be deposit ed in the conductor. The range of thermal perturbations occurring in superconducting magnets is very poorly known. In the worst case scenario, when the conductor is driven to the normal (resistive) regime it will generate Joule heat in addition to the external heat sources above. A positive feedback effect will accelerate the transition to the resistive normal state in which the ensuing resistive dissipation will depend on the percentage of copper included in the conductor and the copper resistivity: a bigger cross-section area of copper or a lower resistivity will result in a smaller resistance and a smaller heat dissipation [14]. If the heat balance between heat generation and removal breaks up, the conductor will either recover to the initial operating condition or the temperature will increase, which contributes to the coil not being able to perform 10

as before, most notably, not being able to reach rated current scenarios. Finally a thermal runaway will make the superconducting coil quench completely, in which case the heat transferred to the neighboring regions and heated liquid helium being expelled, will contribute to the formation and propagation of thermally induced waves.

2.1

The physical description of quench propagation

The physics underlying quenc h propagation are very complex, which includes the coupling effects among strands (copper plus superconductors), stainless steel conduit or jackets, and helium fluid inside the CICC, and the adjacent interaction between layer-to-layer, turn-to-turn, and the impact of unforeseen interstice between epoxy resin or glass fiber and foil of the conductors etc. Because of large coupling effects of CICC, the helium flow inside is under a strong unsteady regime, which is between invsicid high Reynolds number compressible flow and viscous incompressible flow. These features, high Reynolds number, low Mach number and viscous flow, lead to the difficulties in the solution of quench propagation. In the following sections we first give a quick review of CICC, second, we analyze the coupling effects of CICC, and finally describe the thermal and hydraulic evolution of quench propagation. The CICC concept is widely used in superconducting magnets designed and built for other large-scale applications, most notably fusion. In the CICC design, the superconducting cable is placed inside a conduit (or jacket) filled with helium, so that the conductor is not only the main electrical path, but also the helium containment element. This combination of functions allows for more flexibility in the magnet design, with the ensuing potential for simplification and lower cost. There are two types of CICC, one is penetrated by a central, helium filled tube, ITER magnet belongs to this category. Another type, such as the Westinghouse LCP conductor, or the NHMFL 45-T hybrid magnet, does not contain the central pipe. Both types of CICC have similar construction, the strands consisting of copper and superconductor filaments, the stainless steel conduit as a tube for liquid helium flow inside. Figure 2.1 shows the CICC used in the NHMFL 45-T hybrid magnet system.

11

Figure 2.1: Cross sectional views of conductors used for the three subcoils in the superconducting outsert at NHMFL: first stage subcables above and finished CICCs below (provided courtesy of John R. Miller [26])

The cable cross section of coil “A” of the superconducting outsert of the NHMFL 45-T hybrid magnet system is about 79.44 mm 2 , the hydraulic diameter is about 4.45 × 10 − 4 m, however, the wetted perimeter between strands and helium is up to 4.05 × 10 −1 m, and the wetted perimeter between the conduit and helium is about 4.37 × 10 − 2 m. As mentioned in Chapter One, magnets built with CICC are not truly cryostable but rather metastable. Inside the tube the helium inventory available for recovery is limited. Because the time scales for the quench or recovery (tens of milliseconds) and the residence of helium in the conductor (minutes) are different, the replacement of warm helium by fresh cold helium (that could take place in a large bath) may not occur fast enough to promote recovery. Therefore, if a thermal perturbation is strong enough, it may raise the helium temperature beyond the current

12

sharing threshold. The Joule heating of the conductor cannot be extinguished and a quench ensues [2].

Figure 2.2: Schematic of quench propagation in a CICC which is connected to a large helium bath in the end. In the middle region the superconductor has become a normal conductor.

Normal zones grow because heat is transferred along the conductor from resistive region to the still superconducting region. In cable- in-conduit conductors, helium is confined in a CICC and expands when heated, which is different from potted magnets or pool-cooled magnets, so that during the growth of a normal zone, it exhibits powerful longitudinal flows. These flows transport heat along the conductor, and so affect the propagation velocity. And the picture of how a normal zone grows in a CICC is much more complicated than the picture of that in a

13

potted or pool-cooled conductor, and no simple formula for the velocity has been developed yet [2]. At least we can give a qualitative description of normal zone growth based on the coupling effects among the elements in CICC. We know that, immediately upon being heated, the helium opposite the normal conductor begins to expand. At first, the motion of the helium is determined by the heating- induced pressure rise and the inertia of the fluid, friction against the wires playing little role. But this inertial flow regime lasts tens of milliseconds at most. The inertial flow regime plays a role in determining stability because the issue of recovery is decided in tens of milliseconds; but once a quench is in progress, friction of the fluid against the wires quickly becomes the dominating influence. The inertia of the fluid can then be ignored, and the motion of the fluid can be determined from the heating- induced pressure rise and the turbulent friction of the fluid against the wires. The helium which in contact with normalized conductors in CICC is heated to the critical temperature in hundreds of milliseconds by Joule heating alone. After this much time has elapsed the helium being expelled from the normal zone is hot enough to drive cold superconductor with which it comes in contact fully normal. Thus, we expect the helium being expelled from the normal zone to be hot enough to normalize additional conductor in just a few tens of milliseconds. When this happens, the velocity of normal zone propagation equals the local velocity of expansion of the helium. This means that the longitudinal heat conduction is not likely to play a big part in determining the propagation velocity. And this is an assumption in the following modeling [8]. A more theoretical analysis of quench propagation is due to the solution of low Mach number compressible flow. The helium flow in a CICC contains the properties of a relatively uncommon regime of fluid dynamics, the fluid is compressible, and the Mach number is about 0.1, very small compared to that of transonic or supersonic flow in gas dynamics. In this type of internal flow, large length-cross section ratio, finite pressure, temperature, and density gradients, imply the importance of compressibility of helium fluid. There is some critical non-dimensional number, above which, the flow reverses against the direction of the applied pressure gradient causing fluid to leave both the inlet and outlet implying the transition of friction factor and heat transfer coefficients [28].

14

2.2

Modeling of quench propagation

In view of the above analysis and illustration, it would appear that the quench propagation has no rigorous explanation, and the complete theory of quench propagation is unavailable yet. Nonetheless, it is useful as a qualitative picture, for it sheds light on the problem, and provides some limiting cases subject to analytic solution. The comparison of modeling of quench propagation in this section offers us three clues about numerical simulation, each as important as the problem itself. First, the justification of an assumption adopted in modeling, even though it may appear as obvious. Second, the choice of spatial dimensions will depend on the trade-off between accuracy of design and computational cost. Finally, a critical issue, examined in this section, is the choice of mathematical models. What level of complexity is needed to provide sufficient accuracy for magnet design, and what is the impact on cost and implementation? The following sections present most mathematical models available until now.

2.2.1

Zero-dimensional model

The first and simplest approximation which has been proposed, is the zero-dimensional model, in which only local behavior of cable was considered, and a closed, adiabatic system was assumed [2], [3], [45]. In this case all the heat fluxes are assumed to be zero and the helium is at rest. Therefore the mass and momentum equations are not necessary. Only the energy equation needs to be solved. The heat balance is done usually by lumping all the heat capacities of the cable in one. Some assumptions are given as follows, § No temperature gradients in the cross section (the same temperature T); § Adiabatic conditions in the cable; § Zero convection velocity and heat fluxes; § Thermal capacity of the cable is constant

15

Figure 2.3: A simple description of hot-spot calculation model (zero-dimensional), helium is stagnant, and lumped parameters are considered.

Hence we obtain the expression of zero-dimensional model T

t

∫ A ρ C dT = ∫ J

T0

2

ρ Cu ACu dt

(2.2.1)

t0

where the equivalent heat capacity is given by A ρ C = ∑ Ai ρi Ci

(2.2.2)

i

Ai , ρ i , Ci are the cross section area, the density and the specific heat of the i th component (superconductor, copper, stainless steel, helium, and epoxy etc), respectively. The integration form of heat balance (2.2.1) is valid for hot-spot calculation (above the critical temperature) [3]. Because of the strong decrease in the specific heat and density of liquid helium when the temperature rises from superfluid helium to normal helium (especially above 10 K), in practice the process is often assumed to be at constant density. A zero-dimensional model can approximate the transients of short duration and in such cases is in good agreement with experimental data [46]. The advantages of the zero-dimensional model are that it is simple, and it is better approximation for short duration transient problem. However,

16

this simplification can only illustrate a local effect of temperature evolution. For more accurate simulations this model usually gives rise to conservative results with regards to the design of the cable as it neglects the heat sinks around the a normal zone.

2.2.2

One-dimensional model

Most numerical simulations of quench are based on one-dimensional or quasi-twodimensional models, in which helium flow is only taken into account along the channel of the cable. The conductor can be assumed to exchange heat by convection to the flowing helium, the heat diffusion along the length of the cable, parallel to the helium flow, can be retained or neglected. In any case the temperature gradients across the conductor are assumed negligible. For CICC in which copper is used as stabilizer the strands have lower thermal resistance than that due to thermal boundary layer of the helium. Only the spatial dimension along the channel of the cable is included in the governing equations. Basically there are some assumptions applied to 1D model as follows: §

There is heat transfer (convection) between helium and conduit;

§

Thermal diffusion along the length of the cable can be neglected;

§

The convection term of helium can be removed too;

§

The temperature gradient across the conductor cross section is zero, this assumption takes into account the scale across diameter or equivalent diameter to the length of the conductor is very small, usually several order of magnitude;

§

The temperature difference between conductor and helium can be zero;

§

The temperature difference jacket and conductor or helium is different;

§

The boundary condition between conductor and jacket can be adiabatic or nonadiabatic;

§

Friction force was included in 1D model no matter under normal helium or superfluid helium regime;

Other assumptions such as the helium properties (e.g., density, heat capacity, velocity of sound, entropy, enthalpy, etc.) will be discussed in a later section.

17

Figure 2.4: A typical one-dimensional model of quench propagation simulation, usually the cross-sectional length scale of the CICC ( ≈ 0.1 m) is very small compared to its axial length ( ≈ 100-1000 m).

The governing equations for one-dimensional model can be expressed as: §

Helium flow equations in conservation law vector form v v ∂U ∂ F v + =G ∂t ∂x

(2.2.3)

or in a long form

§

∂ ρ ∂ ρv + =0 ∂t ∂x

(2.2.4)

ρv v ∂ ρ v ∂ ρ v2 + p + = −2 f ∂t ∂x Dh

(2.2.5)

∂ ( ρ e) ∂ ( ρ e + p ) v + = h (T jk − THe ) ∂t ∂x

(2.2.6)

Equation of state of helium

e = e ( p ,T )

(2.2.7)

Unlike the ideal gas and barotropic fluid there are analytic expressions relating pressure, temperature and internal energy of liquid helium. The thermal properties data of helium come from NBS reports [47], [48]. §

Well-posed initial and boundary conditions

18

(2.2.4)-(2.2.7) are the conservative form for one-dimensional helium flow problem, before moving ahead we should make some observation about the system of equations. First, the system can also be written as non-conservative form for a different set of dependent variables, however, the advantages of conservation form are that numerical schemes can improve the stability in time stepping, especially across physical discontinuities. Second, strands are thermally coupled with helium, we assume that strands keep the same temperature as he lium, and only the wall temperature (jacket temperature) differs from that of helium. Third, one-dimensional internal channel flow is in turbulent regime, in which turbulent friction effect is compensated by Darcy friction coefficient of laminar or turbulent flow in the momentum equation [49]. This approach applies a great shortcut to the solution of Navier-Stokes equation from laminar flow to turbulent flow, and circumvents the solution of the boundary layer equations of the rough wall region. In this model we omitted the effects such as boiling helium, and phase change, because the helium flow is in supercritical regime during the quench propagation. Hence we obtain a relatively compact one-dimensional model. A variant form of this compact 1D model was given by Shajji and Freidberg [18]. Apart from the system of equations (2.2.1), the complete governing equations also consist of conductor heat balance and conduit heat balance equations. The energy equation for the conductor is given by:

ρ st C st

∂ THe ∂  ∂ THe  k  + S ( x, t , THe ) = ∂t ∂ x  st ∂ x 

(2.2.8)

In the right hand side of equation (2.2.8) heat transfer between strands and conduit was neglected, in fact the value is very small compared to that of strands with helium. The source term is a function of the temperature of the strands, it can be written as:  I( t )   ρe S(x, t, Tst ) = J 2 ρ e =    A Cu  2

(2.2.9)

The energy equation of conduit is given by

ρ jk C jk

∂ T jk ∂t

=

∂ T jk  h jk Pjk ∂   k jk − ( Tw − THe ) ∂ x ∂ x  A jk

(2.2.10)

In the right hand side of equation (2.2.10) no external source term was included, because stainless steel has very high electric resistivity, the Joule heating is very small. Again no heat transfer between conduit wall and strands in the RHS of (2.2.10) is considered.

19

Hence (2.2.4)-(2.2.10) represents a complete one-dimensional model of quench propagation.

2.2.3

Quasi-one -dimensional model

In quasi-one-dimensional model the difference of temperature within the cross-section is also included. Therefore, the coupling terms between helium and strands, helium and conduit, strands and conduit (three-element system) should appear in the right hand side of the corresponding energy equation, respectively.

Figure 2.5: Schematic of quasi 1D CICC in which every component is represented by a point, and there is a heat-link between two components.

20

The basic governing equations consist of continuity, momentum and energy equations for helium, cable energy equation, and jacket energy equation v v v v v U t + F (U )x = D + G

(2.2.15)

m & ρ   0 v   v v  2  v   U =  m&  , F (U ) =  m& / ρ + p  , D =  0   e  m& (e + p ) / ρ  Q       c

(2.2.16)

0   v   G =  − f ρ u u / 2d h    Qh  

(2.2.17)

Qh = β He, st (Tst − THe ) + β He, jk (T jk − THe )

(2.2.18)

∂  ∂T  Qc =α k(T)   ∂x   ∂x 

(2.2.19)

1/ 3

   

Where α = 1 is for superfluid helium, α = 0 for normal helium fluid. The boundary conditions will be discussed later in this chapter. Again, the liquid helium state equation is added to make the system closed:

ρ = ρ (e, p )

(2.2.20)

The equivalent heat transfer coefficients are given by: β He, st = hst pst / AHe

(2.2.21)

β He, jk = hjk p jk / AHe

(2.2.22)

The heat transfer coefficients are very important in determining system behavior, and standard correlations can be found in [50], [51], [52], [53]. All the variables without subscripts or superscripts refer to liquid helium. We simply write the cable energy equation by introducing the equivalent heat transfer between helium and conductor.

e&st = β st,He (THe − Tst ) + s(T, x, t) + Acu∂x (kcu∂ xTst )/ Ast

(2.2.23)

Where e&st = ρst cst∂ tTst β st , He = hst p st / Ast

Ast = Acu + As / c 21

Ast ρstc st = Acu ρcuccu + As / c ρ s / c cs / c The external source term, s(T , x, t ) , is the sum of Joule heating (after quench) and the external heat pulse (to initiate the quench), the units are W/m3 . Similarly we can give conduit energy equation in the following, ρ jk c jk ∂ tT jk = β jk , He (THe − T jk ) + s jk (T , x, t )

(2.2.24)

Where β jk , He = h jk p jk / A jk

In the right hand side of (2.2.24), a source term (Joule heat) has been included, in fact the external heating sources are usually located in the strands and the Joule heating in the jacket is negligible due to its high electrical resistivity at cryogenic temperatures. In the three components of the CICC the jacket energy equation is the simplest and most stable, it is an initial value problem without specific boundary conditions. The quasi-one-dimensional model keeps the advantages of simplicity and compactness. However, it does not consider the thermal boundary layer. Existing codes, such as GANDALF [17] and POCHI [21] are based on quasi-one-dimensional models. Since the liquid helium flow during quench is usually in a turbulent state in CICCs, the coolant temperature is assumed to vary drastically in the thin layer around the strand surface and is a constant in the main flow region, a more precise model was proposed to increase the accuracy in the calculation of transient heat transfer coefficient between liquid helium and strands [22], [50]. In this model the coolant temperature is numerically calculated on twodimensional grids only in this boundary layer, assuming constant helium properties in a direction perpendicular to the strand surface and neglecting the effect from coolant convection and heat conduction in the conductor axial direction. The heat transfer from the strand to the coolant is in transient state, in which heat conduction in the very thin layer of the coolant around the strand is dominant at the beginning of the heating. There is a combination between transient heat transfer with steady heat transfer coefficient. In this calculation constant heat flux from the strand to the coolant is assumed. The governing equations consist of (2.2.17)-(2.2.20), (2.2.23), and (2.2.24). The difference is the thermal boundary layer equations,

ρ He C He

 ∂ 2Θ ∂Θ 1 ∂ Θ  = k He  − 2 ∂t d h / 2 − y ∂ y  ∂ y

(2.2.25)

22

ρ He C He

 ∂ 2Φ ∂Φ 1 ∂ Φ  = k He  − 2 ∂t d h / 2 − y ∂ y  ∂ y

(2.2.26)

(2.2.25) is the thermal boundary layer equation for strands, and (2.2.26) is for conduit wall. The boundary layer temperatures Θ and Φ are closed by supplement of boundary conditions and initial condition. These conditions are given by Θ ( y = 0, t ) = Tst , Θ ( y = δ , t ) = THe Φ ( y = 0, t ) = T jk , Φ ( y = δ , t ) = THe

The boundary conditions at the surface of the strand and the conduit wall can be replaced by Neumann boundary conditions as follows, − k He

∂Θ ∂Φ | y =0 = qst , − k He | y =0 = q jk ∂y ∂y

where the q st and q jk are the heat flux between the strand and helium, jacket and helium, respectively. This 2D model can demonstrate a more detail information near the wall of helium heat transfer and improve the accuracy of the calculation of the temperature, however, more dependent variables are introduced to the system.

2.2.4

Three-dimensional model

Many critical phenomena of helium fluid flow in CICC, such as two-phase change, discontinuities, moving front, and turbulence, are essentially nonlinear, and they also exhibit extreme disparities in time scales. While the thickness of a moving discontinuity is on the order of millimeters, on a length scale of CICC its thickness is approximate zero. Consequently, mathematical models with varying degrees of simplification have to be introduced in order to make computational simulation of flow feasible, and produce viable and cost-effective methods. On the other hand, viscous simulation at high Reynolds number of 3D problem requires vast computational resources. In that context, we will give a quick review of three-dimensional models of the quench propagation problem.

23

Figure 2.6: The cross section of 3D model: helium, strands (copper plus superconductors), the conduit wall, and the insulation between adjacent layers or turns.

A typical schematic diagram of three-dimensional model of CICC is shown in Figure 2.6. Layer-to- layer, turn-to-turn, and cross-section difference are included in this model. The coupling effects such as epoxy, insulation winding, and supplemental factors are also components of this model. On 1D model an adiabatic condition is imposed on the outside boundary except the two ends, this is not the case for three-dimensional model. Looking at the cross-section of a typical CICC, the superconducting strands and the liquid helium cooling channels are embedded in a composite matrix, whic h is formed by the stainless steel jacket (conduit) and the insulation layer (turn-toturn, layer-to- layer). The effect of heat transfer in the transverse direction tends to short-circuit the quench propagation along the conductor which cause the pre-heating and the birth of the

24

normal zones in regions which have not yet been reached by the hot helium front. That is, the heat conduction will initiate the quench in the neighboring conductors, i.e., the perturbation happens at a local region [14]. Considering the geometry of superconducting magnets, the jacket, epoxy insulation and resin layer can be equivalent with a connection between two layers or two adjacent turns. Hence, the three-dimensional quench propagation can therefore be reduced to a set of one-dimensional problems, the mass, momentum, and heat transport in the channels with a coupling node or interface link. These thermal links exist between turn-to-turn, and layer-to-layer. Compared that 3D model will call for a sophisticated approach, therefore an expensive computational cost, 1D model is easy to implement and modify for special magnet systems. Another advantage of a set 1D problems is that spectral element method or finite element method can be successfully. For steady-state non- linear problems it is possible to use a reasonable approximation of the thermal resistance by performing a further subdivision of the thermal link in portions with nearly constant thermal conductivity. In this aspect, the steady state solution can be regarded as asymptotical solution of transient case. The dynamic equivalence will depend on the physical time constant of the transient. There were some experimental tests for quasi-two-dimensional models (e.g., QUELL) [54].

2.3

Boundary conditions

In this study, quasi-1D que nch propagation model is introduced. The system of equations of one-dimensional compressible viscous flow in a CICC must be supplemented by some initial and boundary conditions to make the system well posed. Unlike the Euler equations, we cannot find the characteristics of the system for quench propagation. So there is no way to impose boundary conditions along characteristics, therefore, the only possible approach is based on empirical considerations. So far there is no theoretical guideline to select boundary conditions of nonlinear hyperbolic or convection-dominant systems. In a coil each cable is connected to large inlet and outlet, which provide constant pressure and temperature ambient conditions. The helium flow is usually in subsonic conditions, at Mach numbers of the order of 0.1 during a quench and well below that during normal operation, hence,

25

a simple way to impose the dynamic boundary condition is to assume the static pressure to the pertaining environment, for example, at the end of the tube the static pressure is set equal to the reservoir stagnation pressure or bath pressure. This assumption amounts to the loss of one velocity head for the flow entrance or exit [12]. Similarly, the energy equation can also be imposed at the inlet. The enthalpy can be set equal to the inlet reservoir enthalpy (upwind condition). The enthalpy of a reservoir at the outlet (downwind) does not enter into the problem. We assume that all heat conduction into the ends of the CICC is neglected, including that into the helium. The outer surface of CICC, which is also the outer epoxy surface, can be either adiabatic or isothermal. If it is isothermal, the surface temperature is set equal to that of the reservoir with the highest pressure [12]. The equations of conductor and conduit are of parabolic type. Hence boundary conditions can be simply imposed by a Dirichlet type boundary condition or Neumann boundary condition on the end. Another empirical approach to dealing with the boundary condition is to introduce valve effects on the end. Valve or manifolds are used to decrease the pressure near the node or element at the end. This approach can avoid introducing two-dimensional or three-dimensional expansion of helium in the end by one-dimensional valve to dynamic equivalence of the pressure loss and the change of flow velocity in the ends [10], [13]. A more strict approach to give the boundary conditions is due to a linearized characteristic analysis [17]. Finally, we need to impose Stefan conditions if the helium initially is under superfluid region in which a He II/He I phase-change moving boundary exists. The details about moving boundary conditions are referred to Crank [55].

26

CHAPTER THREE THE ENERGY BALANCE MODEL AND NUMERICAL SOLUTION

In this chapter a simple energy balance model is proposed, numerical solution is given based on this model. The advantage of this model is simplicity yet it can correctly evaluate, analyze, and estimate the temperature evolution in large-scale superconducting magnet systems with subcooled superfluid helium cryostat. Especially the superfluid helium and normal helium fronts are tracked very accurately even though the difficulties due to the exceptional properties of superfluid helium. The numerical simulation is focusing on NHMFL 45-T hybrid magnets systems.

3.1

Thermal analysis of subcooled superfluid helium in CICC

Some large-scale superconducting magnets use superfluid helium as coolant, because superfluid helium has unique properties, such as very high thermal conductivity, free of viscosity. We develop a compact model based on only energy balance in quench propagation under superfluid helium regime. The aim is focusing on the NHMFL 45-T Hybrid Magnet Systems. The initial physical conditions: static superfluid helium (helium II) at 1.8 K, and both ends are clamped to constant bath temperature and stagnation pressure of helium reservoir. The external source terms are the AC losses during the operating current ramping up from zero to the maximum value, and the index heating because of the deviation of superconductor from the ideal state [45]. From experimental observations, the superconducting coil behavior suggests a thermal runaway leading up to a quench. The external heat driving the helium from superfluid state to normal fluid state, and finally arrive at the boiling phase which drastically decrease the heat transfer coefficient between helium and the conduit. Figures 3.1-3.2 show the schematic of a compact model. Some assumptions have been made: 27

1)

Thermal conduction along the metal in the CICC is negligible when the helium is in the superfluid regime.

2)

There is no heat exchange with the liquid helium at the bottom and the top of the coils.

3)

The temperature at the two ends of each layer (the helium bath) is assumed to be constant (1.8 K).

4)

Helium, strands, and jacket are assumed to be at the same temperature at all times. This is a good approximation for the slow ramping up process of interest.

5)

Only helium heat capacity has been taken into account in He II regime.

6)

Liquid helium is static, no convection is included in the model.

7)

Helium thermal conductivity in He I regime has been neglected, and copper thermal conductivity is included once helium transitions above the lambda point.

8)

Layer-to- layer heat transfer is included. The model considers heat conduction through a series of comprising the thickness of two steel jackets, one layer-to-layer insulation (epoxy) thickness (e.g., a one and a half dimension model). The model also includes the Kapitza effect to account for interfacial resistance in low temperature thermal systems (a Kapitza resistance term is included for each liquid-solid or solid-solid interface). The model assumes constant heat transfer coefficient, but heat transfer rate is not constant as actual temperature differences are used.

28

Figure 3.1: Schematic of the one-and-a-half dimensional model used to solve the energy balance for each layer

Figure 3.2: Schematic of the NHMFL 45-T superconducting outsert. There are six layers in the coil “A”, the layer 1 is at the highest magnetic field.

29

The following system of equations is used to mathematically describe the heat balance model presented above. The heat conduction relation for superfluid helium is given by the GorterMellink law: q x = − K ( dT / dx ) 1/ 3

(3.1.1)

where K = ( f − 1 (T , p)) 1/ 3 is the equivalent thermal conductivity of superfluid helium that can be written as a function of temperature and pressure [59]. The energy balance of each layer, adding the corresponding heat generation terms becomes [43]: ρ He cHe

1/3 ∂T ∂   ∂T   = K     + q in + q ∂t ∂x   ∂x  

AC

+ q Across

(3.1.2)

where qin is the index heating, q AC represents the AC losses (only present during current rampup), q Across is the heat transfer between two adjacent layers, ρ c represents the heat capacity of He-II. For the sections where the temperature goes above the lambda transition, the energy balance becomes (Fourier law applies): A′ρ ′c ′

∂T ∂  ∂T  = ACu K Cu + q in + q AC + q Across  ∂t ∂x  ∂x 

(3.1.3)

where KCu is the thermal conductivity of copper (we neglect heat conduction through He-I). A′ρ ′c′ represents the equivalent heat capacity of He I plus copper, which can be written as

follows: A′ρ ′c ′ = AHe ( ρ c )He + ACu (ρ c )Cu

(3.1.4)

The boundary conditions for this model are: T ( 0, t ) = T (l , t ) = Tb

(3.1.5)

T ( s (t ), t ) = Tλ

(3.1.6)

x = s (t ) is the location of the lambda transition, Tb is the clamped helium bath temperature (1.8

K), and Tλ is the superfluid transition (2.176 K under saturated vapor pressure). The index heating is given by Miller et al [42]:

30

E I  J op   qin = 0  A  J c 

n

(3.1.7)

The AC losses, q AC , are computed using standard expressions for hysteresis and cable coupling losses, with the corresponding coefficients benchmarked by both design wire and cable parameters, and global calorimetric data for AC losses in coil A taken during commissioning [56]. The coupling losses in the cable due to transverse field: pc =

2τ & 2 B Ast µ0

(3.1.8)

where τ is the conductor effective time constant and Ast is the total strand cross section in the cable. µ 0 is the permeability of free space. Hysteresis losses in the superconductor due to transverse field are given by the piece-wise formulae [57]: ph =

2 J c B& D eff As / c 3π

ph =

2 3B2 J c B& Deff As / c 2 3π Bp

B > B p = µ 0 J c Deff

(3.1.9)

B < B p = µ 0 J c Deff

(3.1.10)

where J c is the superconductor critical current density, Deff is the superconductor filament effective diameter, and As / c is the total superconductor cross-section in the cable. Eddy current losses in the jacket are very small compared to the other two mechanisms and can be negligible. Hence we have: q AC = pc + p h

(3.1.11)

Finally, for the last term in Eqn. (3.1.2), the equivalent layer-to- layer heat transfer coefficient is given by: l glass−epoxy 1 l 1 = ss + +∑ heq k ss k glass−epoxy i hk

(3.1.12)

31

where k ss and k glass− epoxy are thermal conductivity of stainless steel and insulation glass-epoxy, respectively, hk is the Kapitza conductance at each interface, there are eight Kapitza conductances between two adjacent layers (see [58] for solid/solid Kapitza conductances). l ss , l glass−epoxy are the CICC conduit thickness and the thickness of epoxy between two layer in

coil “A”, respectively. Two fitting formulae are used in the calculation of the layer-to- layer heat transfer coefficient [57]: K epoxy = −.00633 + 0.016T

K ss = 0.0131 + 0.0364T

(3.1.13)

W / mK

W / mK

(3.1.14)

The Kapitza conductance between helium and solid interface is approximated by a cubic expression, which can be written as follows [1], [58], [59], [60],

hk = c Tw2 + T f2 (Tw + T f

(

)

)

W / m2K

(3.1.15)

where c is a coefficient that depends on surface geometry and materials being paired at the interface, and ranges from 100 to 400. In our calculations a value of 200 was chosen for all interfaces. The equivalent heat transfer coefficient between two layers is in the range from 0 (adiabatic condition) to 20 W/Km 2 (good contact). Hence the heat exchange term can be written as:

q& L′′ = − heq (TL +1 − TL ) W / m2 q cross, L = (q& ′L′−1 − q& L′′ ) / Scond , r

(3.1.16)

W / m3

(3.1.17)

where q& ′L′ is the heat flux between layer L and L+1, q cross, L is the heat source term to layer L due to heat transfer from adjacent layers, and Scond , r is the cross section width of cable in the radial

32

direction. We assume that q& ′L′=0 and q& ′L′=6 are equal to zero, in other words, the innermost and outermost layers have no heat transfer with the resistive insert or the inter-coil insulation. There exists two symmetrical moving fronts between He II and He I, across the front there is a discontinuity of heat capacity and other thermal properties, this is a key problem for the numerical method to handle. On the moving front the fluid is locked at the lambda temperature and Stefan condition is satisfied for the heat flow balance [55]. Equations (3.1.2)-(3.1.15) represent a one-dimensional, two- fluid model from which the temperature evolution in the CICC can be calculated. It is a non- linear and stiff system of partial differential equations requiring numerical solution.

3.2.

Numerical algorithm

To solve (3.1.2), a five-point MacCormack explicit scheme is used [61]. This finite difference scheme can be expressed as forward in time, and prediction-correction in space, five points are needed to compute T jn +1 (these being T jn− 2 , T jn−1 , T jn , T jn+1 , T jn+ 2 ). Two additional points (ghost points) are needed at the boundary to start the algorithm: T−1n = Tb

(3.2.1)

T−n2 = Tb

(3.2.2)

Such five-point scheme is particularly suited to improve the stability and accuracy while solving a non- linear system. A hybrid scheme (Crank-Nicolson) is used to discretize (3.1.3) [33]. The corresponding form is a set of linear algebraic equations that can be expressed as a tri-diagonal matrix form. Both finite difference methods introduce a constant time step and uniform mesh, except in the neighborhood of the moving front. The time step was set from 10 −4 to 10 −2 sec, and the spatial mesh size can be varied from layer to layer, but the number of elements in each layer is kept constant, typically 200 elements.

33

There are several methods to handle the moving front [55]. The present model employs a fixed mesh method to track the lambda transition. The temperature T is continuous along the x direction (although its first derivative is not continuous), three points are used to construct Lagrange base functions that will be the weights used to calculate the front. For example, j=0,1,2, front is between x 0 and x1 , then we have: 2

T ( x) = ∑ l j ( x )T ( x j )

(3.2.3)

j =0

where l ( x j ) is a Lagrange-type base function, T ( x j ) is the temperature at mesh point j. Combining the MacCormack scheme and hybrid scheme we can solve this two- fluid model. MacCormack explicit scheme to solve He II parabolic equation can be expressed as follows:

( )

qin = q Ti n =

( )

( ) (T

K Ti n−1 + K Ti n 2(∆x )

1

n

i

− Ti n−1

)

1

(3.2.4)

3

3

Prediction process:

~ Ti = Ti n +

1

(ρ c )

n i

∆t n qi +1 − qin + f i n ∆ t ∆x

(

)

Correction process: ~ ~ ~ K Ti + K Ti+1 ~ ~ ~ q i = q Ti = Ti+1 − Ti 1 2(∆x ) 3

( )

Tin+1 = Ti n +

( ) ( )(

)

(3.2.5)

1 3

(3.2.6)

1 ∆t n qi +1 − qin + q~i − q~i−1 + f in ∆t n ( ρc )i 2∆x

(

)

(3.2.7)

where f i n is the source terms in (3.1.2) and (3.1.3), hence we have,

f =

qin + q AC q + Across AHe ρ He c He ρ He c He

(3.2.8)

We use a semi- implicit scheme (Crank-Nicolson scheme) to solve the classic Fourier law in He I region. The two systems are transformed to non-dimensional equations before applying finite difference schemes. MacCormack explicit scheme has second order accuracy in space and time for steady-state problem, first order in time for unsteady problem. The C-N scheme keeps second order accuracy in space and time.

34

3.3.

Model results and discussions

The modeling is intended to elucidate the importance of index heating and layer-to- layer heat transfer vis-à-vis axial heat transfer through a column of He II (layer end cooling). Cases are considered with and without layer-to- layer heat transfer (heat transfer vs. adiabatic). Also, cases are considered in which the index number n=15 and n=5 (possible conductor damage following an unprotected quench). The experimental runs being modeled are summarized in Table I. In most, but not all, cases the current is ramped up at about 5 A/s. The evolution of AC losses during a 5 A/s ramp is shown in Fig. 3.3, Hysteresis being the dominant term; the time constant for coupling losses was assumed to be 10 ms our calculations (consistent with the range of measured values in short sample specimens). Since the current is changing for the first 2000 s, the AC losses are zero after that. In general faster current ramps mean higher AC losses in the coil, and correspondingly an increase on the risk of a quench. Even if the temperature does not reach the current sharing temperature, the temperature margin will be decreased. If the coil temperature rises up to the lambda point (2.176 K under saturated vapor pressure), the liquid helium becomes normal, and helium heat conductivity is very poor, and the coil stability margin is consequently reduced. The following sub-sections describe the different cases examined with the model.

35

Table 3.1 Record of high-current operations, superconducting outsert

3.3.1. Adiabatic conditions

The first set of cases modeled assumes no heat transfer between layers (adiabatic conditions) [43]. These cases roughly correspond (as it will be seen later) to the ramps following the unprotected quench (Run 12 and later in Table I). From (3.1.8)-(3.1.10) we know that AC losses vary during the current ramp up, the rate of change of magnetic field is proportional to the ramp rate, and critical current density is a function

36

of B and T. Hence dominant hysteresis losses have a peak at about 500 s into the ramp (Fig. 3.3). Index heating is shown in Fig. 3.4 (for L1) and Fig. 3.5 for n=15 and n=5 (damaged conductor). It can be seen that while the current is a small fraction of the rated value, index heating is negligible. However, towards the end of the ramp (and soon after rated current is achieved), the index heating is significant (it will be shown that the temperature at the center is significantly higher than 1.8K, thus decreasing the effective critical current, which leads to increased index heating). It should be noted in Figs. 3.5 and 3.12 that there is a break point in the curve for n=15 (absent in the n=5 case). This can be explained by noting that before the break point the helium temperature is increasing by virtue of AC losses and index heating, and the break point corresponds to the end of the current ramp, at which point the AC losses are zero (so all subsequent temperature increase is due to index heating alone). In the case of n=5 index heating takes off well before the end of the current ramp so that the break point is not present. The temperature evolution for L1 is shown in Figs. 3.6, 3.7, and 3.8. Note that as a result of AC losses the bulk of the fluid in the channel becomes He-I well before the end of the ramp. Initially the temperature in the He-I region is rather uniform as a result of the AC losses being mostly uniform along the x-direction (see Figs. 3.3 and 3.7). After the current reaches its steady-state value (10 kA), index heating is the only energy deposition into the conductor. By then the center point temperature is close to 2.75 K (significantly higher than the design value of 1.8 K), meaning that Jc at the center is lower than the design value and the operating margin has been reduced. Index heating continues to deposit energy even after the conclusion of the current ramp, and the temperature continues to increase (leading to further decrease in Jc, which leads to increased index heating and a thermal runaway situation with the temperature reaching 4.2 K and then current sharing a few minutes after the conclusion of the ramp (see Figs. 3.7-3.11).

37

Figure 3.3: Evolution of total AC losses (hysteresis plus coupling) in Layer 1/Coil A during a 5A/s current ramp.

Figure 3.4: Evolution of index heating in Layer 1/Coil A. Ramp rate is 5 A/s for 2000s, and then current stays constant (n=5).

38

Figure 3.5: Index heating for Layer 1/Coil A during and following a 5 A/s current ramp (layer is adiabatic, insert coil is turned off).

Figure 3.6: Peak temperature in Layer 1/Coil A during a 5 A/s current ramp (layer is adiabatic, insert coil is turned off).

39

Figure 3.7: Temperature evolution of Layer 1/Coil A during a 5 A/s current ramp (layer is adiabatic, insert coil is turned off).

Fig. 3.8: Evolution of the temperature of Layer 1/Coil A during 2.5 A/s current ramp, and insert coil is turned off.

40

Figure 3.9: Temperature distribution in the six layers of coil “A” at t=2033s (end of 5 A/s current ramp). Heat transfer between layers has been turned off. The maximum temperature difference between the innermost and outermost layers is about 2.2K.

Figure 3.10: Peak temperature in Layer 1/Coil A for different current ramp rates (layer is adiabatic, insert coil is turned off). There is no thermal runaway due to index heating for a 2.5 A/s current ramp to 10 kA.

41

Figure 3.11: Temperature, operating current density, and critical current density at the center of Layer 1/Coil A during a 5 A/s ramp to 10 kA (layer is adiabatic, insert coil is turned off, index number n=15). The temperature has exceeded the lambda transition and the actual critical current is below the design value/margin, leading to a thermal runaway due to index heating.

These simulations show that in the absence of interlayer heat transfer a 5 A/s ramp to 10 kA takes the coil to a state in which index heating alone after the ramp is enough to create a thermal runaway in a few minutes, especially with a damaged conductor (low “n”), which is consistent with observations (for example, see Runs 10-12 in Table I). Finally we can compare the temperature evolution of coil “A” for different ramp rates, for example, 5 A/s and 2.5 A/s. Fig. 3.10 shows that 2.5 A/s is a safe ramping rate, even in the absence of layer-to- layer heat transfer. The temperature of L1 increases, but not to the point in which residual index heating will lead to a thermal runaway.

42

3.3.2. Non-adiabatic conditions

There is a possible explanation for the different thermal behavior of the superconducting coils before and after the unprotected quench (Run 11, Table I), and it is the presence of good interlayer heat transfer before the unprotected quench with contact lost due to epoxy delamination after the quench event. The second set of cases simulated corresponds to the one and a half dimensional problem in which layer-to-layer heat transfer was incorporated to the heat balance of each layer (non-adiabatic conditions). Theses cases are reflecting the operations of the superconducting outsert under good contact between turn-to-turn and layer-to- layer (see Runs 7 to 11 in Table I). We will assess the temperature margin and stability of the system under nonadiabatic conditions, and resistive insert coil turned off. Figs. 3.12-3.14 show the calculated thermal behavior of L1, while Fig. 15 indicates that once layer-to- layer heat transfer is included all six layers track each other and remain at about the same temperature. We see from Figs. 3.123.13 that L1 is still prone to index heating runaway for a 5 A/s ramp rate, even if the layers are in good thermal contact (see also Fig. 3.14). Fig. 3.15 shows that the temperature difference between L1, L3 and L6 is quite small, at least during the first part of the ramp, once heat transfer is considered. This indicates that the layer-to-layer heat conduction characteristic time is short compared to the current ramp. As built, there should be good thermal contact between layers in the coil (corresponding to 1.64 mm of stainless steel and 0.6 mm of epoxy resin per CICC turn). AC losses in layers 2, 3, etc. are progressively lower than those in layer 1, so that the temperature in L2 should always lag with respect to L1, providing a sink for some heat removal by conduction. The present model, (2.2)-(2.3) and (2.12), includes an equivalent transverse heat transfer coefficient to solve heat balance equations simultaneously in all six layers, to better approximate the actual operating environment. Consider one layer-to- layer insulation thickness (1.2 mm) and two jacket thickness (3.28 mm), we can estimate that the glass-epoxy thermal conductivity is about 0.025 W/m K at T=2 K, and stainless steel thermal conductivity is about 0.086 W/m K at T=2 K. We know that six interfaces exist so Kapitza conductance should be included in the equivalent heat transfer coefficient. The magnitude of Kapitza conductance is about 770 W / m 2 K . Meanwhile, the

43

approximate value of glass-epoxy heat conductance is 20 W / m 2 K , the value for steel is 27 W / m 2 K , hence the heat resistance between two adjacent layers of coil A is dominated by

insulation and steel, so the equivalent heat transfer coefficient is closer to 20 W / m 2 K (for example, see [Incr90]). This order of magnitude of layer-to- layer heat transfer coefficients can explain why the temperature difference of adjacent layers is not large. We know from the simulation that layer-to- layer has a maximum difference 0.01~0.02 K, which is much smaller than what we had estimated before for adiabatic conditions (see previous sub-section and [43]). Debonding of the epoxy from the conduit as a result of thermal strain sustained during the unprotected quench could then be a possible explanation for the observed quench behavior before and after the event (debonding was observed in mechanical tests of a model coil [Wals98]). When heat transfer between layers is not present, even the small amount of index heating could be enough to upset stability and lead to a thermal runaway. This is especially true for L1 under the worst situation without good heat transfer to the nearest layer and a degraded nvalue.

Figure 3.12: Index heating for Layer 1/Coil A during and following a 5 A/s current ramp (heat transfer between layers is included, insert coil is turned off). Even if heat transfer between layers is included, a 5 A/s current ramp leads to thermal runaway due to index heating.

44

Figure 3.13: Comparison of index heating in Layer 1/Coil A for adiabatic and non-adiabatic conditions (5 A/s current ramp, insert coil turned off).

Figure 3.14: Temperature evolution of Layer 1/Coil A during a 5 A/s current ramp (layer is noadiabatic, insert coil is turned off).

45

Figure 3.15: Temperature distribution in the different layers of coil A during a 5 A/s ramp and at t=500 s. Heat transfer between the layers is included (non-adiabatic conditions) and the insert is turned off.

3.3.3. Influence of the resistive insert coil

The third set of cases shown in the following figures examines the influence of the resistive insert coil on the thermal stability of superconducting outsert. The applied magnetic field before and after the insert coil is turned on is given in Fig. 3.16, which shows the magnetic field distribution on L1 and L6 of coil A at the end of the current ramp (10kA). When the insert is energized, its contribution to the magnetic field on the outsert is negative, that is the applied field on the superconducting layers decreases (Fig. 3.16). As a result, the critical current increases and so does the operating margin. This is reflected in Fig. 3.17 showing that when the current ramp is done with the insert turned on, there is no thermal runaway as the index heating is negligible (similarly, if the outsert is ramped up and then the insert is turned on, the net effect is to turn off index heating, and preventing thermal runaway). Following a ramp at 5 A/s the coil would reach a maximum temperature of about 2.25 K (including layer-to- layer heat transfer, all layers would be at about the same temperature, see Fig. 3.18) After the insert is turned on the critical current density increases, the index heating source 46

becomes negligible, and in the absence of other AC losses (constant current), the coil would eventually return to 1.8K thanks to layer end cooling, as will be discussed in the following subsection.

Figure 3.16: Magnetic field distribution along Layer 1 and 6 of Coil A at the end of a 5 A/s ramp (t=2000 s, rated current of 10 kA) showing the influence of the insert.

47

Figure 3.17: Peak temperature in Layer 1/Coil A during and after a 5 A/s current ramp showing the influence of index number “n”, and lsayer-to-layer heat transfer (adiabatic or non-adiabatic conditions). When the insert coil is turned on there is no thermal runaway due to index heating.

Figure 3.18: There is little effect of heat transfer conditions (adiabatic vs. non-adiabatic) in the peak temperature for Layer 1/Coil A when the insert coil is turned on.

48

Figure 3.19: The initial temperature distribution of Layer 1/Coil A used to calculate the recovery time under adiabatic conditions.

3.3.4. Recovery of the outsert coil temperature

The final case considered with the present model is to assess the thermal recovery time of the outsert coil following a 5 A/s current ramp and subsequent turning on of the insert. A typical temperature distribution in L1 following such ramp scenario is given in Fig. 3.19. The model can now be applied to track the He-I/He-II interface and estimate the time to recover the original helium bath temperature (1.8 K). The thermal recovery process is shown in Fig. 3.20, which indicates that cold-end heat removal through the column of He-II would bring the peak temperature in L1 from 2.43 K to 1.825 K in about one hour. This result is in rough agreement with a previous and more simplified calculation [9], and better matches the experimental observations of coil thermal recovery as seen from the refrigerator operational data following a current ramp. Finally, note that Fig. 3.19 indicates the presence of a strong temperature gradient in the He-I

49

region near the superfluid transition. Our model assumes the fluid to be completely stagnant, and it does not account for convective effects (induced flow). In reality the temperature profile, and evolution of the superfluid interface, are likely affected by fluid transport in the He-I regime even when the bulk remains stagnant. This is an area for future model enhancement.

Figures 3.20: The thermal recovery time of Layer 1/Coil A under adiabatic conditions is approximately 1 h (due to end- layer cooling through a column of superfluid helium).

50

3.4

Conclusions

A transient heat transfer model has been implemented and subsequently improved to study the quench behavior of the NHMFL 45-T hybrid outsert. The model not only tracks the temperature evolution of the column of helium in the CICC, including the transition between He-II and He-I regimes, but also includes layer-to- layer transverse heat transfer. The simulations show tha t if heat transfer between layers is shut off (e.g., as in the case of epoxy debonding), then the innermost layer of coil A will reach 4.2 K and current sharing soon after completion of a ramp to 10 kA at 5 A/s, index heating being the driver for thermal runaway even after the current is kept constant. That ramp rate of 5 A/s is close to the critical value. It is also shown that increased index heating as a result of some cable damage following the unprotected quench could explain the observed coil behavior. The recovery calculation indicates that the stability margin depends on the ramping rate and good heat transfer between layer-to- layer and layer-to-helium bath. It is also shown that for the design ramp rate of 2.5 A/s the coil could reach operating temperature with enough margin to avoid thermal runaway. Future work will include model enhancements to better capture the evolution of the He-I/He-II front, and the inclusion of terms to account for convective effects in He I where temperature gradients are large (enhanced mixing). The merits of this model are its simplicity and ability to correctly track the temperature evolution of the column of helium in the CICC. It also can reflect the moving fronts between superfluid helium and normal helium regimes. Ho wever, the drawbacks are obvious, it cannot predict the maximum pressure, and temperature difference between different components of the CICC, the calculation of recovery is not good enough, and most importantly, it lacks the flexibility to embed convective effects in practical quench propagation evolution. All these will be the impetus to apply the numerical model discussed in the following section.

51

CHAPTER FOUR NUMERICAL SIMULATION OF QUENCH PROPAGATION AT EARLY PHASE BY HIGH ORDER SCHEMES

In this chapter the quasi-1D quench model is used to analyze the supercritical helium flow and coupled heat transfer in CICCs. Liquid helium flow in CICC is in the regime of low Mach number and high Reynolds number flow. The quench propagation is governed by the Euler equations with the addition of large non- linear source terms for modeling turbulent friction forces and heat exchange between helium and conductors or the conduit wall. It is a nontrivial issue to solve this problem very effectively. When discretizing the system of equations of quench propagation, the presence of the nonlinear convective term does not only increase the physical complexity of the fluid flow, but also complicates the design of numerical algorithms. The discrete system becomes nonlinear and nonsymmetrical, and sophisticated numerical techniques have to be applied to deal with these difficulties. The characteristics of the helium flow will change during quench propagation. At the early phase the system is governed by hyperbolic equations, which is similar to acoustic waves. However, it demonstrates more parabolic characters than hyperbolic ones at the later stage. Base on this observation, a highly accurate numerical method is applied to solve the waves equations at the early stage of quenc h propagation. One of our goals in this study is to search for an effective solver for this compressible subsonic flow and coupled heat transfer problem. We choose a high-order method in space to deal with quench propagation simulation by virtue of either accuracy gains or reduction in the number of grid points, which results in savings in computational cost. The same argument can be applied to the issue of temporal accuracy: on the one hand, there exist physical problems that require very accurate approximations in time to perform statistics. On the other hand, it is of importance to maintain a certain level of precision when advancing in time with the “largest possible” time step. Since, in general, the computational cost of a first-order method is hardly any lower than for higher-order methods, the

52

same level of temporal accuracy can be much faster obtained by a high-order method. In general, one is satisfied with third- or fourth-order accuracy in time, such as third- or fourth-order RungeKutta methods. Another issue is, often, the time step is not determined by considerations of accuracy, but by a stability condition, which will be discussed in later sections. Hence, we introduce a recently developed high-order method, optimized dispersion-relationpreserving (DRP) scheme to discretize the system of equations [37]. DRP schemes are easy to implement, its accuracy is adjustable from second-order to higher-order, it can be parallelized, and it incorporates the ability of adjustment of dispersion errors (phase errors in wave propagations). When diffusion is small at the early phase of quench propagation, DRP can effectively decrease the numerical errors (mainly phase errors in wave propagation problems).

4.1

Mathematical model

As discussed above, we will concentrate on a one-dimensional model, in which compressible liquid helium flow is only taken into account along the longitudinal direction of the channel in the CICC. The configuration used in this 1-D model is shown in Figure 4.1. The cable (strands) can be assumed to exchange heat by convection to the flowing helium, the heat diffusion along the length of the cable, parallel to the helium flow, can be retained in the heat balance of the conductor. The temperature gradients across the conductor are assumed negligible. There is temperature difference between each two components of the system: liquid helium, strands and conduit. The coupling effects between adjacent CICC layers or turns are also neglected, so the adiabatic boundary condition is assumed along the side of the channel of the cable. The inlet and outlet are imposed Dirichlet or Neumann boundary conditions. To initiate the quench in CICC, an external heat source is placed in the mid-point location, and a Gaussian profile of the heat pulse is assumed. During a quench, Joule heating occurs only in the strands. Because the conduit (stainless steel) has a high electrical resistivity and low thermal conductivity, we make the assumption that Joule heating is negligible in the conduit (i.e., current does not leave the strands during quench). It is important to note that helium flow is under the turbulent regime during quench propagation, and we circumvent solving the Navier-Stokes equation by modeling the viscous force with a wall friction, through the Darcy friction factor as an empirical function of

53

the Reynolds number.

Figure 4.1: Schematic of one-dimensional CICC, all strands are combined as a component, and there is a thermal coupling between each two components of conduit (jacket), strands, and liquid helium (fluid).

The basic governing equations consist of continuity, momentum and energy for the helium, the dependent variables being density, mass flux, and total energy per unit mass (stagnation enthalpy), we can write down the system in the conservative form:

ρ ρv      0   ∂     ∂  2 ρv ρv + p  +  = F  ∂t  ∂ x 2 2      1 1 ρ u + 2v   ρ v u + 2 v + p v   Qh  where f ρv v F =− 2dh Qh = β He, st (Tst − T ) + β He, jk (T jk − T )

(

)

(

)

54

(4.1.1)

(4.1.2) (4.1.3)

β He, st and β He , jk are the equivalent heat transfer coefficients: β He,st =

pst hHe, st

, β He, jk =

p jk hHe, jk

AHe AHe The equation of state of helium is added to close the system:

p = p( ρ , T )

(4.1.4)

(4.1.5)

All the variables without subscripts or superscripts refer to liquid helium. Besides the governing equations of liquid helium flow, the heat balance for the conductor (strands) is given by: ρ st Cst

∂ Tst ∂  ∂ Tst =  k st ∂t ∂ x  ∂x

h P  hst Pst  − (Tst − T ) − st − jk st Tst − Tjk + S(x, t, Tst ) A st Ast 

(

)

 I(t )   ρe S(x , t , Tst ) = J ρ e =   A Cu 

(4.1.6)

2

2

(4.1.7)

Where S ( x, t , Tst ) is the external heat term, such as Joule heating. Similarly, the energy equation for the conduit is given by: ρ jk C jk

∂ Tjk ∂t

=

∂ Tjk ∂   k jk  ∂x  ∂x

hst − jkPst  h jk Pjk − T jk − T − Tjk − Tst  A jk A jk 

(

)

(

)

(4.1.8)

The heat transfer coefficients used above are very important parameters for modeling the heat transfer between strands/conduit and helium. Some typical heat transfer coefficients correlations can be found in [Bloe86], [52], and [51]. For internal channel flow the Darcy friction factor, f, is given by: f l = 64 / Re (laminar)

(4.1.9)

f t = 0.184 / Re 0.2 (turbulent)

(4.1.10)

We do not include source term (Joule heat) in the right hand side of Eqn. (4.1.8). In fact the external heating sources are usually located in the strands and the Joule heating in the jacket is negligible due to its high electrical resistivity at cryogenic temperatures. Eqns. (4.1.1), (4.1.6), and (4.1.8), consist of the governing system of quench propagation.

4.2

High order finite difference schemes

We introduce a special high order method, spatial approximation with fourth-order to sixthorder, and temporal fourth-order scheme, the Dispersion-Relation-Preserving (DRP) scheme [37] to simulate quench propagation. First, we can discretize the spatial derivatives of the PDEs to

55

obtain a system of ordinary differential equations (ODEs), and then solve the ODE system obtained by use of a time integrator solver, which is similar to the method of lines (MOL) [62], [63]. To implement this high order scheme in Eqns. (4.1.1), (4.1.6), and (4.1.8) we show it by applying the high-order central difference to the first derivative and the second derivatives, viz. (n )

 ∂v 1 ( n)  v  = − vl ∆x  ∂ x l

m

∑a

j= − m

j

vl(+n)j

(n )

m  ∂  ∂ T     k = a j  k l(n+ )j    ∑  j= − m   ∂ x  ∂ x  l

(4.2.1) m

∑a

i =− m

j

T

( n) j + i +l

  

(4.2.2)

Because we apply central difference high-order scheme to spatial first derivatives, artificial (numerical) viscosities are necessary to make the calculation stable. After adding the artificial damping term to the right hand side of the DAE/ODEs, we have a form as follows: r r U t l = Kl

( )

(4.2.3)

Then using fourth-order time integration solver to (4.2.3) to obtain: 3 r r r U l( n+1) = U l(n ) + ∆t ∑ b j Kl(n− j )

(4.2.4)

j =0

where a j are the coefficients of the spatial difference stencil, b j are the coefficients of the time difference stencil. The basic idea of DRP method is that we can adjust these stencil coefficients a j and b j in wavenumber space and frequency space to get optimized effects on special

problems. The objective is to reduce the dispersion errors for wave propagation problems in the finite difference procedure [37], [64]. It is important to note that classic high order finite difference schemes are often constructed by the truncated Taylor series method, and do not optimize the stencil to damp short waves components (high frequency components) which is crucial for long time integration. The details of optimizing the spatial discretization in wavenumber space and the temporal discretization in frequency space are given in [37]. For simplicity, considering a uniform grid with index l in space and n in time, we have: N  ∂v  1  N    ≅  ∑ a−n vl− n + a 0v0 + ∑ an vl +n  n =1   ∂ x  l ∆ x  n=1

(4.2.5)

Optimized 4th order: a 1 = −a −1 = 0.6794 , a 2 = −a − 2 = −0.0897 , a 0 = 0

(4.2.6)

Optimized 6th order: a 1 = −a −1 = 0.7709 , a 2 = −a −2 = −0.1667 , a 3 = −a −3 = 0.02084

56

For comparison, the standard central finite differences discretization to first derivative are given by [36]:

 ∂v  1 (vl+1 − vl −1 ) + O ∆x 2 2nd order:   =  ∂ x  l 2∆ x

( )

∂v 1  = (−v l+ 2 + 8 v l+1 − 8 v l−1 + v l−2 ) + O ∆x 4 4th order:   ∂ x  l 12∆ x

( )

 ∂v  v − 9 vl+ 2 + 45 vl +1 − 45vl −1 + 9vl −2 − vl −3 6th order:   = l +3 + O ∆x 6 60∆x  ∂ x l

( )

Similarly we can find the stencil coefficients of the fourth-order time discretization scheme by optimizing the parameter in the standard stencil coefficients in frequency space through Laplace transform: b 0 = 2.3026 , b 1= −2.4910 , b 2 = 1.5743 , b 3 = − 0.3859

(4.2.7)

Because higher order schemes introduce additional numerical extrema, the computation is subject to loss of stability. For example, the fourth-order time integrator brings in three spurious roots to the system, so we must make sure the system is stable. We will discuss this issue later in Section 4.4. The problem of quench propagation is different from acoustic traveling waves problems in which we can impose outgoing, absorbing, and radiation boundary conditions [65]. However, we may take into account some practical conditions at the two ends of CICC channels, this is because the governing equations are nonlinear, and there is only an empirical relation for the equation of state for helium. Hence, we compare two types of physical boundary conditions for the one-dimensional model. One is the pressure expansion condition of the outlet [10], [59], which can be written as: p = p0 + α 12 ρ v 2

(4.2.8)

where α is the valve expansion loss coefficient. Another choice is to impose the Neumann condition [21]: ∂v =0 ∂x

(4.2.9)

After numerical experiments using these two physical boundary conditions, we found that Eqn. (4.2.9) is more stable than (4.2.8), so we impose it at the outlet to close the system of

57

equations. Because the helium is expelled from both ends, we do not consider inlet boundary conditions. Symmetrical boundary conditions are imposed in the middle point: ∂p =0, v=0 ∂x

(4.2.10)

Because there is not sufficient number of grid points to implement the high-order scheme, a special numerical process is necessary near the boundary, such as lower-order scheme, or an upwind method, or both. The most conventional method of boundary point treatment is to apply lower-order schemes near the boundaries. The typical replacement is first-order or second-order upwind approximation is given by [36], [66]: 1st upwind approximation in the out let:

 ∂v  1   = ( vl − vl −1 ) + O(∆x ) , l = N ∂ x ∆ x  l

(4.2.11)

2nd order upwind approximation in the outlet:

 ∂v  1   = ( vl− 2 − 4 vl −1 + 3 vl ) + O ∆x 2 , l = N ∂ x 2 ∆ x  l

( )

(4.2.12)

Even though the high-order scheme is degraded near the boundary, it is beneficial to the numerical stability, and it was found from many numerical experiments that the results with lower-order scheme near the boundary are far superior to those obtained by the lower-order scheme everywhere in the computation domain. However, the overall accuracy decreases compared to what is expected from a high-order scheme. If the accuracy is sufficient, a highorder scheme is beneficial from the point of view of computing efficiency.

4.3

Adaptive mesh technique

The solution of quench propagation in CICC exhibits steep fronts that move in time over the spatial domain. Because of large temperature and pressure gradients, the moving fronts need a special mesh process to obtain reasonable accuracy. Our goal is to resolve details around the moving front, the spatial discretization grid must be very fine around the interface to increase resolution of the solutions but can be significantly coarser far from the interface. There are two types of techniques to resolve large gradients at the propagation front: Moving Mesh Method

58

[67] and Adaptive Mesh Refinement (AMR) [68], [69]. Since the high-order scheme described above is optimized in unique grid space, we applied the AMR technique instead of moving mesh method. Another reason to introduce AMR is because our time integration is also high order and explicit, so the overly restrictive small time step of the finest grid in steep gradient regions is not applied over the entire computation domain. The basic AMR algorithm [68] is a cumbersome and complex method, and difficult to manage with a complicated hierarchical data structure. When the quench front propagates we can clearly define the moving interface at each time step. Therefore, the adaptive approach can be simplified using a user-controlled adaptive technique. We use the idea in [17] and make a modification in the fine grid generation. Since we can locate the position of the quench propagation front in one coarse cell, we can choose a simple mesh density in this small subdomain to obtain a finer grid. Instead of setting a Gaussian profile of mesh density distribution [17], we rather simply refine the mesh size by a factor of 2, and the level of fine grid is controlled by user input. A special artificial damping is added on the interface to make the numerical solution stable. The damping stencil is variable. Suppose the finer grid is generated, we may use the same highorder scheme like a uniform grid (single domain) as before. For the points nearby the finer/courser interface, the missing values of some stencil are found by interpolation, Figure 4.2 shows that three cases in which two need interpolation to complement their stencil.

59

Figure 4.2: Three cases (a)-(c) are shown on how to choose special stencils to compute first derivatives at node i near the interface of two different mesh-size sub domains. The cross indicates a mesh point included in the stencil.

For time marching, if single step method is used, the step size is restricted by the size of the smallest mesh. The coefficients of spatial derivatives, time integration, and the artificial damping are all adjusted to attain optimization (error tolerance). Front-tracking method is used to calculate the moving normal zone front [55]. A detailed usage of this method to track the superfluid helium to normal helium interface is given in [70].

4.4

Numerical stability analysis

It is important to note that it is not sufficient to have a good approximation of the spatial or time derivatives alone to simulate unsteady wave propagation accurately. Besides the fourthorder to sixth-order spatial discretization, an explicit fourth-order time marching approach is used in this paper, all these make numerical stability an important issue. We applied the

60

optimized time discreztization [37], which is constructed by applying the Laplace transform to the difference of the time derivative, and finding the minimum integrated error in time between the finite difference and the ODE. Because the first derivative in time is discretized by a fourthorder approximation, the finite difference marching scheme will contain spurious numerical solutions. We demonstrate the stability margin for optimized time marching schemes based on an ordinary differential equation, and will draw the stability regions in the complex ω plane: du = F (u ) dt

(4.4.1)

For simplicity we discretize scalar equation (4.4.1) to a two-level expression by weight average of the derivatives, u

n +1

 du  ≈ u + ∆t ∑ b j    dt  j= 0 1

(n− j )

n

n −1   du  n  du   ≈ u + ∆t b0   + b1     dt   dt   n

(4.4.2)

Using the Taylor series expansion of the right hand side of Eqn. (4.4.2) to get second order truncation error, suppose u is a continuous function,  du     dt 

n −1

n

n

n

 d2u   d 3u   du  =   − ∆t  2  + ∆t 2  3  + L  dt   dt   dt 

So b 0 + b1 = 1

(4.4.3)

− b1 = 0.5

(4.4.4)

In a two- level scheme we can only obtain second-order accuracy in time, hence this scheme is already optimized. Then we apply the Laplace transform to (4.4.2), and assume that u is replaced by continuous variables. Hence:

~  1  1  i ω j ∆t  du ~   ∑ b je iω j ∆t  ω ~ e − 1 u ≈ ∆t  ∑ b je = − i ∆ t (4.4.5)    u  j=0  dt  j= 0  or − i ω ~ u ≈ −i ω ~ u , where ω is the frequency of the differential equation, and we can introduce

(

)

−i ω ∆t

a numerical frequency,

( (

)

i e− iω ∆t − 1 ω = ∆t b0 + b1 ei ω j ∆t

)

(4.4.6)

Even if ω ∆t is real, ω ∆t may be a complex number. To find the relationship between ω ∆t and ω ∆t , change (4.4.6) to a high degree algebraic equatio n as,

61

(

)

2

iω ∆t

)

i ω ∆t

b1 e

i  iω ∆t i  +  b0 + − =0 e ω ∆t  ω ∆t 

(4.4.7)

or

(

0 .5 e

2

i  i ω ∆t i  − 1.5 + + =0 e ω ∆t  ω ∆t 

(4.4.8)

For a given ω ∆t this equation has two roots, one is the physical root, the other is the spurious root. The solution can be written as

e iω ∆t = 1.5 +

i i 1 ± 2.25 + − ω ∆t ω ∆t (ω ∆t ) 2

(4.4.9)

Figure 4.3 shows the two roots ω ∆t in complex plane when ω ∆t is a real variable, from which we get the stability margin is ω ∆t ≈ 0.67 . Similarly Figure 4.4 shows the two roots ω ∆t in the complex plane in which ω ∆t is a pure imaginary variable, from which we obtain that the stability margin is ω ∆t ≈ −0.70 . Figure 4.5 shows the stability region for ω ∆t for the two- level and four- level time marching schemes in the complex plane.

Figure 4.3: The real parts and imaginary parts of the two roots of ω∆t as a function of ϖ∆t in the complex plane (ϖ∆t is a real variable).

62

Figure 4.4: The real parts and imaginary parts of the two roots of ω∆t as a function of ϖ∆t in the complex plane (ϖ∆t is a pure imaginary variable).

Figure 4.5: The stability region of the second-order and fourth-order time marching schemes in complex plane of ω ∆t.

63

Based on the stability region of time marching schemes demonstrated in Figure 4.5, we can find out the maximum time step for explicit time integration. Instead of finding eigenvalues λ of the governing equations, we resort to the method in [71] to find the dispersion relation between frequency and wavenmber, then solve out the maximum time step. The common way is to deal with this problem is to linearize system of equations (4.1.1) by dropping the second-order terms and the source terms. We can introduce special units and recast the pressure equation and energy equation as the non-dimensional 1D linearized Euler equations (pp. 132-135 in [2]): ∂u ∂ p ∂ p ∂u + = 0, + =0 ∂t ∂x ∂t ∂ x

(4.4.10)

Suppose that U = (u p )T , Ε = ( p u )T . For simplicity the second-order central difference in spatial derivatives and the two-level discretization in time are used and we obtain the discretized form of Eqns. (4.4.10): 1

Uln +1 ≈ Unl + ∆t ∑ b jΚ (ln − j)

(4.4.11)

Ε nl+1 − Εln−1 2∆x

(4.4.12)

j= 0

K ln ≈ −

Apply Fourier-Laplace transform to (4.4.11) and (4.4.12), which result in frequency space and wavenumber space relation, the dispersion relation [71]:

 e i ∆xα − e − i ∆xα ~ ~ Κ = −i α Ε ≈ − 2∆x  e

−i ω ∆t

~ U+

~ ~  Ε = −i α Ε 

(4.4.13)

 1 ~ i ~ ~ −i ω ∆t Uinit 1 − e ≈ U + ∆t  ∑ b je i jω ∆t  Κ 2π ω  j=0 

(

)

(4.4.14)

So that we have the dispersion relation: ω ~ ~ ~ − i ω U + iα Ε = U init 2π ω

where ω = i

e −i ω ∆t − 1 1

∆t ∑ bj ei j ω ∆t

(4.4.15)

 e i ∆xα − e − i ∆xα , α = −i  2∆x 

 ~  , and the initial condition Uinit by the Laplace 

j =0

transform. Eqn. (4.4.15) can be simplified as ω = ±α . From Figure 4.5 we can get the following inequality

64

ω ∆t = α ∆t = (α ∆x )

∆t ≤ 0.67 ∆x

(4.4.16)

Figure 4.6 shows that the margin for α ∆x in wave number plane is below 1.57(≈π/2) for second-order central difference scheme. We can substitute it into (4.4.17) to get the maximum time step ∆t max : ∆t max ≤

0.67 ∆ x 0.67 ∆ x ≈ ≈ 0.42∆ x (α ∆ x )max π /2

(4.4.18)

Because the stability of fourth-order explicit method involves the interplay of spatial discretiztion of nonlinear terms, an analytical Courant-Friedrichs-Lewy (CFL)-like stability criterion is not readily available [72]. However, based on the Fourier-Laplace transform above we can determine the stability margin by multiplying by a constant the value obtained from analyzing the linearized system of equations.

Figure 4.6: The numerical dispersion errors as indicated by the dependence of the numerical wave number α on the wave number α of the partial differential equation, the straight line is the ideal non-dispersion error case.

65

4.5

Numerical damping and testing

It is well known that high-order schemes for numerical integration suffer from oscillation, particularly near steep gradients. The central difference scheme presented before is nondissipative and subject to grid-to-grid oscillation, we need to add artificial viscosities to make the calculation stable. The first introduction of artificial (numerical) viscosity is due to Neumann and Richtmyer [36]. It is composed of linear and quadratic terms, in which the linear term is used to damp background (numerical) noise, however, the quadratic term is applied to compress the oscillation in the region of a sharp discontinuities, e.g., a shock, where the velocities would go through large gradient. Another choice is Jameson’s numerical viscosity [36]. In this paper we introduce the numerical viscosity in a different way, and will test its validity later in this section through a simple wave equation. Artificial damping for numerical noise (grid-to-grid oscillation):

Di = −

c R∆ ∆x

m

∑d U( ) j

j =− m

n i+ j

(4.5.1)

Artificial damping for discontinuities (variable selective damping):

Di = −

vstencil Rstencil ∆x

m

∑d U( ) j

j = −m

n i+ j

(4.5.2)

where R∆ and Rstencil in Eqns. (4.5.1) and (4.5.2) are the numerical Reynolds numbers.

v stencil = v max − vmin

in a stencil. We note that these definitions of artificial viscosities are

considerably more general than those given by Neumann and Richtmyer or Jameson’s formulae. We also note that the latter can be placed into the above formats by performing some algebraic manipulation. We can choose a reasonable numerical Reynolds number to make the time step as large as possible. The correct choice of numerical damping coefficients and numerical Reynolds numbers is not a trivial problem. If the numerical damping is too large, the real physical dissipation will be smeared, and the numerical results will finally deviate from the correct solution. For large external disturbance, such as Joule heat induced by quench, variable coefficient damping, expressed by (4.5.2) in a stencil, is preferred.

66

The length of the damping stencil in the computing domain is subject to change. For example, no numerical viscosity is needed at the two ends, and three points damping stencil for the penultimate points, and five points damping stencil for the third grid near the two ends, and five or seven points stencil for other interior points of the domain. It is important to note that the formulation of numerical viscosity for a non- uniform grid is also different from that of a uniform grid. The difficulty is that the numerical viscosity is problem-dependent, the values chosen for quench simulation are different from those used for acoustic waves. In the rest of this section two test cases will be used to illustrate the usage of numerical damping for this high-order scheme. The first problem is a typical parabolic equation with initial and boundary values, the second is a step initial signal traveling problem, and the difference between numerical results and exact solutions is shown in detail. First, we apply the method to a simple diffusion equation as follows:

 ∂ 2 u ∂ 2u  ∂u  = κ  2 + ∂t ∂ y 2  ∂x

t =0 u =e

− ( ln 2 )

(x

2

+y 9

(4.5.3) 2

) (4.5.4)

The exact solution of the initial value problem for this diffusion equation is: −

u (x , y , t ) =

1 e 4(ln 2 )κ t 1+ 9

(

)

− x 2 + y2 9 +4 κ t (ln 2)

(4.5.5)

Let ∆x = ∆y be the length scale and (∆x )2 / κ be the time scale. After choosing t =

(∆x ) 2 t′ κ

and ∆x = ∆y, x = x ′∆x, y = y ′∆y , we drop the prime and we obtain the nondimensional diffusion equation and solve it in a finite domain [− 50,50 ]× [− 50,50] . The boundary conditions are given to be consistent with the initial condition. We compare the numerical solutions with the exact solution by plotting u(x,0,t) at t=40, 80. The comparisons are shown in Figures 4.7, which shows a good agreement between computation and the analytic solution.

67

Figure 4.7: Comparison of numerical solution of Eqns. (4.5.3)-(4.5.4) by 4th order scheme with exact solution at t=40, 80. Dirichlet boundary conditions are assumed at the ends.

The second example is a simple wave equation with a constant speed, ∂u ∂ u + =0 ∂t ∂x

(4.5.6)

with initial condition, t=0 u = H(x + 50 ) − H( x − 50) , where H(x) is the unit step function. Considering the initial conditions and the starting of the scheme, we have, u (l0 ) = H (x l + 50) − H (x l − 50 )

(4.5.7)

u l(n ) = 0.0

(4.5.8)

n 0  ∂x 

(A.2.1)

where α is a user-chosen parameter. Unfortunately this function is not analytical integrable and unsuitable for multi-phase problems, hence, the equidistribution principal has to be approximated using quadrature. So we consider equidistributing the monitor function below [89]

M (x ) = 1 +

µ1 µ

2 2

(x − s)

2

+1

.

(A.2.2)

121

where the parameters µ1 and µ 2 are positive constants to control the smoothness and clustering of the grid around the front s(t). The function (A.2.2) is analytic integrable and the grid points

{x } is the solution of the following scalar nonlinear equations [89] j

xj +

µ1 µ  j  sinh −1 (µ 2 (x j − s ) ) − 1  1 −  sinh −1 (µ 2 ( x L − s ) ) µ2 µ2  N 

j  µ  (x R − x L ) + 1 sinh N µ2

−1



(µ 2 ( x R − s ) ) = 0

(A.2.3)



for j=1,2,…,N-1. Figures A.3 and A.4 show the moving grid points track the He II/He I phase change front.

Figure A.3: A typical moving grid trajectories for the algorithm used with nondimensional scales. µ1 = 10 , and µ 2 = 1000 are used to control the clustering and smoothness of the grid points distribution along the He II/He I phase change front.

122

Figure A.4: A typical moving grid trajectories for the algorithm used in [70]. µ1 = 50 , and µ 2 = 1000 are used to control the clustering and smoothness of the grid points distribution along the He II/He I phase change front. More grid points are concentrating near the moving interface.

A.3 Adaptive mesh refinement (AMR) and other methods

Local grid refinement is different from the global rezone methods. Given a base grid, an error tolerance, and an error indicator (or monitor) for the grid cells, local refinement can be done as follows: First, compute the error estimate for each cell by the error indicator; second, flag those cells whose error estimate is larger than the error tolerance; finally, refine the flagged cell into several small ones, and group the newly refined cells into a new base grid. Go back to step 1 until the error tolerance is satisfied in each cell or the maximum refinement level is reached.

123

The local refinement is just adding new nodes to the base grid level by level until enough accuracy is achieved in each cell. This type technique is widely used in finite element method, which is called h-method. AMR is very simple and useful approach. In this thesis, AMR and moving grid method will be used to improve the accuracy near the steep gradient region of the solutions, and at the same time reduce the grid points needed in the uniform grid system. The choice of the other techniques is still an open topic.

124

APPENDIX B NOMENCLATURE

A B c C Dh E e f I h H J K L M Nu Re Pr Pe P p Wb q t T U v V ρ ρe S s

Area Applied magnetic field Speed of sound Specific heat Hydraulic diameter Total energy per unit volume Total energy per unit mass Friction factor Current Heat transfer coefficient Helium enthalpy Current density Thermal conductivity Length Mass Nusselt number Reynolds number Prandtl number Peclet number Wetted perimeter Pressure of helium Weber number Heat flux Time Temperature Internal energy per unit volume Velocity of helium Voltage Density Copper resistivity Heat source term Helium entropy

125

Subscripts He

Helium

Cu

Copper

St

Strands

Jk

Jacket (conduit)

In

Insulation

K

Kapitza conductance

s

Steady state

t

Transient state

Eq

Equivalent

He-st

Between helium and strands

He-jk

Between helium and jacket

Jk-st

Between jacket and strands

s/c

Superconductor

b

Bath values

w

wall

126

REFERENCES

[1] Seeber, B., Handbook of Applied Superconductivity, IOP Publishing Ltd., 1998. [2] Dresner, L., Stability of superconductors, Plenum Press, New York, 1995. [3] Wilson, M.N., Superconducting Magnets, Clarendon Press, Oxford, 1983. [4] Dresner, L., Propagation of normal zones in composite superconductors, Cryogenics 16, pp 675-681, 1976. [5] Dresner, L., Quench detection by fluid dynamic means in cable- in-conduit superconductors, Advances in Cryogenic Engineering 33, pp 167-174, 1988. [6] Dresner, L., Thermal expulsion of helium from a quench cable-in-conduit conductor, Proceeding 9th Symposium on Engineering Problems of Fusion Research, pp 618-621, Chicago, 1981. [7] Shajii, A., Freidberg J.P., Quench in superconducting magne ts II. Analytic solution, Journal of Applied Physics 76, pp 3159-3171, 1995. [8] Dresner, L., The growth of normal zones in cable- in-conduit superconductors, 10th Symposium on Engineering Problems of Fusion Research, pp 2040-2043, 1983. [9] Dresner, L., Protection considerations for force-cooled superconductors, 11th Symposium on Engineering Problems of Fusion Research, pp 1218-1222, Texas, 1985. [10] Arp, V. Stability and thermal quenches in force-cooled superconducting cables, Proceedings of Superconducting MHD Magnet Design Conference, MIT, Cambridge, pp. 142-157, 1980. [11] Marinucci, C., A numerical model for the analysis of stability and quench characteristics of forced- flow cooled superconductors, Cryogenics 23, pp 579-586, 1983. [12] Wong, R. L., Program CICC, flow and heat transfer in cable- in-conduit conductors equations and verification, Lawrence Livermore National Laboratory internal report, UCID-21733, 1989 (unpublished).

127

[13] Luongo, C.A., Loyd R.J, Chen F.K., and Peck S.D., Thermal- hydraulic simulation of helium expulsion from a cable- in-conduit conductor, IEEE Transaction on Magnets 25, no. 2, pp 1589-1595, 1989. [14] Bottura, L., Quench analysis of superconducting magnets. A numerical study, PhD thesis, University of Wales, 1991. [15] Bottura, L., and Zienkiewicz O.C., Quench analysis of large superconducting magnets. Part I: model description, Cryogenics 32, no. 7, pp 659-667, 1992. [16] Bottura, L., and zienkiewicz, Quench analysis of large superconducting magnets. Part II: model validation and application, Cryogenics 32, no. 8, pp 719-728, 1992. [17] Bottura, L., A numerical model for the simulation of quench in the ITER magnets, Journal of Computational Physics 125, pp 26-41, 1996. [18] Shajii, A., Freidberg J.P., Quench in superconducting magnets I. Model and numerical implementation, Journal of Applied Physics 76, pp 3149-3158, 1994. [19] Cheng, X., Numerical analysis of thermally induced transients in forced flow of supercritical helium, Cryogencis 34, no.3, pp 195-201, 1994. [20] Ito, T., Takata Y., and Kasao D., Stability analysis of a forced- flow cooled superconductor, Cryogenics 32, no. 5, pp 439-444, 1992. [21] Koizumi, N., Takahashi Y., and Tsuji H., Numerical model using an implicit finite difference algorithm for stability simulation of a cable- in-conduit superconductor, Cryogenics 36, no.9, pp 649-659, 1996. [22] Koizumi, N. et al, Quasi-two-dimensional numerical model for stability simulation of a cable- in-conduit conductor, Cryogenics 39, pp 495-507, 1999. [23] Wang, Q.L. et al, thermohydraulic simulation on CIC conductor with adaptive mesh finite volume method for KSTAR Tokamak superconducting magnet, IEEE Transactions on Applied Superconductivity 11(1), pp 2070-2073, 2001. [24] Wang, Q.L et al, Simulation on helium expansion and quench in CICC for KSTAR with moving mesh methods, IEEE Transactions on Applied Superconductivity 12(1), pp 15911594, 2002. [25] Miller, J.R., An overview and status of the NHMFL 45-T hybrid project, Teion Kogaku (Japanese Cryogenic engineering), Vol. 31, no 5, pp 240-249, 1996. [26] Miller, J.R., The NHMFL 45-T hybrid magnet system: past, present, and future, IEEE Transactions on Applied Superconductivity, 13 (2), pp. 1385-1390 Part 2, JUN 2003.

128

[27] Lax, P.D., Hyperbolic systems of conservation laws and the mathematical theory of shock waves, Region conference series in applied mathematics 11. Philadephia: SIAM, 1973. [28] Shajji, A., Freidberg, J.P., Theory of hybrid contact discontinuities, Physical Review Letters 76, no. 15, pp 2710-2713, 1996. [32] Osher, S., and Sethian J.A., Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations, Journal of Computational Physics 79, pp 12-21, 1988. [33] Ames, W.F., Numerical methods for partial differential equations, 3rd ed., Academic Press, 1992. [34] Zienkiewicz, O.C., and Taylor R.L., The finite element method. Volume 3: Fluid dynamics, 5th ed., 2000. [35] Canuto, C., Hussaini M.Y., Quarteroni A., and Zang T.A., Spectral Methods in Fluid dynamics, Springer-Verlag, Berlin, 1988. [36] Hirsch, C., Numerical computation of internal and external flows, Vol. 1, Wiley, New York, 1988. [37] Tam, Chris.K.W. et al, Dispersion-relation-preserving finite difference schemes for computational acoustics, Journal of Computational Physics 107, pp 262-281, 1993. [38] Patera. A.T., A spectral element method for fluid dynamics; laminar flow in a channel expansion, Journal of Computational Physics 54, pp 468-488, 1984. [39] Maday, Y., and Petera A.T., spectral element methods for the Navier-Stokes equations, in A.K. Noor, editor, State-of-the-art surveys in computatiobal mechanincs, ASME, New york, 1988. [40] Ho, Lee-wing, A Legendre spectral element method for simulation of incompressible unsteady viscous free-surface flows, PhD thesis, Massachusetts Institute of Technology, 1989. [41] Prestemon, S.O., Development of a spectral-element code for the study of heat and momentum transfer in turbulent flows through rough channels, PhD thesis, Florida State University, 2001. [42] Miller, J.R., Eyssa Y.M., Sayre S., and Luongo, C.A., Analysis of observations during operation of the NHMFL 45-T hybrid magnet system, , Cryogenics 43(3), pp. 141-152, MAR 2003.

129

[43] Mao, S., Luongo C.A., and Miller J.R., Investigation of index heating as source of quench in the NHMFL 45- T hybrid superconducting outsert, IEEE Transactions on Applied Superconductivity 12(1), pp 1520-1523, 2002. [44] Bottura, L., Transient thermal analysis of forced-flow condctors, Proceedings of the 12th ICEC Conference, pp 182-186, Southampton, 1988. [45] Iwasa, Y., Case studies in superconducting magnets: design and operational issues, Plenum Press, New York, 1994. [46] Minervini, J.V., and Bottura, L., Stability analysis of NET TF and PF conductors EC report, IEEE transactions on magnets 24, pp 1311-1314, 1988. [47] Arp, V., and McCarty R.D., Thermophysical properties of Helium-4 from 0.8 to 1500 K with pressures to 200 MPa, NITS Tech Note 1334, US Government Printing Office, Washington, USA, 1998 (revised). [48] McCarty, R.D., Thermophysical properties of Helium-4 from 2 to 1500 K with pressures to 1000 atmospheres, NBS Tech Note 631, Washington, USA, 1972. [49] White, F.M., Viscous fluid flow, 2 ed., MacGraw-Hall, 1991. [51] Giarratano, P.J., and Steward W.G., Transient forced-convection heat transfer to helium dump a step in heat-flux, Journal of Heat transfer-Transactions of the ASME 105(2), pp 350357, 1983. [52] Giarratano, P.J., Arp V.D, and Smith R.V., Forced convection heat transfer to supercritical helium, Cryogenics 11, no. 5, pp 385-393, 1971. [53] Krafft, G., Heat transfer below 10 K, in Cryogenic engineering, B.A. Hands ed., pp 171-192, Academic Press, 1986. [54] Anghel, A., etc., The quench experiment on long length (QUELL): Final Report, 1997. [55] Crank, J., Free and moving boundary problems, Clarendon Press, Oxford, 1984. [56] Miller, J.R., Walsh, R.P., Haslow, M.R., Kenny, W.J., and Miller, G.E., Characterization of high-current Nb3Sn cable-in-conduit conductors vs applied sheath strain, Adv. Cryo. Eng., vol. 44, pp. 967-974, 1998. [57] Bonito Oliva, A. et al., Calculation of AC losses and temperature increase in the superconducting coils during the ramp up of the hybrid magnet, NHMFL Report (unpublished) July 1993. [58] Frost, W., Heat transfer at low temperatures, Plenum Press, New York, 1976.

130

[59] Van Sciver, S.W., Helium Cryogenics, Plenum Press, New York, 1986. [60] Barron, R.F., Cryogenic Heat Transfer. Taylor & Francis, 1999. [61] Smith, G.D., Numerical solution of partial differential equations, 3rd ed., Oxford University Press, Oxford, 1985. [62] Gear, C.W. Numerical Initial Value Problems in Ordinary Differential Equations, PrenticeHall: Englewood Cliffs, New York, 1971. [64] Lele, S.K., Compact finite difference schemes with spectral- like resolution, Journal of Computational Physics, 103(1), pp. 16-42, 1992. [65] Engquist, B. and Majda, A. Radiation boundary conditions for acoustic and elastic wave calculations, Commun. Pur. Appl. Math., 32(3), pp. 313-357, 1979. [66] Warming, R.F. and Beam, R.M., Upwind second-order difference schemes and applications in aerodynamics flows, AIAA J., 14, 1241, 1976. [67] Dorfi, E.A., and Drury L. O’c., Simply adaptive grids for 1-D initial value problems, Journal of Computational Physics 69, pp 175-195, 1987. [68] Berger, M.J. and Oliger, J., Adaptive mesh refinement for hyperbolic partial differential equations, Journal of computational Physics, 53, pp. 484-512, 1984. [69] Li, Shengtai, Adaptive mesh methods and software for time-dependent partial differential equations, PhD thesis, University of Minnesota, 1998. [70] Mao, S., Luongo C.A., and Miller J.R., Analysis of the NHMFL 45-T hybrid magnet thermal behavior, Cryogenics 43 (3), pp. 153-163 MAR 2003. [71] Whitham, G.B., Linear and Nonlinear Waves, pp 161-165, Wiley, 1974. [72] Roache, P., Computational Fluid Dynamics, Hermosa Publishers, Albuquerque, 1972. [73] Ando, T., Propagation velocity of the normal zone in a cable- in-conduit conductor, Advances in Cryogenic Engineering 35A, pp 701-708, 1988. [74] Marinucci, C., Bottura, L., Vecsey, G., and Zanino, R., The QUELL experiment as a validation tool for the numerical code Gandalf, Cryogenics, 38(5), pp. 467-477, 1998. [75] Gottlieb, D., and Orszag S.A., Numerical Analysis of Spectral Methods: Theory and Applications, SIAM, Philadelphia, 1977. [76] Kopriva, D.A., Spectral Methods in PDEs, Course Note of MAD5745, Florida State University, Math Department, Fall 2002.

131

[77] Trefethen, L.N., Spectral Methods in MATLAB, SIAM, Philadelphia, 2000. [78] Arp, V., Private Communication, April 2004. [79] Roe, P.L., Approximate Riemann solvers, parameter vectors, and difference-schemes, Journal of Computational Physics, 43(2), pp. 357-372, 1981. [80] Press, W.H., Flannery, B.P., Teukolsy, S.A., and Vetterling, W.T., Numerical recipes, Cambridge University Pres, Cambridge, 1989. [81] Shu, C-W., Efficient implementation of essentially non-oscillatory shock-capturing schemes, II, Journal of Computational Physics, 83, pp. 32-78, 1989. [82] Fischer, P.F., Spectral Element solution of the Navier-Stokes Equations on High Performance Distributed-Memory Parallel Processors, PhD thesis, Massachusetts Institute of Technology, 1989. [83] Rønquist, E. M., Optimal Spectral Element Methods for the Unsteady Three-dimensional Incompressible Navier-Stokes Equations, PhD thesis, Massachusetts Institute of Technology, 1988. [84] Couzy, W., Spectral Element Discretization of the Unsteady Navier-Stokes Equations and Its Iterative Solution on Parallel Computers, PhD thesis, Ecole Poly-technique Fédérale de Lausanne, 1995. [85] Morton, K.W., Numerical Solution of Convection-Diffusion Problems. Chapman-Hall, London, 1996. [86] Hawken, D.F., Gottlieb J.J., and Hansen J.S., Review of some adaptive node- movement techniques in finite-element and finite-difference solutions of partial differential equations, J. Comput. Phys. 95, pp 254-302, 1991. [87] De Boor, C., Good Approximation by Splines with Variable Knots, II. Springer-Verlag, Berlin, 1973. [88] Hyman, J.M., and Naughton M.J., Static rezone methods for tensor-product grids, in Proceedings of SIAM-AMS conference on large scale computations in fluid mechanics, SIAM, Philadelphia, 1984. [89] Mackenzie, J.A., and Robertson, M.L. The numerical solution of one-dimensional phase change problems using an adaptive moving mesh method, Journal of Computational Physics 161, pp. 537-557, 2000. [90] Yaskin, L.A., Jones M.C., Yeroshenko V.M., Giarratano P.J., and Arp V.D., A correlation for heat transfer to supercritical helium in turbulent flow in small channels, Cryogenics 17, pp 549-552, 1977.

132

BIOGRAPHICAL SKETCH

Personal Shaolin Mao was born in Wuxie, a small city nearby the Yangtze River, Hubei Province, P. R. China. He got his Ph.D. degree in Mechanical Engineering from Florida State University in fall semester, 2004. He got his M.S. degree in Engineering Mechanics from Tsinghua University, P. R. China in July 1992. He also got his B.S. degree in Aerospace Engineering from National University of Defense Technology, P. R. China in July 1989. He has seven years of industrial experiences in central air conditioning system, air conditioner, control systems, video game software development, education software development, and electronic products design and development in south China from 1992 to 1999. His main interests cover development of large-scale engineering software, numerical simulation, scientific computation, fluid mechanics and heat transfer application, etc.

Publication •

Mao, S., Luongo, C.A., and Kopriva, D.A., Discontinuous Galerkin Spectral Element Simulation of quench Propagation in Superconducting Magnets, to appear in IEEE Trans. Appl. Supercond., 2004.

•

Mao, S. and Luongo, C.A., Numerical Simulation of Quench Propagation in Magnets Using High Order Schemes, in prepare, 2004.

•

Mao, S. and Luongo, C.A., Numerical Analysis of a Stefan Problem in Large Scale Superconducting Magnets System, in prepare, 2004.

•

Mao, S. and Luongo, C.A., Numerical Simulation of Moving Boundary Problems in Superconducting Magnets, Moving Boundary 7, pp. 359-369, 2003, WIT Press.

•

Mao, S., Luongo, C.A. and Miller, J.R., Numerical Simulation of Thermal Stability in NHMFL 45-T Hybrid Magnet Using Front-Track Method, The Proceeding of the Second M.I.T. Conference on Computational Fluid and Solid, pp. 1438-1441, 2003, Elsevier.

133

•

Mao, S., Luongo, C.A. and Marinucci, C., Quench simulation in Superconducting Cables Using Optimized DRP Scheme, IEEE Trans. Appl. Supercond., 13(2), pp. 1692-1695, 2003.

•

Mao, S., Luongo, C.A. and Miller, J.R., Thermal-Hydraulic Analysis of the NHMFL 45T Hybrid Outsert Coils, Cryogenics, 43 (3), pp. 153-163, 2003.

•

Mao, S., Luongo, C.A. and Miller, J.R., Investigation of Index Heating as Source of Quench in the NHMFL 45-T Hybrid Superconducting Outsert, IEEE Trans. Appl. Supercond., 12(1), pp. 1520-1523, 2002.

•

Mao, S., The Design of the Music Algorithm for Video Game System, in Lecture Notes in Scientific Computation, Edited by Luoluo Li, Zhongying Chen and Yongdong Zha ng, International Culture Publishing Inc., 167-172, Beijing. (ISBN 7-5026-1176-2) 2000.

134