Numerical study on flux-jump occurrence in a cup-shaped MgB2 bulk for magnetic shielding applications

MgB2 is one of the most promising materials for superconducting bulk applications. However, thermomagnetic instabilities can arise in the material because of its low heat capacity and thermal conductivity as well as its high critical current density. Being able to predict these phenomena, can guide and optimize MgB2-based devices for magnetic flux shielding or trapping applications. In this work, the flux-jump occurrence in an MgB2 cup-shaped shield is numerically studied using the finite element method by means of the commercial software COMSOL 6.0 Multiphysics®. To this aim, we developed a 2D axial-symmetric model coupling the heat diffusion equation and the magnetic equations based on a magnetic vector-potential ( A⃗ ) formulation. The comparison of the computed shielding curves with the experimental ones evidenced a good agreement between the two sets of data at different temperatures and positions along the shield’s axis. The as-validated model was then exploited to investigate possible optimization routes via the improvement of both the thermal conductivity of the material and the thermal exchange between the device and the cooling stage.


Introduction
Superconducting (SC) materials are nowadays an innovative solution for a wide range of applications, mostly as wires, tapes, and thin films [1,2]. However, in the last few years, also bulk materials have found a more widespread use thanks to the * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. improvement of both their growth techniques and the cooling technologies required to exploit their potential [3]. SC bulks are high-performing candidates to replace the traditional solutions for applications such as permanent magnets or magnetic shields [4,5], the latter with good results in several fields ranging from biomedical engineering to high-sensitivity instrumentation protection [6][7][8].
In this framework, MgB 2 is one of the most promising SC materials for bulk applications [3,9], due to the low-cost and non-toxic precursors not including rare earth elements, the low weight density, and the long coherence length. In particular, the last characteristic implies that clean grain-boundaries do not prevent the flow of high current densities, thus opening the employment of large polycrystalline manufacts-fabricated by in situ or ex situ sintering processes [10][11][12][13][14] or by infiltration processes [9,15,16]-for large scale applications.
Although a lot of efforts have been made in order to fabricate more and more efficient shields, both improving the material properties and optimizing the manufact shape [17][18][19], the SC shield performances can still be weakened by different uncontrolled causes, such as quench phenomena or other thermomagnetic instabilities [20]. In the case of MgB 2 samples, despite the improvements in the fabrication processes, thermomagnetic instabilities such as flux jumps still occur at low operating temperatures. A flux-jump consists of a sudden penetration of the magnetic field inside the superconductor coupled with a rapid increase of its temperature, which can even induce the SC material to switch from SC to normal state. In general, the trigger of a flux jump is described as an iterative process [21,22]: (a) an initial amount of heat δQ 0 is released (for instance due a small variation of the external field or temperature with a resulting magnet flux penetration/motion in the sample), (b) as a consequence the initial SC temperature T 0 is locally increased by δT 0 ; (c) since the critical current density J c (T) is a decreasing function of temperature, there is a reduction of the density of the screening current; (d) due to the motion of magnetic flux into the SC, an electric field perturbation δE 0 is induced, (e) which, in turn causes an additional heat release δQ 1 and so on. The theoretical investigation of the flux-jump phenomena usually considers the accomplishment of a local adiabatic condition, which is assumed fulfilled when D t ≪ D m , being D t and D m the thermal and the magnetic diffusion coefficients, respectively [21,22]. For type II SCs the adiabatic condition D t ≪ D m are assumed to be satisfied [22] and in MgB 2 , this condition is verified at low temperatures due to the low heat capacity and thermal conductivity of this compound [14,[23][24][25][26], as estimated by Chabanenko et al [27].
Thermomagnetic instability-driven flux avalanches in MgB 2 have been revealed by several techniques both in films, tapes, and bulks [28][29][30][31][32][33][34]. It is well known that these instabilities of the critical state can nucleate, for instance, due to small temperature fluctuations or variations in the external magnetic field [35,36]. Their occurrence probability is affected by the ramp rate of the external applied magnetic field and by the working temperature [37,38]. Moreover, Aldica et al [39] showed that fabrication parameters, such as the heating temperature in spark plasma sintering (SPS) process, can influence the flux jump occurrence, as well. The flux-jump occurrence drastically deteriorates the intrinsic capacity of a superconductor to shield magnetic fields [14,40], becoming the major factor limiting the shielding performance at low temperatures [20].
In order to increase the critical state stability against fluxjumps, a deep investigation of the occurrence conditions as well as of the flux avalanche evolution is needed. To this aim, being able to know and model the in-field behavior of a SC shield is crucial to predict its properties and optimize its layout and thermal anchoring for the specific application. Since relying on suitable simulation tools can be a successful approach to deeply understand the intrinsic electromagneticthermal behavior of the SC materials [41], in the last few years, computational approaches were applied with the aim to reproduce and predict the flux jump phenomena [40,[42][43][44].
In our previous works [14], the shielding properties of a single-capped MgB 2 tube (hereafter named cup) were experimentally investigated in the range of operational temperatures (T OP ) from 20 to 35 K, highlighting the presence of flux jump phenomena up to 32.5 K. The sample was obtained by means of a SPS technique and the occurrence of this phenomenon up to so high temperature could be ascribed to the addition of hexagonal boron nitride (hBN) powder to the starting MgB 2 powder (MgB 2 + 10%wt hBN). This addition improves the SPS bulk machinability [45], but, as it happens for graphene addition [46], at the cost of a major flux jumps occurrence, which can be explained by a reduction of the thermal conductivity. However, the computational approach that we previously employed to reproduce the experimental data was focused only on the electromagnetic behavior of the superconductor, thus resulting unsuitable to predict the flux-jumps occurrence [47].
To overcome this shortage, in this paper, the numerical procedure based on the magnetic vector-potential ( ⃗ A) formulation described in [48] is coupled with a thermal model. The resulting approach can predict the decay of the shielding properties due to the flux-jump occurrence inside the cup-shaped MgB 2 bulk when an external field is applied parallel to the sample axis. Then, a comparison between the experimental and the computational data is carried out. Finally, we give a prediction of feasible enhancements in the shielding capability to be achieved by improving the MgB 2 thermal conductivity and the heat exchange with the cryocooler cooling stage. Using the numerical approach based on the ⃗ A-formulation instead of the widely used ⃗ H-formulation allowed us to compare the results presented in this work with those obtained in our previous works [13,14,47] and to take advantages of the less computational-time-consuming properties of the ⃗ Aformulation [19].
The paper is organized as follows. Section 2 presents the finite-element method used for the numerical calculations. In particular, sections 2.2 and 2.3 deal with the electromagnetic and the thermal equations, respectively. The results obtained by the coupled formulation described in the first part of the paper are listed in section 3, where experimental and computational data are compared for different working temperatures and predictions of attainable enhancements are given. The main outcomes are summarized in section 4.

Numerical modeling
To reproduce the effects of the magneto-thermal instabilities on the shielding properties of an MgB 2 cup-shaped bulk, we combined the magnetic vector potential formulation discussed in [48] with the heat balance equation [43,44]. The model was implemented by the finite element software COMSOL 6.0 Multiphysics® [49] through the magnetic field interface (mf ) coupled with the heat transfer module (ht). To validate this numerical approach by comparing experimental and computed data, we assumed that the sample here investigated by numerical modeling has the same size and physical properties as the SC cup-shaped shield whose experimental characterization is reported in [14]. A sketch of the sample is drawn in (figure 1(a)). As in [14], the sample was assumed subjected to the external magnetic field ⃗ H app (t) applied parallel to the cup's axis (z-axis, henceforth) after the zero-field cooling. Taking advantage of the cylindrical symmetry of the sample, we used an axisymmetric 2D approach, establishing a cylindrical coordinate system (r, ϕ, z) with the origin placed on the base center of the sample and the z-axis parallel to the applied field direction.

Boundary conditions and modeling constraints
To calculate the effect of the thermomagnetic instability on the local magnetic field inside the shield, the electromagnetic equations (section 2.2) were solved for both the SC (Ω SC ) and the external domain (henceforth named vacuum domain, Ω vacuum ) ( figure 1(b)). In contrast, the thermal equation (section 2.3) was solved only for the Ω SC domain. The value of the applied field as well as the thermal exchange between the superconductor and the cooling stage was accounted for by means of the boundary conditions. In more detail, at a large distance from the MgB 2 shield, the magnetic flux density ⃗ B was assumed to be equal to µ 0 ⃗ H app and to rise with a ramp rate of 0.035 T s −1 , which turned out to be the less consuming computational time ramp rate that allowed us to reproduce the experimental data, i.e.: In the experiment, the shield was in vacuum and cooled through thermal contact with the cooling stage of a cryogenfree cryocooler [13]. To improve this contact all the external surfaces of the sample (Γ 2 , Γ 3 , Γ 4 boundaries in the crosssection scheme in figure 1(b)) were wrapped with an indium layer. Consequently, in the calculations we assumed heat exchange only across the surfaces covered by indium and a constant and homogeneous temperature, T OP , of the indium layer, i.e: whereû n is the unit vector normal to the sample surface and Υ is a heat transfer coefficient determined through iterative adjustments [44,50].

Electromagnetic equations
The electromagnetic behavior of the SC and the surrounding domains is described by Maxwell's equations, which are formulated using the magnetic vector potential ⃗ A: However, as a consequence of the chosen axisymmetric approach, the vector potential has only one component, A ϕ , and equations (5) and (6) become: beingû r andû z the unit vectors along the r-and z-directions, respectively. The relationship between the local electric field and current density (E − J) in the SC domain was modeled by using the hyperbolic tangent function reported in [48], which is a smooth approximation of the stepwise E − J behavior predicted in the critical state model. Accordingly, the relation J = σE presented in the Ampere Law's module in Comsol is customized introducing the electrical conductivity of the superconductor σ SC which takes the form: Here, E 0 is a computational parameter that defines the sharpness of the transition function [48]. Thus, its physical meaning does not correspond to the criterion usually adopted to determine the value of the critical current density from the experimental I-V curves [51] and its value (reported in table 1) was chosen to achieve a good compromise between the computational performance and the accuracy in reproducing the experimental data.
Furthermore, to ensure an appropriate value of the electrical conductivity during the transition between the SC-and normal-state, a term representing the normal state conductivity σ norm was added in parallel to that of the superconductor, as proposed in [43]. We assumed σ norm = 10 8 S m −1 , but it should be noted that we found very minor differences in the magnetic flux density outputs for σ norm values ranging between 10 7 S m −1 and 10 8 S m −1 . The electrical conductivity, σ was then modified as follows: In the equation above, the critical current density, J c (B, T), depends on the local magnetic field and temperature through the following general form [32,43]: where and α, β,α,β, a and b are fit coefficients obtained by the experimental J c0 (T) and B 0 (T) values, in turn achieved by fitting the experimental J c − B curves measured at different temperatures [14]. Their values are listed in table 1.
Although γ is usually assumed to be constant, its value increases with increasing the temperature [13,52]. Therefore, to properly fit the experimental J c (B) curves, the following parabolic behavior of γ as a function of T was assumed: where l and g are still fitting parameters, listed in table 1.

Thermal equation
The changes in the local value of temperature in the SC domain were calculated using the following heat diffusion equation: where Q = ⃗ E · ⃗ J gives the volumetric heating rate and κ(T), C(T) and ρ m are the thermal conductivity, specific heat and mass density, respectively. Due to the lack of data on the thermal properties of our sample, the temperature-dependent behavior of κ(T) and C(T) were included in the model via a piecewise cubic interpolation of the experimental data of the sample HIP#38 reported in [52]. It is worth mentioning that, after a comparison among different curves of κ(T) and C(T) reported in literature-independently of the fabrication techniques [24, 52-60]-, we chose to start from the κ(T) and C(T) curves that best allowed us to reproduce the experimental Normal state conductivity (equation (11)) 10 8 S m −1 α Parameter of J c0 (T) (equation (13)) 4.841 × 10 9 A m −2 α Parameter of J c0 (T) (equation (13)) 1 a Parameter of J c0 (T) (equation (13)) 1.507 β Parameter of B 0 (T) (equation (14)) 1.658 T β Parameter of B 0 (T) (equation (14)) 5 b Parameter of B 0 (T) (equation (14)) 1.694 l Parameter of γ(T) (equation (15)) 0.024 K −1 g Parameter of γ(T) (equation (15)) 0.003 K −2 ρm MgB 2 mass density (measured [14]) 2500 kg m −3 Υ Heat transfer coefficient (equation (3)) 3000 W m −2 K −1 data. Moreover, to also considering that the addition of the hBN powder degrades the thermal properties of the material, the κ(T) values from [52] were divided by 5 to further improve the agreement.

Numerical simulation results and comparison with experimental data
The magnetic flux density measured along the axis of the SC cup-shaped shield as a function of the increasing applied field and reported in [14] evidenced the occurrence of a flux jump at µ 0 H appl = 1.0 T and T OP = 30 K. Figure 2 shows the comparison between these measured magnetic flux density values (symbols) and those calculated with the electromagnetic-thermal coupled approach (dashed lines) assuming the same operating temperature and the parameters listed in table 1. For the sake of completeness, the output achieved with the only electromagnetic model is also reported (solid lines). In this case, only the dependence of J c on the local magnetic field was considered and the parameters J c0 , B 0 , and γ were assumed constant and equal to the outputs of equations (13)-(15) if T = T OP = 30 K. In addition, it is worth mentioning that the values of the magnetic flux density were evaluated at the same fixed positions along the cup's axis, as those of the Hall probes in the experimental setup [14].
As can be seen in figure 2, the coupling of the electromagnetic and thermal equations allows us to well reproduce the effect of the occurrence of the flux jump in the superconductor, which was not predicted when only the electromagnetic equations were considered. In the latter case, a homogeneous field penetration is indeed expected.
To better highlight the ability of the electromagneticthermal model in predicting the shielding property breakdown due to thermomagnetic instabilities, we also calculated the shielding factor (SF) along the axis of the cup, here defined as the ratio between the applied magnetic field µ 0 H app and the local magnetic flux density B z . The SF curves obtained both from measurements (symbols) and numerical calculations (dashed and solid lines) are plotted in figure 3. It is again evident that the numerical calculation performed using the electromagnetic-thermal approach (dashed lines) predicts the abrupt deterioration of the shielding properties after the occurrence of the flux jump at µ 0 H app ≈ 1 T, in remarkable correspondence with the measured values. Conversely, this drastic deterioration of the cup's shielding capability is not reproduced using the electromagnetic-only model (solid lines), which predicts a smooth drop.
The occurrence of a flux jump is a signature of the breaking of the equilibrium between Lorentz and pinning force with the consequent production of a flux avalanche and of Joule heating. Therefore, the correlation of the evolution of the bulk temperature with the magnetic flux density curve and the analysis of the sample temperature maps in the time interval around the flux jump occurrence can give an insight on how to prevent or at least minimize this phenomenon. Figure 4 correlates the mean temperature of the sample with the magnetic flux density calculated in correspondence to position z 1 assuming T OP = 30 K. The superconductor temperature is obtained as the average of the variable T within the SC domain using the Mean Probe function implemented in Comsol. It can be seen that the average temperature of the bulk (left red axis) reaches its maximum, 40.5 K, when the flux jump occurs. At the same time, the magnetic field completely penetrates the superconductor, resulting in a higher value of B z (right blue axis).
The inset in figure 4 shows a detail of the evolution of the temperature peak and the local magnetic field B z as a function of the applied magnetic field. To locally correlate the evolution of temperature, magnetic flux density, and current density during the jump we focused on five different time points, i.e. five different values of the applied magnetic field, near the peak: 0.9765 T (t = t a ), 0.9800 T (t = t b , when the flux jump starts), 0.9817 T (t = t c , when the flux jump ends and the sample has lost its shielding ability), 0.9835 T (t = t d ), and 0.9870 T (t = t e , when the starting temperature is nearly recovered). Figure 5 shows the temperature distribution in the SC domain (color maps) and the ⃗ B(r, z) vector plots at the above-mentioned time  instants. The corresponding current density maps are depicted in figure 6. Figure 5(a) reveals that the temperature starts rising and, consequently, the flux jump originates in a zone along the inner wall of the cup, near its open edge. As expected, in this zone there is a relevant gradient of the magnetic flux density, as highlighted by the ⃗ B vector plot superimposed to the temperature map. Accordingly, in the same region, the current density reaches its maximum value ( figure 6). Considering the evolution of the shielded field, it can be seen that before the flux-jump occurrence the flux lines do not fully penetrate the superconductor ( figure 5(b)) and that the strong increase of   temperature (figure 5(c) induces an instantaneous field penetration in the inner region of the shield. After the flux jump took place, the sample starts cooling down thanks to the heat exchange with the cryocooler cooling stage, gradually reaching the operating temperature T OP = 30 K (figures 5(d) and (e)). However, the shielding properties are not recovered anymore.
As for the temperature and the local magnetic field, an abrupt change in current density value can be observed when the flux jump occurs (t = t c , figure 6(c)). Actually, as a consequence of the flux jump, the value of the current density decreases by more than three orders of magnitude. It starts to increase again close to the external wall of the shield when the temperature of the bulk stabilizes back to 30 K (t = t d and t = t e , figures 6(d) and (e)).
To complete the validation of our modeling approach, it was also applied to calculate the magnetic flux density along the shield's axis at 25 and 20 K, then comparing computational and experimental data. As obtained for T OP = 30 K, the coupling between the electromagnetic and the thermal model makes it possible to well reproduce the occurrence of the flux jumps detected during the experimental characterization also at these working temperatures (figures 7 and 8), including the double flux jump occurrence at 20 K and µ 0 H app = 1.82 and 2.52 T. For the sake of comparison, the curves achieved using the only electromagnetic model are plotted as well.

Influence of the thermal conductivity and heat exchange on the shielding properties
The as-validated modeling approach allows studying possible solutions and improvements to avoid or lessen the thermomagnetic instability occurrence or to shift them to higher applied fields. It was shown that a greater thermal exchange with the cooling system can increase the stability against flux jumps [38]. Moreover, using precursors that ensure better thermal properties of the MgB 2 bulk could prevent the occurrence of thermomagnetic instabilities as well [16,36]. For this reason, focusing on T OP = 30 K, we numerically evaluated how enhancements in MgB 2 thermal conductivity, and/or in the heat exchange with the cooling stage can affect the cup's SF. The latter point was addressed by adding a heat exchange (cooling) also across the inner lateral surface of the cup, i.e. imposing the boundary condition (equation (3)) also on the Γ 5 boundary ( figure 1(b)). Figure 9 compares the B z (figure 9(a)) and SF (figure 9(b)) curves calculated in section 3 in correspondence to positions z 1 and z 3 -i.e. considering κ(T) values from [52]   • κ(T) values from [52] divided by 1 and the heat exchange as detailed in equation (3) • κ(T) values from [52] divided by 5 and the heat exchange (see equation (3)) on boundaries Γ 2 ∪ Γ 3 ∪ Γ 4 ∪ Γ 5 • κ(T) values from [52] divided by 1 and the heat exchange (see equation (3)) on boundaries Γ 2 ∪ Γ 3 ∪ Γ 4 ∪ Γ 5 .
As can be seen, assuming higher value of κ(T) shifts the flux jump occurrence to slightly higher values of the applied field. Hence, the SF worsens at 1.02 T assuming κ(T)/1 instead of 0.98 T with κ(T)/5 at z 1 position (black curves). On the other hand, the addition of a heat exchange across the inner wall of the shield totally prevents the  [52] divided by 5 and by 1, respectively, and heat exchange across boundaries Γ 2 ∪ Γ 3 ∪ Γ 4 ∪ Γ 5 . Black and green curves refer to z 1 and z 3 positions, respectively. magnetothermal instabilities for both the κ(T) values taken into account.

Conclusion
In this work, we investigated the reliability of a computational approach in predicting the effects of the flux-jump occurrence on the magnetic shielding performance of a bulk superconductor shield. To this aim, we coupled an electromagnetic model exploiting an ⃗ A-based formulation with the heat diffusion equation. The model output was then compared with the corresponding experimental data attained on a cup-shaped MgB 2 shield [14].
At first, we focused on an operating temperature of 30 K and showed that the coupled numerical model quantitatively reproduces the experimental SF curves and their fall due to flux jump occurrence. Next, we correlated this abrupt worsening of the shielding performance inside the SC cup with the evolution of the local magnetic field, temperature, and current density in the superconductor. In doing so, we determined the locations in the sample where flux avalanches are expected to nucleate. This thermomagnetic model was then applied to reproduce the shielding curve at 20 and 25 K, reaching a remarkable agreement between the experimental and calculated results.
Finally, we exploited this numerical procedure to investigate the influence of the thermal exchange with the coolant and the thermal conductivity of the material on the shielding properties. It emerges that an improvement of the thermal conductivity has a more limited impact on the SF, whereas a better thermal contact could even prevent the flux jump occurrence.
It is worth underlining that the use of a thermomagnetic model is a prerequisite to realistically predict the shielding performance of a high-J c MgB 2 bulk up to quite high temperatures, as also evidenced by the comparison of the outputs of the thermomagnetic model and of an electromagnetic-only model. With this in mind, additional refinements of the model as well as further investigations on the material properties are still underway, also including hybrid SC-ferromagnetic layouts.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.