Numerical modelling of mechanical stresses in bulk superconductor magnets with and without mechanical reinforcement

The magnetic field trapping capability of a bulk superconductor is essentially determined by the critical current density, Jc(B, T), of the material. With state-of-the-art bulk (RE)BCO (where RE = rare earth or Y) materials it is clear that trapped fields of over 20 T are potentially achievable. However, the large Lorentz forces, FL = J × B, that develop during magnetisation of the sample lead to large mechanical stresses that can result in mechanical failure. The radial forces are tensile and the resulting stresses are not resisted well because of the brittle ceramic nature of (RE)BCO materials. Where fields of more than 17 T have been achieved, the samples were reinforced mechanically using resin impregnation and carbon-fibre wrapping or shrink-fit stainless steel. In this paper, two-dimensional (2D) axisymmetric and three-dimensional (3D) finite-element models based on the H-formulation, implemented in the commercial finite element software package COMSOL Multiphysics, are used to provide a comprehensive picture of the mechanical stresses in bulk superconductor magnets with and without mechanical reinforcement during field-cooled magnetisation. The chosen modelling framework couples together electromagnetic, thermal and structural mechanics models, and is extremely flexible in allowing the inclusion of various magnetisation processes and conditions, as well as detailed and realistic properties of the materials involved. The 2D model—a faster route to parametric optimisation—is firstly used to investigate the influence of the ramp rate of the applied field and any heat generated in the bulk. Finally, the 3D model is used to investigate the influence of inhomogeneous Jc(B, T) properties around the ab-plane of the bulk superconductor on the developed mechanical stress.


Introduction
Bulk superconductors, acting as trapped field magnets, can trap magnetic fields over ten times larger than the maximum fields produced by conventional permanent magnets, which are limited practically to rather less than 2 T [1]. It has been shown that (RE)BCO (where RE = rare earth or Y) bulk superconductors can trap fields greater than 17 T: 17.24 T at 29 K [2] and 17.6 T at 26 K [3] have been demonstrated. The magnetic field trapping capability of a bulk superconductor is essentially determined by the critical current density, J c (B, T), of the material. With state-of-the-art bulk (RE)BCO materials it is clear that trapped fields of over 20 T are potentially achievable. However, the large Lorentz forces, F L = J × B, that develop during magnetisation of the sample lead to large mechanical stresses that can result in mechanical failure, with unreinforced samples typically failing for magnetic fields greater than 7-9 T [4,5]. The radial forces are tensile in nature [6] and the resulting stresses are not resisted well because of the brittle ceramic nature of (RE)BCO materials [7]. To achieve the record trapped fields >17 T, the samples were reinforced mechanically using resin impregnation and carbon-fibre wrapping [2] and shrink-fit stainless steel [3].
Extending numerical models developed to date to investigate the magnetisation of bulk superconductors, which have been primarily focused on the electromagnetic and thermal analyses [8], there has been a great deal of interest recently in simulating and analysing the mechanical properties of bulk superconductors using numerical tools [9][10][11][12][13][14]. These studies have been predominately based on the A-formulation used by Fujishiro et al [9][10][11][12][13], assuming a constant temperature, but the H-formulation was recently used to analyse the stresses developed during pulsed field (a) 2D axisymmetric model of a disc-shaped bulk superconductor with a mechanical reinforcement ring in the r-z-Φ cylindrical coordinate system, making additional use of symmetry along the r-axis to model half the cross-section. (b) 3D model in the x-y-z coordinate system, where geometric symmetry is used to model only 1/8th of the bulk: 1/4 of the bulk around the ab-plane and 1/2 of the bulk along the c-axis (z-axis). In this work, a diameter of 30 mm and thickness of 15 mm are assumed as the dimensions of the bulk, and the reinforcement ring is assumed to be 10 mm wide and the same thickness as the bulk.  [24,25] for fields up to 10 T over a temperature range of 40-85 K, extended here to 20 T using the equation presented by Jirsa et al in [31]. The data is input into the model using a two-variable, direct interpolation, as described in [32,33]. Figure 3. Time-dependence of the temperature, T(t), as well as the applied magnetic field, B z (t), during the FCM process as simulated in these models. magnetisation (PFM) by Wu et al [14], which was coupled with a thermal model to include the influence of heat generated during the PFM process.
In this paper, two-dimensional (2D) axisymmetric and three-dimensional (3D) finite-element models based on the H-formulation, implemented in the commercial finite element software package COMSOL Multiphysics [15], are used to provide a comprehensive picture of the mechanical stresses in bulk superconductor magnets with and without mechanical reinforcement during field-cooled magnetisation (FCM). The modelling framework couples electromagnetic, thermal and structural mechanics models, and is extremely flexible to include various magnetisation processes and conditions, as well as detailed and realistic properties of the materials involved. In section 2, the modelling framework is introduced and the results of the 2D and 3D models are compared to demonstrate the consistency of the two models, as well as to show the effect that the reinforcement has on the hoop stress in the bulk during FCM. In section 3, the 2D model-a faster route to parametric optimisation-is firstly used to investigate the influence of varying the ramp rate of the applied field and any heat generated in the bulk. Finally, the 3D model is used to investigate the influence of inhomogeneous J c (B, T) properties around the ab-plane of the bulk superconductor on the developed mechanical stress.

Modelling framework
Numerical techniques based on the finite element method (FEM) have been applied to many superconducting material problems using a variety of formulations [8]. Each of these formulations are equivalent in principle, i.e. the choice of formulation should not result in a different physical meaning of the solution, but the solutions of the corresponding partial differential equations can be very different [16]. For more detailed information, including techniques not based on FEM, the reader may refer to recent review papers on numerical methods for calculating AC losses in high-temperature superconducting (HTS) materials [17], the modelling of bulk superconductor magnetisation [8] and HTS applications [18]. Trapped field at the centre of the bulk and at the centre of the top surface at z=+0.5 mm above the bulk during the ramp down of the applied field in the FCM process (t=1250 s), for the 2D axisymmetric and 3D models, under the following magnetising conditions: T op =50 K, B app =20 T, dB down /dt=25 mT s −1 . Table 1. List of assumed material properties for the coupled electromagnetic-thermal-mechanical modelling of mechanical stresses in bulk superconductor magnets with and without mechanical reinforcement.

Parameter Description
Value In-field, temperature-dependent critical current density Interpolation (see figure 2)  Poission's ratio (stainless steel) 0.28 [10][11][12] A bulk superconductor is usually fabricated in the form of a cylindrical disc and this geometry lends itself well to simplification using a 2D axisymmetric model (r, z, Φ in a cylindrical coordinate system), assuming that its properties are homogeneous in the Φ-direction, i.e. around the bulk's ab-plane. However, a 3D geometry (x, y, z coordinates) is required when including any inhomogeneous material properties around the ab-plane, e.g. differences in J c between growth sector boundaries (GSBs) and growth sector regions (GSRs) [19], and for shapes without cylindrical symmetry,  such as rectangular-shaped bulks [20]. It can also be useful to exploit symmetry, where possible, in a 3D model by applying appropriate boundary conditions to model, for example, one half (2D or 3D) or one quarter of the bulk. This can reduce the number of mesh elements required and improve both computational memory requirements and speed, while retaining the necessary 3D effect. Figure 1(a) shows a 2D axisymmetric model of a disc-shaped bulk superconductor in the r-z-Φ cylindrical coordinate system, making additional use of symmetry along the r-axis to model half the cross-section. The cross-section (Φ-direction; rz-plane) corresponds to the ab-plane of the bulk and the z-axis corresponds to the c-axis. Figure 1(b) shows a 3D model in the x-y-z coordinate system, where geometric symmetry is used to model only 1/8th of the bulk: 1/4 of the bulk around the ab-plane and 1/2 of the bulk along the c-axis (z-axis). In this work, a diameter, D, of 30 mm and thickness, H, of 15 mm are assumed as the dimensions of the bulk, which is a typical size of those fabricated by our research group, and the reinforcement ring is assumed to be 10 mm wide and the same thickness as the bulk.
The bulk's electromagnetic properties are simulated using the H-formulation, implemented in the commercial software package COMSOL Multiphysics 5.3a [15]. This framework has been used previously by the authors to simulate bulk (RE)BCO materials under various magnetisation conditions [19,[21][22][23][24][25][26], as well as FCM of MgB 2 [27] and iron-pnictide [28] bulks. In the H-formulation, the governing equations are derived from Maxwell's equationsnamely, Ampere's (1) and Faraday's (2) laws: the permeability of free space, and for the superconducting and air sub-domains, the relative permeability can be assumed as simply μ r =1.
The E-J power law [29,30] is used to simulate the highly non-linear resistivity of the superconductor, where E is proportional to J n and n=20 is assumed as a typical value for HTS materials and a good approximation of Bean's critical state model [26].
The results of the numerical simulation depend strongly on the assumed J c (B, T) characteristics of the material, and in this paper, the measured J c (B, T) characteristics of a representative bulk HTS (15 wt% Ag-containing GdBa 2 Cu 3 O 7−δ ), presented in [24,25] for fields up to 10 T over a temperature range of 40-85 K, is extended to 20 T using the equation presented by Jirsa et al in [31]: Figure 2 shows the assumed J c (B, T) characteristics, which are input into the model using a two-variable, direct interpolation, as described in [32,33]. Since FCM is being simulated, there are three distinct steps in the numerical model [27]: (1) applying a ramped external magnetic field parallel to the c-axis of the bulk up to a maximum magnitude, B app , by setting appropriate magnetic field boundary conditions [34] such that H z (t)=H app (t/t ramp ) for tt ramp , where H app =B app /μ 0 and t ramp is the duration of the ramp, while the temperature is held at T>T c =92 K (T=300 K in this paper); (2) cooling the bulk to an appropriate operating temperature, T op =T c , while the external field is held at B app ; and (3) once the temperature has stabilised, ramping the external field down from B app to zero. Hence, J c must also be defined for TT c =92 K. Here, J c is assumed to be 1×10 6 A m −2 for T92 K. This is also true for when J c <1×10 6 A m −2 for T between 60 K and 85 K in figure 2.
Since the temperature of the bulk changes significantly during FCM, the electromagnetic model is coupled with a thermal model, based on the following thermal transient equation: The models in this section are coupled through the J c (B, T) characteristics described earlier and isothermal conditions are assumed while ramping down the field, i.e. Q=0, because the magnetisation process is slow. It is possible to further couple the electromagnetic and thermal models by including any heat generated during the magnetisation process by assuming This was recently performed in [14] to analyse the stress during PFM and will be used in section 3.1. Figure 3 shows the time-dependence of the temperature, T(t), as well as the applied magnetic field, B z (t), during the FCM process as simulated in these models. It should be noted that in section 3.1, the effect of heat generation during ramp down of the field is investigated by including Q=E·J and varying the ramp rate.
Finally, the mechanical stresses during FCM are derived from the following, principal governing equation, an expression of Newton's second law in direct tensor form: where σ s is the Cauchy stress tensor and F L is the Lorentz force (described below). Thus, the structural transient behaviour is assumed to be quasi-static and the second order timederivatives of the displacement variables, u-the inertial terms, ρd 2 u/dt 2 -are zero. The strain-displacement relationship is given by where ε is the infinitesimal strain tensor. The bulk is assumed to be an isotropic, linear elastic material, obeying Hooke's law,   such that σ s =C(E, ν): ε; C is the fourth-order stiffness tensor.
Here the mechanical properties of the bulk are defined by its Young's modulus, E bulk =1×10 11 Pa, and Poisson's ratio, ν bulk =0.33 [9][10][11][12][13], and a density, ρ bulk =5900 kg m −3 , is assumed [35]. The stainless-steel reinforcement ring, where present, is also assumed to be an isotropic, linear elastic material with E sus =1.93× 10 11 Pa, ν sus =0.28 and ρ sus =8000 kg m −3 [10][11][12] and, by default, a perfect mechanical contact between the ring and the bulk is assumed. The Lorentz force, F L =J×B, that develops during the FCM process, derived from the electromagnetic model, is used as the input to calculate the mechanical stresses. Table 1 lists all of the assumed material properties for the coupled electromagnetic-thermal-mechanical modelling, including applicable references to other work in the literature.

Preliminary results (2D axisymmetric and 3D models)
In this section, the results of the 2D axisymmetric and 3D models are compared, both with and without mechanical reinforcement with a stainless-steel ring, to demonstrate the consistency of the two models, as well as to show the effect that the reinforcement has on the hoop stress in the bulk during FCM.
2.2.1. Without ring reinforcement. Figure 4 shows the trapped field at the centre of the bulk and at the centre of the top surface at z=+0.5 mm above the bulk during the ramp down of the applied field in the FCM process, for the 2D axisymmetric and 3D models, under the following magnetising conditions: T op =50 K, B app =20 T, dB down / dt= 25 mT s −1 . At the end of the FCM process (t=1250 s), the trapped field at the centre of the top surface (z = +0.5 mm), B t , is approximately 7.6 T, and the trapped field at the centre, B c , is approximately 11.4 T. Figure 5 shows a comparison of the magnetic flux density, |B|, and current density, |J|, distributions throughout the cross-section of the bulk at the end of the FCM process (t=1250 s), for the 2D axisymmetric and 3D models, under the same magnetising conditions. Figure 6 shows a similar comparison of the Lorentz force in the radial direction, F L,r , and hoop stress, σ j , distributions. The results show clear consistency between the two models and the cross-sectional plots give a clear indication of the location of maximum Lorentz force, as well as maximum hoop stress, which is highest at the centre of the bulk. Figure 7 shows the time-dependence of the hoop stress, σ j , across the centre of the bulk (z=0) at discrete points every 4 T from during the ramp down of the applied field in the FCM process. The results clearly show the dynamic nature of the generated hoop stress, which is directly related to the Lorentz force generated by the distribution of the simultaneously reducing magnetic field, B, and increasing induced current, J, which takes a maximum when the field has ramped down to around 4 T under the given magnetising conditions and material properties. The influence of the reinforcement on the dynamic nature of the hoop stress, as well as its maximum value, is investigated in the next section.

With ring reinforcement.
In this section, the same magnetising conditions are assumed (T op =50 K, B app =20 T, dB down /dt=25 mT s −1 ), but a stainless-steel ring of width 10 mm (inner radius, r inner , of 15 mm, to match the radius of the bulk; outer radius, r outer , of 25 mm) and the same thickness as the bulk (15 mm) is added around the periphery of the bulk to provide mechanical reinforcement. The resultant plots of the trapped field (see figure 4), as well as |B| and |J| (see figure 5) and F L,r (see figure 6) are the same as for the unreinforced case in section 2.2.1; however, the resultant hoop stress changes because of the contribution from the thermal compressive stress applied to the bulk by the reinforcement ring when cooled during FCM (before ramping down the applied field). This occurs due to the difference in the thermal coefficient of expansion between the two materials. A comparison of the hoop stress distributions throughout the cross-sections of the bulk and reinforcement ring at the end of the FCM process (t=1250 s), for the 2D and 3D models, is shown in figure 8. Again, the results show clear consistency, and it is also clear that the bulk is in compression (negative σ j ), whereas the reinforcement ring is in tension (positive σ j ), at the end of the FCM process. Figure 9 shows the time-dependence of the hoop stress, σ j , across the centre of the bulk and reinforcement ring (z=0) at the beginning (20 T) and end (0 T) of the ramp down of the applied field in the FCM process when the bulk is mechanically reinforced. The results at 4 T, corresponding to the maximum σ j in the unreinforced case (see figure 7), are Figure 11. Trapped field at the centre of the top surface at z=+0.5 mm above the bulk (without reinforcement) during the ramp down of the applied field in the FCM process for ramp-down rates of 2.5, 25 and 250 mT s −1 , for T op =50 K and B app = 20 T. The heat source, Q=E · J, is introduced to include the influence of any heat generated and the bulk is cooled from its bottom surface. also shown for comparison. It should be noted that σ j in this case includes both contributions from σ j (FCM), the hoop stress related to the Lorentz force, and σ j (COOL), the thermal compressive stress applied to the bulk by the ring when cooling to T op . The results suggest that for these particular magnetising conditions, there is adequate reinforcement, as the centre of bulk (where the largest hoop stress occurs in the unreinforced case) is in compression, as is the rest of the bulk. Figure 10 shows the time-dependence of the maximum hoop stress, σ j,max , anywhere in the cross-section/volume of the bulk, with and without mechanical reinforcement, as the field is ramped down. For the unreinforced bulk, σ j,max corresponds to the maximum value at the centre of the bulk as shown in figure 7: σ j,max = 63.3 MPa at approximately 4 T. However, for the reinforced bulk, σ j,max (4.9 MPa at approximately 6 T) is not located at the centre of the bulk, but at the edge of the bulk near the surface, which can be seen in figure 8. This was also observed in numerical models for reinforced ring-shaped bulks in [12], and careful optimisation of the ring geometry compared to the bulk is necessary to avoid local stresses and improve the uniformity of the effect of the mechanical reinforcement, which could be carried out easily using this numerical framework.

Influence of system parameters
The 2D axisymmetric model provides a faster route to parametric optimisation in comparison to the 3D model, if the assumption that the material properties are homogeneous around the bulk's ab-plane holds. In the following sections, the effect of varying the ramp rate of the applied field and any heat generated in the bulk is investigated using the 2D axisymmetric model, and then the 3D model is used to investigate the influence of an inhomogeneous J c (B) distribution around the bulk's ab-plane.
3.1. Influence of ramp rate and heat generated 3.1.1. Without ring reinforcement. In the simulations so far, isothermal conditions have been assumed, i.e. Q=0 in equation (4), neglecting the influence of any heat generated during the ramp down of the applied field. In general, the ramp-down rate in the FCM process is usually slow and the heat generated is nowhere near as extreme as that experienced during PFM [14]. However, heat generated during any magnetisation process limits the trapped field capability of the material and it is important to understand how much heat may be generated, as well as its influence on the mechanical stresses in the material. In this section, the heat source, Q=E·J, is introduced to the thermal model (equation (4)), and the influence of three ramp down rates (2.5, 25 and 250 mT s −1 ) on the FCM process is investigated. It is assumed, as done in section 2.2, that T op =50 K and B app =20 T.
In the previous analyses, the temperature of the bulk was set by a temperature constraint across the centre of the bulk; however, when including the heat generated during FCM, one must apply a more realistic assumption regarding cooling and here it is assumed that the bulk is cooled from its bottom surface (a typical scenario when magnetised in solenoidal high field magnets), necessitating the removal of the symmetry along the centre of the bulk and modelling of the whole cross-section of the bulk (in the 2D axisymmetric model). It is also possible to cool the sample from the periphery of the bulk for a split coil magnetisation fixture [24]. For simplicity, perfect thermal contact between the bulk  and the cryocooler is also assumed; however, as shown in [24,25,35], it is possible to adapt the model to include more realistic assumptions regarding this thermal contact, as well as the finite cooling power of the cryocooler. Figure 11 shows the trapped field at the centre of the top surface at z=+0.5 mm above the bulk during the ramp down of the applied field in the FCM process for ramp-down rates of 2.5, 25 and 250 mT s −1 , for T op =50 K and B app =20 T. As the ramp-down rate is increased, there is a reduction in the trapped field due to a reduction in J c , which is most noticeable for the fastest ramp rate. Figure 12 shows the temperature of the bulk at the end of the FCM process at t=for 450, 1250 and 8450 s for the three ramp rates, 2.5, 25 and 250 mT s −1 , respectively. There is almost no temperature rise for the slowest ramp rate (2.5 mT s −1 ), and for the intermediate ramp rate (25 mT s −1 ), there is a temperature rise of around 3 K towards the upper surface of the bulk. In the case of the fastest ramp rate (250 mT s −1 ), there is a significant temperature rise of up to 18 K in the same region.
The impact of this temperature rise on the magnetic flux and current density distributions is shown in figures 13 and 14, respectively, which accounts for the reduced trapped fields in figure 11, except for the 2.5 mT s −1 case. In this case, although the temperature rise is negligible, because the ramp down occurs over a longer time scale (t=450→8450 s), there is more flux creep, resulting in a slightly lower trapped field at the end of the FCM process. The 250 mT s −1 rampdown rate results in a significantly distorted current density/ trapped field distribution, and although there is a small temperature rise for the 25 mT s −1 case, which deviates slightly from the previous assumption of isothermal conditions, a significant impact is not observed.
The impact on the developed hoop stress, σ j , within in the bulk is shown in figure 15 at various positions within the bulk, z=0 and ±6 mm, for the three ramp-down rates, including 25 mT s −1 under the assumption of isothermal conditions. The faster ramp rates, resulting in a higher temperature rise and subsequently reduced J c , see lower stresses developed overall, which agrees well with the analysis in the previous section. In addition, there is a larger distortion of the stresses between the top and bottom regions of the bulk. In the next section, the impact of the reinforcement ring on these findings is explored.

3.1.2.
With ring reinforcement. The stainless-steel reinforcement ring is now added around the periphery of the bulk and the same magnetising conditions as the previous section are assumed: T op =50 K, B app =20 T and three ramp-down rates (2.5, 25 and 250 mT s −1 ), including the heat source, Q=E·J. It is also assumed again that the thermal contact between the bulk, ring and cryocooler is ideal. Figure 16 shows the trapped field at the centre of the top surface at z=+0.5 mm above the bulk during the ramp down of the applied field in the FCM process for ramp-down rates of 2.5, 25 and 250 mT s −1 , for T op =50 K and B app =20 T. As seen in the previous section, as the rampdown rate is increased, there is a reduction in the trapped field due to a reduction in J c , which is most noticeable for the fastest ramp rate. However, there is notably less reduction in trapped field with the reinforcement ring, compared to without one. Figure 17 shows the temperature of the bulk with reinforcement at the end of the FCM process at t=for 450, 1250 and 8450 s for the three ramp rates, 2.5, 25 and 250 mT s −1 , respectively. There is almost no temperature rise for the slowest ramp rate (2.5 mT s −1 ), and for the intermediate ramp rate (25 mT s −1 ), there is a temperature rise of around 1.2 K towards the upper surface of the bulk (3 K without reinforcement). In the case of the fastest ramp rate (250 mT s −1 ), the temperature rise is less than 10 K (18 K without reinforcement) in the same region. The presence of the reinforcement ring provides an additional thermal pathway that allows cooling through the ab-plane of the bulk, for which the thermal conductivity is several times larger than along the c-axis [9].
The impact of the reduced temperature rise on the magnetic flux and current density distributions is shown in figures 18 and 19, respectively, and it is clear that there is less distortion of the current density/trapped field distribution. There is then significantly less distortion of developed hoop stress, σ j , within in the bulk, which is shown in figure 20 at various positions within the bulk, z=0 and ±6 mm, for the three ramp-down rates, including 25 mT s −1 under the assumption of isothermal conditions. It should be noted again that σ j in this case includes both contributions from σ j (FCM), the hoop stress related to the Lorentz force, and σ j (COOL), the thermal compressive stress applied to the bulk when cooling to T op . There is only a perceptible difference in σ j for the fastest ramp rate (250 mT s −1 ), meaning that the reinforcement ring not only provides compression of the bulk during the FCM process to avoid mechanical fracture, but provides an additional thermal pathway to remove the heat generated during the magnetisation process, which is particularly useful in the case of PFM.

Inhomogeneous J c distribution around the ab-plane
Although the 2D axisymmetric model provides a faster route to parametric optimisation, this assumes that the material properties are homogeneous around the bulk's ab-plane. In this section, the 3D model is used to investigate an example of the influence of inhomogeneous material properties around the ab-plane, such as those that occur during the growth process of c-axis seeded, single grain (RE)BCO bulk superconductors [19]. The 3D model is also useful for shapes without cylindrical symmetry, such as rectangular-shaped bulks [20].
Following the same method presented in [19,36], the J c (B, T) characteristics are varied using a cosine function around the ab-plane, supposing some difference of ±10% in J c between the GSBs and GSRs, such that Given the relevant material properties are known, by measuring a number of sub-specimens from a bulk superconductor, any position-dependent J c could be introduced in the model.
The same magnetising conditions as the previous sections are assumed: T op =50 K, B app = 20 T and a ramp- At the end of the FCM process (t=1250 s), the trapped field at the centre of the top surface (z=+0.5 mm), B t , of the inhomogeneous bulk is approximately 7.4 T, and the trapped field at the centre, B c , is approximately 11 T (see 7.6 T and 11.4 T, respectively, for the homogeneous bulk described in section 2.2.1). Figure 21 shows a 3D picture of the hoop stress, σ j , distributions for the homogeneous (left) and inhomogeneous (right) J c assumptions at the end of the FCM process and figure 22 shows a 2D surface plot of σ j across the centre of the bulk (xy-plane, where the c-axis is oriented along the zaxis), where the supposed GSBs and GSRs are highlighted by dotted and dashed white lines, respectively. It is clear that the inhomogeneous ab-plane J c distorts σ j with lower stress in low-J c regions and higher stress in high-J c regions. Figure 23 shows a 2D plot of the σ j distributions for the homogeneous and inhomogeneous J c assumptions at z=0 mm, i.e. across the centre of the bulk. The plots for the inhomogeneous bulk correspond to the supposed GSBs and GSRs highlighted in figure 22. This further emphasises the difference between the stresses developed in an inhomogeneous bulk compared to a homogeneous one and shows how the 3D model is useful for incorporating more realistic assumptions when a 2D axisymmetric model is insufficient. Figure 24 shows the 2D surface plot of σ j across the centre of the bulk (xy-plane) and figure 25 shows a 2D plot of σ j across the centre of the bulk at z=0 mm along the supposed GSBs and GSRs at the end of the FCM process, in the case that the bulk is reinforced by the stainless-steel ring. In order to better highlight the differences within the homogeneous and inhomogeneous bulks, the ring is not included in these plots. Here σ j includes both contributions from σ j (FCM), the hoop stress related to the Lorentz force, and σ j (COOL), the thermal compressive stress applied to the bulk when cooling to T op . The trends for σ j are identical to those shown in figures 22 and 23 for the unreinforced case; however, there is an offset owing to σ j (COOL), such that σ j is entirely negative (and less than 60 MPa). Thus, the 3D model provides a useful tool to analyse the influence of inhomogeneous material properties and design adequate reinforcement so that mechanical fracture is avoided in high magnetic fields.

Conclusion
In this paper, 2D axisymmetric and 3D finite-element models based on the H-formulation, implemented in the commercial finite element software package COMSOL Multiphysics, are used to analyse the mechanical stresses in bulk superconductor magnets during FCM with and without mechanical reinforcement using a stainless-steel ring.
1. Under identical magnetising conditions and assuming the same material properties, the 2D axisymmetric and 3D finite-element models produce consistent results: the 2D model provides a faster route to parametric optimisation, but the 3D model is required for bulk superconductor shapes that do not have cylindrical symmetry or for inhomogeneous properties in the Φdirection, i.e. around the ab-plane. The models can be used to examine the stresses developed during the magnetisation process, which depend on a number of important factors, including the J c (B, T) characteristics, sample geometry, ramp-down rate and so on. It is shown how the use of a mechanical reinforcement ring, with a higher coefficient of thermal expansion than the bulk superconductor, significantly reduces the mechanical stresses developed during FCM owing to σ j (COOL), the thermal compressive stress applied to the bulk by the ring when cooling to the operating temperature, T op . 2. The 2D axisymmetric model, assuming homogeneous properties around the ab-plane, is used as a fast optimisation tool to investigate the influence of any heat generated during the ramp down of the applied field by including the heat source, Q=E·J. A faster ramp-down rate of the applied field results in a larger temperature rise, resulting in an uneven temperature distribution within the bulk, depending on the method of cooling. This will distort the magnetic flux and current density distributions, which, in turn, affect the developed stresses. The presence of a reinforcement ring not only provides compression of the bulk during FCM to avoid mechanical fracture, but also provides an additional thermal pathway to remove the heat generated during the magnetisation process, which is particularly useful to consider in the case of PFM. 3. The 3D model is used to investigate the influence of inhomogeneous J c (B, T) properties around the ab-plane and it shown how the developed hoop stress is distorted, with lower stress in low-J c regions and higher stress in high-J c regions. When the bulk is reinforced with a stainless-steel ring, the trends for σ j are identical to those in the unreinforced case; however, there is an offset owing to σ j (COOL), such that σ j is entirely negative under the conditions analysed.
The modelling framework couples together electromagnetic, thermal and structural mechanics models to provide a complete and detailed picture of the mechanical stresses developed in a bulk superconductor during the FCM process and a useful and flexible design tool to adequately reinforce the bulk to avoid mechanical fracture at high magnetic fields, taking into account many of the practical situations faced when carrying out such high field experiments.    1. The electromagnetic model is implemented using the 'Magnetic Field Formulation' interface in COMSOL's AC/DC module. 2. In order to make use of the symmetry in the 3D model (see figure 1(b)), additional boundary conditions must also be set for the sides (yz-plane, x=0; xz-plane, y=0) and bottom plane of the entire geometry: (xyplane, z=0): this is achieved using the 'Magnetic Insulation' node for the former (setting the tangential component of the magnetic potential to zero) and the 'Perfect Magnetic Conductor' node for the latter (setting the tangential component of the magnetic field to zero). 3. The thermal model is implemented using COMSOL's 'Heat Transfer in Solids' interface and the temperature, T(t), is set by applying a constraint across the centre of the bulk in the z=0 plane. 4. The electromagnetic and thermal models are coupled with COMSOL's 'Solid Mechanics' interface in its Structural Mechanics module to include analyses of the mechanical stresses during FCM. A 'roller' constraint (equivalent to a symmetry condition) is applied to the bulk (and ring when applicable) in the z=0 plane, making the displacement in the direction normal to this plane zero. 5. The Lorentz force, F L =J×B, is implemented as a force per unit volume using the 'Body Load' node where F r =B z ·J j and F z =-B r ·J j (2D axisymmetric model) or F x =J y ·B z -J z ·B y , F y =J z ·B x -J x ·B z and F z =J x ·B y -J y B x (3D model). 6. A thermal coefficient of expansion is included for both the bulk and reinforcement ring (see table 1) using COMSOL's 'Thermal Expansion' multiphysics interface to simulate the thermal compressive stress applied    to the bulk by the reinforcement ring when cooled during FCM because of the difference in this coefficient between the two materials [9].