Determining Stability Margins in Adiabatic Superconducting Magnets with 3-D Finite Element Analysis Author: Arnaldo Rodríguez González Department of Mechanical Engineering University of Puerto Rico at Mayagüez Supervisor: Tengming Shen, Ph.D. Superconducting Materials Department Fermi National Accelerator Laboratory SIST Program 2012 Table of Contents Abstract ................................................................................................................................................................ 2 Introduction ......................................................................................................................................................... 2 MATLAB prediction program & the MIITS method ........................................................................................... 4 FEA simulation ..................................................................................................................................................... 8 Pre-solving method ......................................................................................................................................... 8 Solving method .............................................................................................................................................. 10 Tested cases ....................................................................................................................................................... 12 Single-strand case ......................................................................................................................................... 12 Multi-strand case........................................................................................................................................... 13 Results ................................................................................................................................................................ 15 FEA vs. MIITS .................................................................................................................................................. 15 Minimum Quench Energy ............................................................................................................................. 17 Conclusions ........................................................................................................................................................ 18 Future use .......................................................................................................................................................... 19 Acknowledgements ........................................................................................................................................... 19 References ......................................................................................................................................................... 20 1 Abstract Superconducting magnets play a key role in the development of experiments at Fermilab; understanding the operating stability of these can allow us to utilize more potent magnets for future experiments (like the proposed Muon Collider), optimize the design of magnets in more immediate experiments (like Mu2e), and research the future use of more exotic materials (like high-temperature superconductors). This summer, I developed a 3-D parametric FEA program in ANSYS Mechanical APDL that simulates quench in superconducting magnets, and I also developed a parametric MATLAB program that predicts thermal behavior in magnet quench using the MIITS method. These programs can provide useful quenching parameters (like minimum quench energy and normal zone propagation velocity) for different cases of quench, leading to the previously mentioned objective of magnet design optimization. To test the programs, preliminary cases were run and the data produced was compared and analyzed. The results of these analyses, as well as the program operating methods, are discussed in this project. Introduction Superconductivity, discovered in 1911 by Heike Kamerlingh Onnes, is a quantum mechanical phenomenon where a material shows zero electrical resistivity and expulsion of internal magnetic fields below a characteristic critical temperature T. This effect is extremely useful for high-field magnet c applications: since no Joule heating in the electromagnet is generated regardless of operating current, very high-intensity fields can be generated due to high operating currents in the magnet itself. This application was pioneered here at Fermilab with the Tevatron, and the technology has branched out into various other fields, like MRI (Magnetic Resonace Imaging) machines and mass spectrometry. However, it does come with certain operating risks, the biggest possibly being quench. 2 Quenching in superconductors occurs when a section of the superconductor snaps out of the superconductive state and become normally resistive. Because of this, Joule heating is generated in this section, which then expands the normal zone and causes a chain reaction in which the entire conductor becomes resistive. This, if left, unchecked, can easily damage or destroy the magnets and any surrounding parts. A quench event can be initiated if any of the three following things occurs: 1) The temperature of any section of the magnet exceeds the critical temperature T . c 2) The current of any section of the magnet exceeds the critical current I c . 3) The magnetic field through any section of the magnet exceeds the critical field H . c Interestingly enough, these values are coupled: operating the superconducting magnet at a higher temperature, for example, yields a smaller I . This coupling causes behavior that not only takes into c consideration quench due to critical temperature, but also due to other parameters (like critical current) which change due to changes in operating conditions. In most composite magnets, this leads to current sharing. Current sharing is a phenomenon that occurs when a composite superconductor (generally described as superconducting filaments encased in a normally resistive matrix) surpasses the critical current for its operating temperature. Due to this, a certain portion of the current leaps out of the superconductor and flows through the matrix, creating Joule heating. Since the values of critical current, critical field and critical temperature are coupled, and current sharing is mediated by the critical current, it’s possible to determine a current sharing temperature T that specifies when current sharing behavior starts for a cs constant operating current, operating magnetic field, and initial temperature. The equation that determines this is below. 3 (cid:5).(cid:12)(cid:13) (cid:5).(cid:12)(cid:13) (cid:9) (cid:9) (cid:19)(cid:20) (cid:1)(cid:2)(cid:3) =(cid:1)(cid:2)(cid:5)(cid:6)1− (cid:10) − (cid:15)(cid:16)(cid:1)(cid:2)(cid:5)(cid:6)1− (cid:10) (cid:17)−(cid:1)(cid:18) (cid:9)(cid:2)(cid:5) (cid:9)(cid:2)(cid:5) (cid:19)(cid:2)(cid:5) (Eq. 1) This current sharing temperature depends on multiple values: it depends on the critical temperature at zero field , the critical field at zero degree , the critical current at zero field , the operating (cid:1)(cid:2)(cid:5) (cid:9)(cid:2)(cid:5) (cid:19)(cid:2)(cid:5) current , the operating field , and the initial operating temperature . This T variable is the variable cs (cid:19)(cid:20) (cid:9) (cid:1) directly used in the simulation to determine when current sharing behavior starts. Of note is the fact that the phenomenon is dependent on the temperature of the superconducting filament/s: but the heating itself is dependent on the temperature of the matrix, since the resistivity is temperature- dependent. MATLAB prediction program & the MIITS method To estimate the thermal behavior of the magnet for various cases in a quick fashion, a prediction program was developed in MATLAB that uses conditions similar to the simulated case. The prediction program uses numerical integration among other techniques to quickly approximate the behavior of the maximum hot spot temperature in the conductor as a function of time, centering the solution process around the MIITS method. It does the calculation based on the following equation: (cid:26)(cid:1) (cid:28)(cid:29)(cid:24)(cid:1)(cid:25) (cid:30) (cid:21)(cid:2)(cid:22)(cid:23)(cid:2)(cid:22)(cid:24)(cid:1)(cid:25) = (cid:19) (cid:20)(cid:31)(cid:24)t(cid:25) (cid:26)(cid:27) (cid:21)(cid:29) (Eq. 2) This equation describes adiabatic heating in a constant-current mode per unit length for a homogeneous conductor with a volumetric heat capacity that is a weighted average of the heat capacities that make it 4 up. For the following case, the composite heat capacity is determined by the following equation, where is the ratio of matrix material to superconductor: ! ! 1 (cid:23)(cid:2)(cid:22) = (cid:23)(cid:29)(cid:6) (cid:10)+ (cid:23)(cid:3)(cid:2)(cid:6) (cid:10) !+1 !+1 (Eq. 3) x 106Volumetric heat capacities for copper, niobium titanium, and weighted composite average 3.5 Copper 3m)) 3 Niobium Titanium K Weighted Average ( J / 2.5 ( y cit 2 a p a c at 1.5 e h etric 1 m olu 0.5 V 0 0 25 50 75 100 125 150 175 200 225 250 275 300 Temperature (K) Figure 1: Heat capacities for Cu, NbTi and weighted average with γ = 1.3 Inserting the term (Eq. 4) and understanding that (Eq. 5), the equation % #(cid:29) = (cid:19)(cid:20)(cid:31) ⁄(cid:21)(cid:29) (cid:21)(cid:29) ⁄(cid:21)(cid:2)(cid:22) = %&’ is transformed into: (cid:26)(cid:1) ! (cid:30) (cid:23)(cid:2)(cid:22)(cid:24)(cid:1)(cid:25) = (cid:6) (cid:10)(cid:28)(cid:29)(cid:24)(cid:1)(cid:25) # (cid:29) (cid:26)(cid:27) !+1 (Eq. 4) Rearranging the equation so that temperature-dependent terms are together and performing appropriate integration, the result is: 5 )* (cid:23)(cid:2)(cid:22)(cid:24)(cid:1)(cid:25) ! (cid:30) ( (cid:26)(cid:1)= (cid:6) (cid:10)# (cid:29)(cid:27) (cid:28)(cid:29)(cid:24)(cid:1)(cid:25) !+1 ) (Eq. 5) The right side of the equation has been considerably simplified by assuming that the initial time is zero and the current density in the matrix is constant over time; therefore, the MIITS method prediction does not take into account current sharing behavior at all. The left-hand side of this equation is very important: this parameter, known as the Z function where = (Eq. 6), describes a characteristic curve which shapes the thermal behavior of )*/01(cid:24))(cid:25) +(cid:24)(cid:1),,(cid:1)(cid:25) .) 23(cid:24))(cid:25)(cid:26)(cid:1) the conductor for any applied current. Thermal curves differing from one another for various cases with same conductor materials are only due to the scaling factor , which depends on the operating % (cid:30) 4%&’5# (cid:29) current and on the ratio of matrix to superconductor of the strand. x 1016 Z value for copper and weighted composite average as a function of temperature 14 13 Copper 12 Composite 11 10 4m) 9 s / 8 2A 7 e ( 6 u val 5 Z 4 3 2 1 0 0 25 50 75 100 125 150 175 200 225 250 275 300 Temperature (K) Figure 2: Z plots for copper and weighted composite average with γ = 1.3 6 This prediction program, like the simulation, is parametric: the program asks for key values as initial input and then performs the calculation based on these. This allows for multiple predictions to be made based on a changing single parameter, like operating current, and for their behavior to be compared. An example of predictions based on changing parameters is shown below, in this case for a changing reduced current (operating current to critical current ratio). Predicted maximum hot spot temperatures for different operating current modes vs. time 250 25% Io/Ic 225 50% Io/Ic 200 75% Io/Ic 175 K) e ( 150 ur at 125 er mp 100 e T 75 50 25 0 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 Time (s) Figure 3: MIITS predictions for varying operating currents However, the prediction program has its limitations: as previously mentioned, it does not take into account either current sharing behavior or thermal diffusion. Whether or not the current sharing phenomenon changes the behavior of the temperature drastically in comparison to the prediction is discussed in the results section. Also, the prediction program assumes the distribution of material in the conductor is uniform in its cross-section. This can cause deviations in the comparison between simulated temperature and predicted temperature due to an increasing temperature difference between the matrix and the superconductor as a result of differing heat capacities. 7 FEA simulation To be able to account for these behaviors and conditions, a parametric FEA simulation was developed in ANSYS Mechanical APDL 14 that simulates quench in superconducting magnets of variable geometries and under variable conditions. The equation that governs the quench phenomenon in the FEA simulation is a power density equation, defined as: (cid:26)(cid:1) (cid:30) (cid:23)(cid:24)(cid:1)(cid:25) = ∇∙(cid:24)8∇(cid:1)(cid:25)+# (cid:28)(cid:24)(cid:1)(cid:25)+ℎ:;:<:=> (cid:26)(cid:27) (Eq. 7) This equation incorporates both thermal conductivity and the initial heat pulse into the solution, unlike the MIITS program. It also allows for variations in the cross-section of the strand, since the equation describes power per unit volume instead of per unit length as in the MIITS case. In fact, one can derive Equation 4 (sans the ratio term, which wouldn’t be needed in this case) from Equation 7 by removing the thermal conduction and initial heat pulse terms. However, Equation 7 by itself cannot incorporate current sharing behavior, so that behavior had to be taken into account directly by “hard-wiring” it into the FEA program. This meant that a (relatively) simple coupled-field analysis with ANSYS would not do: the electrical phenomena in the superconductor had to be manually included into the program through macros in the solution step, while ANSYS took care of the more complicated thermal phenomena. The methods used to accomplish this, as well as other important program methods, are discussed below. Pre-solving method The first thing that occurs in the simulation program is the assignment of values to important parameters in the simulation: these include the geometry parameters, the quench section length (which determines the disturbance spectrum of the quench), the initial heat energy input, the critical 8 temperature and the operating current among others. Secondary values (like current sharing temperature) are calculated automatically based on these primary input variables. The pre-solving method is as parametric as possible: soft alterations to the geometry and to other numerical parameters in the simulation do not require any additional changes to the code. However, due to the use of spatial “differentials” in the solving process (which are explained in more detail in the solving method section), it is required that finite element nodes be located at regular intervals on the z- axis. This is achieved via performing line divisions on the longitudinal line geometries prior to the meshing. (An important condition worthy of note that this creates is that the geometry is then required to be located in Cartesian coordinates and to extend longitudinally over the z-axis. If a geometrically faithful representation of a solenoid were to be built, for example, it would require significant alteration to the program.) Figure 4: Example mesh for a given superconductor geometry. (Notice regular line intervals) 9