Fast and efficient critical state modelling of field-cooled bulk high-temperature superconductors using a backward computation method

A backward computation method has been developed to accelerate modelling of the critical state magnetization current in a staggered-array bulk high-temperature superconducting (HTS) undulator. The key concept is as follows: i) a large magnetization current is first generated on the surface of the HTS bulks after rapid field-cooling (FC) magnetization; ii) the magnetization current then relaxes inwards step-by-step obeying the critical state model; iii) after tens of backward iterations the magnetization current reaches a steady state. The simulation results show excellent agreement with the H-formulation method for both the electromagnetic and electromagnetic-mechanical coupled analyses, but with significantly faster computation speed. Solving the FEA model with 1.8 million degrees of freedom (DOFs), the backward computation method takes less than 1.4 hours, an order of magnitude or higher faster than other state-of-the-art numerical methods. Finally, the models are used to investigate the influence of the mechanical stress on the distribution of the critical state magnetization current and the undulator field along the central axis.


Introduction
Research and development work on short-period and high-field staggered-array high-temperature superconducting (HTS) bulk undulators [1,2] is ongoing in a European project for the construction of compact free electron lasers (FELs) [3,4]. This new technology utilizes a 10 T level superconducting solenoid magnet to realize field-cooling (FC) magnetization of a series of staggered-array ReBCO bulks at a temperature around 10 K. With this concept, a sinusoidal undulator field of amplitude 2 T with a period length of 10 mm along the central beam axis can be obtained [2,5]. One key challenge for this technology is that the mechanical properties of ReBCO bulk superconductors are ceramic-like: friendly towards compressive stress and unfriendly towards tensile stress. Thus, some form of external mechanical reinforcement is usually used to compress the ReBCO bulk to trap high magnetic fields [6][7][8][9]. Trillaud et al. (2018) showed the critical current density Jc of ReBCO bulk superconductor will degrade when the Lorentz force-induced mechanical stress is of the order of the fracture strength [10]. This indicates that the critical current density Jc should be a function of both the magnetic flux density B and the mechanical strain ε when a ReBCO bulk superconductor traps a high magnetic field and experiences the associated large Lorentz force. Regarding the short-period and high-field staggered-array HTS bulk undulator, estimation of the magnetization current that follows a Jc(B,ε)-determined critical state model [11,12], without time-dependent flux creep effects, is of great interest for the purpose of optimizing the first and the second integrals of the undulator field along the central axis.
Recently an iterative algorithm method was proposed to compute a Jc(B,θ)-determined critical state model for ReBCO tape stacks [22]. It avoids using unnecessary iterative steps to obtain a resistivity matrix [23,24] but still requires hundreds of iterative steps to obtain adequate results. This paper introduces a new backward computation method to accelerate modelling the Jc(B,ε)-determined critical state magnetization current in the periodical HTS bulk undulator. It takes only tens of backward iterations to reach a steady-state solution, which can be an order of magnitude or higher faster than other state-of-the-art numerical methods. Figure 1 shows the one-period 2D FEA model of the periodical HTS bulk undulator created in ANSYS 18.1 Academic. For the electromagnetic analysis, the magnetic flux density Bx is applied to the outer air subdomain boundaries @y = ± 315 mm to provide the background magnetic field; a flux normal boundary (default in ANSYS) is applied to the boundaries on the sides @x = ± 315 mm to model the periodicity/symmetry. For the mechanical analysis, a displacement constraint is applied the x-direction @x = 0 and x = ± 5 mm and the y-direction @y = ± 8.5 mm to avoid movement due to the action of the Lorentz force. The pre-stress, if any, is applied to the top and bottom sides of the HTS bulks as mechanical reinforcement to compensate   and then reduced from 10 T to zero over a short time (50 seconds; step 2). The eddy current solver is turned off during the first step and turned on during the second step to calculate the induced magnetization current. The resulting mechanical stress is analyzed by importing the Lorentz force and applying a pre-stress. Afterwards, a backward loop computation of the relaxation of the magnetization current is carried out as follows a) Obtain the magnetization current JT and the equivalent mechanical strain εeq for each HTS element, and update Jc; b) For each HTS element, force the magnetization current JT to Jc•JT/|JT| if |JT| > Jc or the element has been penetrated (Each HTS element has a "label" with the default value of zero; once the HTS element is penetrated its "label" becomes 1); c) Carry out the transient electromagnetic analysis with a small time increment (Δt = 0.5 s); d) Carry out the 2D plane strain mechanical analysis.

FEA model and backward computation
During the backward iterations the resistivity of the superconductor is set to a fixed low value (1x10 -15 Ωm) and the A-V formulation is used for fast and efficient electromagnetic analysis.
The entire process follows these Maxwell's equations and the modified critical state model for which Jc @10 K is expressed as where kc,m is the mechanical degradation factor describing the Jc degradation due to the mechanical stress. The assumed values of Jc1, Jc1, BL, Bmax and y are 1.0x10 10 A/m 2 , 8.8x10 9 A/m 2 , 0.8 T, 4.2 T and 0.8, respectively. These values refer to the Jc data [25] of the ReBCO bulk @ 40 K and are scaled to 10 K from the first experimental result of our ten-period GdBCO bulk undulator tested at the University of Cambridge [2]. The mechanical degradation factor [10] is a function of the equivalent mechanical strain ε eq where γ = 0.1, β = 0.025, α = 10% and εc = 6.0x10 -4 (σc = 90 MPa, E = 150 GPa).

Computational results from the backward computation method
The simulation results for three different cases are compared and plotted in figure 3 (multimedia view: load step 1-5, ramp Bx from zero to 10 T; load step 6, damp Bx from 10 T to zero; load step 7, start backward iterations). In figure   The simulation results confirm two facts: a) applying a pre-stress on the bulk HTS can enhance its mechanical performance for the purpose of trapping high magnetic field; b) the applied pre-stress can affect the distribution of the magnetization current in the bulk HTS, thus reducing the undulator field along the central axis.

Validation by the electromagnetic-mechanical coupled H-formulation
In this section, the electromagnetic properties of the bulks are simulated using the H-formulation, implemented in COMSOL Multiphysics (version 5.4) using the 'Magnetic Field Formulation' interface in COMSOL's AC/DC module. Both COMSOL and the H-formulation are currently used by dozens of research groups worldwide to model bulk superconductors [26,27] and other superconductivity-related problems [28,29].  ∇ × H = J (7) The permeability µ = µ0, and equations (7) and (8) are combined with the E-J power law (9), used to simulate the nonlinear   Ez] are the current density and electric field, respectively, which are assumed to be parallel to each other such that E = ρJ. E0 = 1 μV/cm is the characteristic electric field and n defines the steepness of the transition between the superconducting state and the normal state; we assume here that n = 100 to reasonably approximate the critical state model [27,33].
Isothermal conditions are assumed while ramping down the field, so no thermal model is included.
The electromagnetic model is coupled with COMSOL's 'Solid Mechanics' interface as described in [34]. The Lorentz force, FL = J x B, is implemented as a force per unit volume using the 'Body Load' node, where Fx = -Jz·By and Fy = Jz·Bx. The displacement constraints are added using the 'Prescribed Displacement' node. The pre-stress is applied using the 'Boundary Load' node such that Fpre = -pn, where p is the applied pressure.  Figure 6 shows the mechanical stress in the periodical HTS bulk undulator at "t= 100 s". Overall, the simulation results are highly consistent with those shown in figure 4.

Discussion
The backward computation method has been proven successful to model the critical state magnetization current in the HTS bulk (c) The backward computation method solves the eddy current in the HTS subdomains by defining a fixed low resistivity.
Solving an equation representing the nonlinear resistivity, such as the E-J power law, is not required.
In order to demonstrate the high efficiency of the backward computation method, we conducted two identical simulations  Table 1. For the electromagnetic-only analysis, the backward computation is ~4 times faster than the H-formulation; for running the electromagnetic-mechanical coupled analysis, the backward computation is ~14 times faster than the H-formulation. The rapid coupling analysis suggests the backward computation can solve any modified critical state models without reducing computation speed. It should be pointed out that the COMSOL H-formulation was run on a PC with an To better understand the efficiency of the backward computation, we conducted two additional simulations by reducing the element size (increasing the total number of DOFs) in the HTS undulator model. As shown in Table 2, the total number of DOFs increases to 452,034 and then 1,799,210 when the element size is "0.125 mm x 0.125 mm" and "0.0625 mm x 0.0625 mm", respectively. It took approximately 15 minutes and 82 minutes, respectively, to run the electromagnetic-mechanical coupled simulation. We have compared these results with other state-of-the-art techniques for the electromagnetic analysis of HTS materials, like the H-formulation [35], the T-formulation [36], the variational method [15], the T-A formulation [37], and the recently proposed iterative algorithm method [22,38]. As shown in figure 7, the backward computation shows a surprising order of magnitude computation speed than all the other methods. It should be noted, however, that the listed H-formulation, T-formulation and T-A formulation were implemented for other applications (e.g., AC loss or screening-current-induced fields) and that benchmarking this particular problem would provide a true comparison. Nevertheless, solving such a large-scale HTS electromagnetic problem with 1.8 million DOFs within 1.4 hours is remarkably fast and was achieved using a normal PC.

Conclusion
We have demonstrated that the backward computation method can model the critical state magnetization current in a staggered array HTS bulk undulator quickly and efficiently by running benchmark simulations using a mechanical-coupled H-formulation in COMSOL. The algorithm of the backward iterations is realized by utilizing the function of multi-frame restart analysis and the A-V formulation in ANSYS 18.1 Academic. The highly efficient computation, even with millions of DOFs, is because a nonlinear resistivity equation is not required and no eddy current is solved in the non-superconductor regions. These advantages, along with the backward concept itself, make this new method superior to many other numerical methods used to model the critical state magnetization current. Finally, we show that applying a pre-stress to the HTS bulks could enhance their mechanical performance when trapping high magnetic fields, but could result in a reduction in Jc in the outer layer of the HTS bulks, thus reducing the induced undulator field. This important information will help guide future optimization of the integral undulator field along the meters-long central axis.