An accurate and efficient three‐dimensional high‐order finite element methodology for the simulation of magneto‐mechanical coupling in MRI scanners

Transient magnetic fields are generated by the gradient coils in an magnetic resonance imaging (MRI) scanner and induce eddy currents in their conducting components, which lead to vibrations, imaging artefacts, noise, and the dissipation of heat. Heat dissipation can boil off the helium used to cool the super conducting magnets and, if left unchecked, will lead to a magnet quench. Understanding the mechanisms involved in the generation of these vibrations, and the heat being deposited in the cryostat, are key for a successful MRI scanner design. This requires the solution of a coupled physics magneto‐mechanical problem, which will be addressed in this work. A novel computational methodology is proposed for the accurate simulation of the magneto‐mechanical problem using a Lagrangian approach, which, with a particular choice of linearisation, leads to a staggered scheme. This is discretised by high‐order finite elements leading to accurate solutions. We demonstrate the success of our scheme by applying it to realistic MRI scanner configurations.


FIGURE 1
Magnetic resonance imaging scanner: model Magnetom Vida 3 T. Image courtesy of Siemens Healthineers combinations of these, a gradient in any desired direction can be obtained. This is then used to construct an image of the body.
The interaction between the resulting magnetic fields and the conducting radiation shields gives rise to an important problem in MRI scanner design. The transient magnetic field generated by the gradient coils induces eddy currents in conducting components. These eddy currents, in turn, produce electromagnetic stresses, causing the radiation shields to deform and vibrate. [1][2][3] The vibrations are undesirable as they lead to a deterioration in image quality with image artefacts as well as the associated noise, which can cause patient discomfort. The eddy currents, on the other hand, lead to heat being dissipated and deposited in the cryostat, which can cause helium boil off and potentially result in a costly magnet quench. * Understanding the mechanisms involved in the generation of these vibrations and the heat being deposited in the cryostat is key for a successful MRI scanner design. This requires the solution of a coupled physics magneto-mechanical problem, which will be addressed in this work.
The latest developments in MRI design have resulted in this problem attracting even more interest. Recently, the Magnetom Terra 4 manufactured by Siemens Healthineers, has become the first 7 T MRI cleared for clinical use. 5 This improvement in field strength leads to a higher image resolution making easier the diagnosis of subtle pathologies which were until now difficult to identify, but it also brings the problem of larger eddy currents and stronger vibrations. Low helium magnets have also come into production in the last years. Typical MRI scanners require thousands of litres of helium to maintain the superconducting coils at a temperature of 4 K. The availability of helium is becoming challenging, leading to rising prices and a big impact in the production cost of MRI scanners. On the other hand, magnets using low helium technology use only around 1 − 20 litres to operate. 6,7 For such designs, the minimisation of the heat deposited into the cryostat becomes even more important, as a small amount of helium boil-off could potentially result in a magnet quench.
For these reasons, and considering the high cost of experiments, as well as the large number of iterations in the design process required to develop a new MRI scanner, manufacturers are interested in the development of accurate numerical simulation tools that are able to predict the complex coupling between the electromagnetic and mechanical fields. Attempts at modelling the magneto-mechanical effects with commercial software, focusing on the structural design of superconducting coils for high field MRI scanners, have been made. 8,9 However, the currently available commercial software tools capable of solving the coupled physics problem of interest, such as ANSYS, 10 NACS,11 or Opera FEA, 12 are not specifically designed for MRI scanners and do not offer tailored solution methodologies. In general, even with fine discretisations, they are unable to correctly predict the coupled physics in an MRI scanner with sufficient accuracy and, therefore, manufacturers are interested in specialist solutions.
The application of finite difference time domain methods for the prediction of induced eddy currents in MRI scanners, but neglecting mechanical effects, was studied in other works. [13][14][15] Methods for the rapid design of MRI coils based on three-dimensional (3D) electromagnetic simulations have also been considered. 16 In the work of Rausch et al, 17 a low-order temporal finite element scheme for the solution of 3D magneto-mechanical problems was presented. This *A quench refers to the sudden loss of superconductivity when its temperature is raised. was extended in the work of Rausch et al 18 to consider acoustic effects. In the previous work by our group, we considered solutions in an axisymmetric configuration, using a hp-finite element software, a stress tensor formulation and a novel linearised approach using an alternating current-direct current (AC-DC) splitting, which allows for its time harmonic solution in the frequency domain. 1,2 an Eulerian setting, where the assumption of small displacements and velocities was made. This work was then extended to consider the fully nonlinear problem in the time domain in the works of Bagwell et al. 3,19 The use of hp finite elements, 20,21 as opposed to low-order finite element discretisations, was employed as it is known to enhance the resolution of fine scale features such as small skin depths at high frequencies, 22 allow volumetric locking phenomena in nearly incompressible materials to be overcome, [23][24][25][26] and improve the resolution of acoustic wave propagation. 27 The drawbacks of our earlier axisymmetric approach include that it only allows the rotationally symmetric z gradient coil to be considered and that it assumes small velocities and accelerations. These drawbacks will be addressed in this work.
The main novelty of the work is the development of a new computational framework for the simulation of 3D coupled magneto-mechanical problems arising in MRI scanners. Our work draws on a new theoretical Lagrangian formulation for the treatment of coupled magneto-mechanical problems, which no longer requires the assumption of small velocities and accelerations. This article complements a forthcoming work, which will derive this approach, by providing a full computational treatment of the resulting transmission problem that is necessary to obtain numerical solutions. This includes the approximation of the discrete fields using high-order H 1 and H(curl) conforming finite elements 28 and the use of a regularised approach to the eddy current problem, 22,29 avoiding the need to solve a saddle point problem to ensure that the Coulomb gauge is satisfied. As first presented in the work of Bagwell et al, 2 the authors employed an AC-DC splitting of the fields, but, importantly, when this linearisation is applied to the aforementioned Lagrangian formulation, it results in a novel staggered computational scheme, which lends itself to a computationally efficient and accurate solution methodology to the magneto-mechanical problem. 30 We are then able to utilise efficient Jacobi preconditioners previously developed for the solution of single physics problems 22,29 for the rapid solution of large systems arising from the discretisation of the problem. The approach is applied to the solution of challenging academic and industrially relevant problems. In particular, the simulation of MRI configurations including longitudinal as well as transversal gradient coils is considered. This paper is organised as follows. In Section 2, we provide a brief discussion of the Lagrangian and Eulerian descriptions of the fields and describe the key assumptions and steps made in deriving the coupled transmission problem for the Lagrangian treatment of magneto-mechanical problems. In Section 3, the AC-DC splitting and linearisation of the problem is applied to the new Lagrangian formulation of the problem. This section also presents the weak form of the problem and the treatment of surface integrals. The discretisation of the system for both the DC and AC stages and the staggered solution strategy is discussed in Section 4. Finally, a series of numerical results for the solution of challenging decoupled and coupled problems is included in Section 5. The document finishes with concluding remarks in Section 6.

EULERIAN AND LAGRANGIAN APPROACHES
In this section, we briefly describe the general features of Eulerian and Lagrangian approaches and present the assumptions made in our new approach. A transmission problem, which is the starting point for our computation procedure, is also included.

Eulerian formulation
In the Eulerian formulation, we express quantities in terms of the spatial coordinates x for a fixed frame of reference that is not a function of t. See Figure 2. This frame is, of course, the usual one used for Maxwell's equations. We write E, H, D, and B for the electric field intensity, magnetic field intensity, electric flux density and magnetic flux density, respectively, and e for the volume charge density in this configuration. It is also generally accepted that for linear materials the electromagnetic constitutive laws hold in an Eulerian setting where is the permittivity, P is the polarisation, is the permeability, M is the magnetisation, and the subscripts 0 and r denote the free space and relative values, respectively.

FIGURE 2 General motion of a deformable body
The momentum balance equation, describing the mechanical deformations, can also be written in an Eulerian form, but it is more complicated than the simpler form usually employed in the Lagrangian frame. Even if the displacements u are small, the Lagrangian and Eulerian forms do not necessarily coincide unless the velocities are also small. Assuming they do coincide, then a linear elastic body with small deformations is described by the constitutive relationship where m (u) is the Cauchy stress tensor, (u) ∶= ((∇u + (∇u) T )∕2 the small strain tensor, and C = C ijk e i ⊗ e j ⊗ e k ⊗ e is a rank 4 constitutive tensor with entries with ij the Kronecker delta, e i the ith unit basis vector, and and G the Lamé parameters. In magneto-mechanics with simple materials, the body force can be written as the divergence of a Maxwell stress tensor and the description of this, in the Eulerian frame, coincides with its description in the Lagrangian frame for small displacements and velocities. This was the approach described in other works [1][2][3]19 for magneto-mechanical coupling, which we henceforth call our Eulerian approach.

Lagrangian formulation
In a total Lagrangian form, the governing equations are expressed in terms of the reference configuration X. In the current or deformed configuration, after a given time t, the particles will in general have different coordinates x. The motion can be therefore expressed by means of a mapping as x = (X, t); see Figure 2. The Lagrangian form of Maxwell system is similar to that of the Eulerian form, except that the differential operators act with respect to X. In this configuration, the physical variables are given a subscript 0 and, hence, we write E 0 , H 0 , D 0 , and B 0 for the electric field intensity, magnetic field intensity, electric flux density, and magnetic flux density, respectively. In addition, we set e 0 as the volume charge density and J 0 to be the solenoidal current source. Importantly, J 0 has to be defined in the reference configuration as does the ohmic current E 0 , where is the electrical conductivity. When Maxwell's equations are transformed from the Lagrangian frame to the Eulerian frame, E 0 transforms to become the sum of an ohmic-type current J o = E and a Lorentz current B × v in the Eulerian form of Maxwell's equations. On the other hand, the constitutive laws (1) must be transformed from the Eulerian form to the reference frame for the total Lagrangian form.
The momentum balance equation in a total Lagrangian frame is expressed in terms of the displacements u, the first Piola-Kirchhoff tensor  , the body force f 0 , and the accelerations a = v t = 2 u t 2 . It takes a simple form with differential operators acting with respect to X. Note that the form of  depends on the constitutive model, it is certainly a function of the deformation, and it can be a function of D 0 and B 0 as well as other fields in general. [31][32][33] In the updated Lagrangian form, the Maxwell and linear momentum equations are expressed in terms of the current coordinates x = (X, t). In this configuration, the particles X acquire the temporary coordinate labels x = (X, t) at time t. The form of both the Maxwell system and the linear momentum equations in terms of these transformed variables is simple.
Although displacements are generally small in MRI scanners, the same is not necessarily true of velocities (especially in the case of stronger coupling). For small deformations, but not necessarily small velocities and accelerations, the total Lagrangian and updated Lagrangian descriptions coincide and, henceforth, we simply call this our Lagrangian approach. This approach will be pursed in this paper and is based on the following key assumptions and steps.
1. The eddy current approximation cannot be applied in the Lagrangian frame and instead must be applied in the Eulerian frame. After applying the eddy current approximation, the simplified Maxwell equations are then transformed to the Lagrangian frame, but without assuming small displacements. 2. The electromagnetic constitutive laws (1) hold in an Eulerian frame and are transformed to the Lagrangian frame. 3. In the Lagrangian frame, the displacements are assumed to be small so that total and updated Lagrangian descriptions coincide. This also means that elastic bodies have the constitutive behaviour m (u) ∶= C ∶ (u) and the body force can be rigorously described in terms of the divergence of a Maxwell stress tensor in terms of B 0 and H 0 . 4. Gauging is applied to the Lagrangian eddy current equations obtained from the transformation in 1 and employing the constitute behaviour in 2 assuming the displacements (but not the velocities) are small in 3. 5. For simplicity, we assume Biot-Savart coils so that the support of the current sources are not treated as conductors.
However, in general, this choice is not a limitation of our approach.
Applying this approach to a magneto-mechanical problem set on an unbounded domain n=1 Ω C,n , denotes the union of N disjoint elastic conducting bodies each with ≠ 0 and possibly different from 0 and Ω c C ∶= R 3 ⧵ Ω C denotes the nonconducting region ( Figure 3), and then we show in a forthcoming work that this problem can be described by the transmission problem. Find (A, u)(t) ∈ (R 3 × R 3 )(0, T] such that curl where I is the rank 2 identity tensor, [·] Ω C denotes the jump over Ω C , and we have chosen to set the initial conditions for the fields to be zero, corresponding to a system at rest at t = 0. The physical electric and magnetic fields are coupled to A and v through which can be applied once (4) is solved.

AC-DC SPLITTING AND PROBLEM LINEARISATION
For our desired application, we can decompose the current source as J 0 (t) = J DC +J AC (t) (Figure 4), where J DC corresponds to the static current source of the main magnetic coils and J AC (t) to the transient current source of the gradient coils. 2,3 Limiting ourselves to linear materials, such that ≠ (A), this decomposition allows us to introduce the following static problem. curl

Variational formulation 3.1.1 DC stage
To present the variational formulation of the transmission problem, we define The variational formulation of the DC stage then reads as follows: To overcome the requirement for weak solutions to satisfy ∇ · A DC = 0 in R 3 , we defineX and look for solutions.
In the aforementioned statement, is a small regularisation parameter that has been introduced to avoid imposing the Coulomb gauge and offers considerable computational advantages. This approach is discussed in detail in the works of Ledger and Zaglmayr 22 and Zaglmayr 29 ; it is shown that where * denotes the dual space. Note that, in the DC stage, we have a one way coupling, as A DC is independent of u DC , but u DC depends on A DC . This allows us to solve the system in a staggered manner, as opposed to the approach presented in the work of Bagwell et al, 3 which involved linearisation and then application of a monolithic Newton-Raphson solver. In this new approach, once discretised, we first solve (8a) followed by (8b). The discretisation process will be described in Section 4. (4), which satisfy the residuals

AC stage
for all real valued ( A, u) ∈ (X R 3 × Y ( )). The residuals vanish when (A, u)(t) are the true solutions. Following the same procedure as in the Eulerian approach (see the works of Bagwell et al 2,3 ), the problem can be linearised about the static solutions leading to the system that is linear in the time dependent updates ΔA and Δu. The definitions of the directional derivatives are as stated in the work of Bagwell et al. 2 This motivates the time harmonic representation of the fields, ie, where i ∶= √ −1 and = 2 f represents the angular frequency of the driving current in the gradient coils in the case of a harmonic excitation. In addition, A A A AC ∶= ΔA A A, U U U AC ∶= ΔU U U and J J J AC represent the complex amplitudes of the corresponding time-varying (update) fields ΔA, Δu, and J AC , respectively. The total time dependent fields can be recovered as Then, the problem can be formulated in the time harmonic domain using the regularised approach as follows: Find In the aforementioned statement, the directional derivatives are defined as is the linearised Maxwell stress tensor. Note that, as opposed to our previous Eulerian approach, 2,3,19 vanishes in the Lagrangian approach. Thus, once discretised, the problem can be solved in a staggered manner by solving (13a) followed by (13b). Furthermore, the time harmonic physical electric and magnetic fields for the AC stage are defined by the linearised and regularised version of (5) where B B B AC 0 = curlA A A AC and B DC 0 = curlA DC and the complete physical fields for both stages are

Surface integral treatment
In both the DC and AC stages, surface integrals are present, which involve e and −1 S, respectively. In the case of the former, integration by parts yields In the case of the AC stage, it can be shown that div ) (Appendix) and, thus, we also have for all real valued u ∈ Z since the current sources are away from the interface. When we replace A DC by A DC and A A A AC by A A A AC , a small approximation error is introduced and we assume the effect on the aforementioned results is negligible. These results will be of practical value for our discrete implementation presented in the next section.

DISCRETE SYSTEM
In this section, we briefly describe the discretisation of the weak forms for the DC and AC stages of the Lagrangian approach. To do this, we replace the bounded domain with a finite computational domain Ω = Ω C ∪ Ω NC , which is truncated sufficiently far from the object. As an approximation to the static decay of A DC and A A A AC in Ω c C , we set n×A DC = and n × A A A AC = 0 on Ω and ensure that Ω is sufficiently far from the object. Throughout, we focus on the partition of Ω by unstructured tetrahedral elements, due to the availability of automatic unstructured meshing algorithms for generating meshes around complex configurations, and, here, we employ the NETGEN mesh generator 35 for this purpose. To illustrate the discretisation process, we begin with the DC stage and then also summarise the results for the AC stage.

DC stage 4.1.1 Electromagnetics
We introduce a nonoverlapping unstructured tetrahedral partition of the domain Ω = where Ω (e) is the region corresponding to tetrahedral element (e). Then, we introduce a discrete Galerkin approximation and look for solutions This approach greatly reduces the dimension of the resulting system. 22 Here, P global is the number of electromagnetic degrees of freedom in the problem and N 1 , N 2 , … , N P global are H(curl)-conforming (see the work of Nédélec 36 )-type basis functions, for which we choose the particular set introduced by Schöberl and Zaglmayr 28 and Zaglmayr. 29 Their construction is hierarchic and the lowest order q = 0 elements are often called edge elements, since, for q = 0, there is one degree of freedom (DOF) per edge (six DOF per tetrahedron) and, for higher order discretisations (q ⩾ 1), there are additional degrees of freedom associated with edges, faces, and element's interiors. For further details of the implementation, see other works. 22,28,29,37

Mechanics
In a similar manner to Section 4.1.1, a nonoverlapping unstructured tetrahedral partition of the domain Ω C = ⋃ E m e=1 Ω (e) ⊂ Ω is introduced. A discrete Galerkin approximation is employed and we look for and L i ∈ H 1 (Ω C ) are the hierarchic H 1 -conforming basis functions proposed by Schöberl and Zaglmayr. 28 The basis functions are scalars, as opposed to the vectorial H(curl)-conforming basis functions, and Q global is the number of mechanical degrees of freedom in our problem. The order of the mechanical approximation is denoted by p ⩾ 1 and the mechanical degrees of freedom comprise of low-order vertex and high-order edge, face, and interior DOFs. The p = 1 basis functions correspond to the standard linear finite element hat functions. For further details of the implementation, see other works. 28,29,37

System
By assembling blocks, the discrete systems are as follows: Find A DC such that and then u DC such that Note that, here, and in the following, we use Roman fonts for matrices. The different blocks of the system are where r = ∕ 0 , 0 = 4 × 10 7 H/m, m (u DC ) = C ∶ (u DC ) and the symmetry of m (u DC ) has been used. The multiplication by 0 is used to improve the scaling of the equations. The matrix K DC AA + C DC AA is real and positive semidefinite and it is of the form obtained in magnetostatic problems discretised by the hp-finite element basis of Schöberl and Zaglmayr. 28 By exploiting the low-order-edge high-order-edge-face-interior splitting simple block Jacobi preconditioning with conjugate gradients becomes effective 22,29 for the solution of (20). In the case of (21), the matrix K DC uu is real positive definite and is of the form obtained in elasticity problems discretised by the hp-finite element basis of Schöberl and Zaglmayr. 28 The block Jacobi preconditioner which exploits the low-order-vertex high-order-edge-face-interior splitting of the basis functions is used in combination with conjugate gradients for its solution. In the aforementioned statement, the superindices V, E, F, and I denote the vertex, edge, face, and interior degrees of freedom, respectively, and the tilde is used to indicate that the block is block diagonal, ie, only the degrees of freedom associated to a single entity (edge, face, or cell), are used.

AC stage
The details of the discretisation process for the AC stage is similar to that described in Section 4.1; for further details, we refer to the work of Seoane. 37

System
Mixed approaches arise in the discretisations of linear elasticity problems to avoid mechanical locking, 23,24 Navier-Stokes equations for describing the pressure and velocity fields, 38 and in the discretisation of the Maxwell system if a Lagrange multiplier is used to enforce the divergence constraint. 39 In each case, they lead to systems that combine elements of different types and different orders. In the latter two cases, a saddle point problem results and, in such cases, it is not only important that the correct discrete spaces are chosen, but also that the order of the elements be chosen with care so as to ensure that the (discrete) LBB condition is satisfied. Following similar steps to those described in Section 4.1, we approximate U U U AC by U U U AC hp and A A A AC by A A A AC ,hq using H 1 and H(curl) conforming elements, respectively, and obtain the discrete dynamic system in the Lagrangian approach in the form where different blocks of the system can be shown to be 37 [ These can be expressed in terms of elemental contributions and assembled to form the global matrices. The system (24) combines two element types, but in common with mixed approaches to elasticity with ≠ 1∕2, it does not represent a saddle point problem and, therefore, we can be flexible with the degree approximation p for U U U AC hp and q for A A A AC ,hq .

Staggered algorithm combining AC and DC stages
Unlike the Eulerian approach, 2,3 the solution vector  AC can be computed independently of  AC in the AC stage and we propose the simple staggered solver Algorithm 1, which combines the DC stage (20) and (21) and the AC stage (24) as well as the physical field representations (16). Notice that steps 1, 2, 4, and 5 are simple linear algebra steps and steps 3 and 6 follow from using the discrete representation of the fields using the coefficients computed from the solution of the linear systems. This is further emphasised by the different choice of fonts in steps 3 and 6 compared to the others.
Compared to the Eulerian approach, Algorithm 1 offers considerable advantages when applied to large 3D problems since, in the AC stage, rather than the solution of a single complex indefinite system, the problem has been reduced to the solution of two smaller systems, the larger of which can be solved by preconditioned iterative approach. Specifically, the matrix K AC  + C AC  is complex symmetric and is positive semidefinite. It is of the same form as obtained in the discretisation of eddy current problems by hp-finite elements. 22 Therefore, we apply the same preconditioned GMRES technique as described in this work for computing  AC . The matrix K AC   − 2 M AC   is real and indefinite and it is of the form obtained in the discretisation of elastic wave propagation by hp-finite elements. For a sufficiently fine grid, an iterative approach to its solution could be applied, but, given its relatively small dimension compared to K AC  + C AC  ; we instead employ a direct solver for the computation of  AC , which is less restrictive.
The AC stage of the Eulerian approach, which was solved monolithically in the works of Bagwell et al, 2,3 can be solved instead using a fixed point strategy. 37 A flow chart comparing the solver structure for the Lagrangian and Eulerian formulations is provided in Figure 5, where the advantages of the Lagrangian formulation can be clearly identified.

Mechanical damping
In practical MRI scanners, the conducting components are connected together by materials that absorb and dampen the vibrations. These effects are not taken in to account in the mathematical model described previously and, in general, the amount of damping is difficult to quantify precisely. As a first approximation to damping effects, the simple case of (A) (B) proportional damping is considered. In this model, a damping block C AC (24), where M is the mass proportional damping coefficient, which can be determined from 40 where is a dimensionless factor, called the damping ratio, which is used to attenuate the frequency . The damping ratio describes how the amplitude of the oscillations in a system decay to an equilibrium position after a perturbation. According to the damping ratio, a system can be undamped ( = 0), underdamped ( < 1), critically damped ( = 1), and overdamped ( > 1). If the system is critically damped, it will decay to equilibrium without oscillation in the fastest admissible way. In an overdamped system, an equilibrium position will be reached without oscillation, but in a longer time than in the case of critical damping. If the system is underdamped, it will oscillate to reach equilibrium. The cases of an undamped and underdamped system will be considered in the next sections.

Decoupled physics
A series of examples are presented to illustrate the performance of the DC and AC stages of our scheme for decoupled physics problems. For further examples, see the work of Seoane. 37 We will measure the accuracy of the approaches using the error measures

DC stage, elasticity: hollow cylinder subject to pressure field
In order to benchmark the DC (elasticity) stage, we consider the hollow cylinder Ω C ∶= {(r, , z) ∶ a ⩽ r ⩽ b, 0 ⩽ < 2 , −L∕2 ⩽ z ⩽ L∕2} subject to internal and external pressure. A cross-section of the scheme of the problem is presented  Figure 6. This problem has an analytical solution in terms of u exact , which can be found in the work of Bower. 41 The particular case with internal radius a = 1 m, external radius b = 2 m, length L = 5 m, Young modulus E = 1 × 10 3 Pa, and Poisson's ratio = 0.33 is considered here. The internal and external pressure are set to p a = 10 4 Pa and p b = 10 7 Pa, respectively. This problem is simulated by solving a simplified version of (21) using appropriate boundary conditions (28b) Figure 7 shows where the boundary conditions are applied and one of the meshes used for the computations. We consider a sequence of unstructured tetrahedral meshes with 1366, 3144, 7953, 29 637, and 53 230 elements and discretisations with uniform element orders p = 1, 2, 3, 4, 5 applied, in turn. Quadratic representation of the geometry is used throughout. In each case, the relative error ||u DC hp − u exact || H 1 (Ω C ) ∕||u exact || H 1 (Ω C ) is measured. Figure 8 shows this relative error against the number of degrees of freedom for p-refinement, where each line represents a single mesh and the points increasing polynomial degree, and Figure 9 those for h-refinement, where each curve is for a fixed order and the points represent mesh refinement.
In Figure 8, we can observe the exponential rate of convergence of the relative error with the number of degrees of freedom under a p-refinement analysis. This trend is confirmed by the algebraic rate of convergence of the error, which is obtained when the error is plotted against the number of degrees of freedom to the power of 1∕3, and agrees with the predicted asymptotic rate. 42 After an initial preasymptotic region, Figure 9 shows the algebraic rate of convergence expected for the h-refinement.
Despite the coarse nature of some of the grids employed, accurate representation of the solutions can still be achieved with high-order elements. For the purpose of plotting, coarse meshes are subdivided and solutions evaluated on the split mesh and then passed to ParaView 43 for plotting. Figure 10 shows contour plots of the displacements in the cylinder obtained by a mesh of 7953 tetrahedra of order p = 4. A, f = 50Hz; B, f = 1000Hz

AC stage, eddy current: conducting sphere in a uniform magnetic field
In order to benchmark the AC (eddy current) stage, we consider the conducting sphere Ω C ∶= {x ∶ |x| ⩽ R} placed in an unbounded region of free space Ω c C , where H H H AC 0 tends to a time harmonic uniform field far from the object. An illustration of the problem is presented in Figure 11. This problem has an analytical solution 44 in terms of A A A exact . The case of the sphere with R = 1 cm, = 1.5 0 , 0 = 4 × 10 −7 H/m, and = 6 × 10 6 S/m set in a background field with lim |x|→∞ H H H AC 0 = 0 e z A/m and a frequency f = 50 Hz is considered, unless otherwise stated. This problem is simulated by solving a simplified version of (24) with J J J AC = on a truncated domain where Ω is the surface of the sphere of radius 2 cm. On the boundary Ω, the condition n × A A A AC ,hq = n × A A A exact , where A A A exact is the analytical solution for this problem proposed by Smythe, 44 is imposed. The regularisation parameter = 10 −5 is applied in all cases. In Figure 12, we show the convergence of the preconditioned GMRES solver for two different frequencies, ie, f = 50 Hz and f = 1000 Hz, when a p-refinement analysis is performed using a mesh of 2319 tetrahedral elements. The chosen tolerance for the iterative solver is set as 10 −5 . It can be observed that, for the worst case, corresponding to f = 1000 Hz and q = 4, only 120 iterations are required, which shows the efficiency of the preconditioning technique for the computation of the solution vector A A A AC .
To investigate the accuracy of the approximation, we consider a sequence of unstructured tetrahedral meshes with 195, 313, 509, 1083, 6654, 9957, and 19 791 elements and use discretisations with uniform element orders q = 0, 1, 2, 3, 4 applied, in turn. Initially, we set the geometry to be approximated by polynomials of degree g = 4. In each case, the relative error is measured. Figure 13 shows this relative error against the number of degrees of freedom for p-refinement and Figure 14 those for h-refinement. In a similar manner to the DC (elasticity) stage example, we observe exponential rate of convergence of the relative error with the number of degrees of freedom under a p-refinement analysis in Figure 13. In addition, after an initial preasymptotic region, Figure 14 shows the algebraic rate of convergence expected for the h-refinement.
We have investigated the role played by the geometry approximation under p-refinement. To do this, we fix a mesh of 195 tetrahedra and consider the accuracy of the solution for geometry approximations of degrees g = 0, 1, 2, 3, 4 in Figure 15. In this figure, we can observe the need for accurate geometry resolution when computing solutions on coarse grids.
The skin depth is a measure of the penetration of the fields into the conducting object, and is defined as the distance at which the fields decay to 1∕e of its value in the conductor's surface. 44 Therefore, for higher frequencies (and hence smaller skin depths), we will have larger gradients in field close to the conductor's surface, which makes the accurate numerical simulation of the problem more challenging. This effect is illustrated in Figure 16, where we show a contour plot of the computed eddy currents J J J o hq = −i A A A AC ,hq in the sphere, in particular, in the planes z = 0 and x = 0. To obtain these results, we have employed sufficiently fine discretisations to capture the solutions accurately.

Coupled physics
We now present a series of more challenging coupled examples to illustrate the performance of the staggered Lagrangian approach advocated in Algorithm 1 compared to the Eulerian approach using a fixed point algorithm described in the work of Seoane. 37

Conducting and elastic sphere in a background magnetic field
The case of a conducting elastic sphere Ω C ∶= {x ∶ |x| ⩽ R}, placed in an unbounded region of free space Ω c C , where H H H AC 0 tends to a time harmonic uniform field far from the object, is now considered. Unlike the previous two examples, this does not have an analytical solution. The sphere has a radius R = 1 cm and material properties = 0 , = 6 × 10 7 S/m, = 7800 kg/m 3 , E = 10 8 Pa, and = 0.3. For computational purposes, the domain was truncated at a finite distance from the conducting sphere. In particular, the truncation was chosen to be at a radius of R 0 = 2 cm from the centre of the sphere. The geometry was further simplified by exploiting the symmetry of the problem about e z , which reduces Ω to a hemisphere with outer radius R 0 . An illustration of the reduced domain is provided in Figure 17. On the truncation boundary, the approximate condition n×A A A AC ,hq = n×A A A exact , where A A A exact is the analytical solution for the corresponding nonelastic sphere proposed by Smythe,44 as also used in Section 5.1.2, is imposed. On the symmetry boundaries, the conditions are imposed, which follow from the symmetric nature of the problem. The regularisation parameter was chosen to be = 10 −6 . As discussed in Section 4.2.1, we are free to choose different combinations for the element orders p and q provided that the discretisation for each field is sufficient to resolve accurately both the electromagnetic and mechanical fields. To illustrate this, a frequency f = 50 Hz is considered, and the aforementioned sphere is discretised by a mesh of 2319 tetrahedral elements. The converged results for the order combinations p = q, p = q − 1, and p = q + 1 are then studied.
In particular, Figure 18 Figure 19 for the aforementioned mesh and elements of order p = 4 and q = 5 for all frequencies in the sweep. This combination of p, q is based on the observation that higher order was needed to accurately resolve the eddy currents compared to the displacements at high frequencies. At lower frequencies, it represents an overkill solution. The order of approximation for the geometry was set to g = 2. The result obtained with the Eulerian formulation is also included to compare both formulations. Note that, for the Eulerian approach, a small mechanical damping has been applied to the system to remove numerical artefacts leading to nonphysical peaks. To apply the damping, we select f = ∕2 = 2960 Hz, corresponding to the resonant peak, and apply (26) with = 5.3 × 10 −3 to obtain M = 200, which we use for all frequencies in the sweep. The axisymmetric software, 2,3 using a monolithic instead of a fixed point solver for the Eulerian formulation, is also shown as a further means of comparison. This software being applicable, in this case, is due to the symmetry of the problem. Furthermore, the previous Eulerian axisymmetric solver has also been adapted to perform the Lagrangian approach, described by Algorithm 1, and these results are included in the figure. For both the Eulerian and Lagrangian approaches, the axisymmetric software is in good agreement with the results of the new 3D solver. In this problem, and for the considered frequency range, no significant differences between both approaches are observed.
(A) (B)  In order to better illustrate the physical behaviour of the problem, a contour plot of |Re(J J J o hpq )|in the conducting sphere and |Re(B B B AC 0,hq )| in the surrounding air is presented in Figure 20. Here, as in (30a), we use J J J o = E E E AC , approximated as J J J o hpq = (−i A A A AC ,hq +i B DC 0,hq ×U U U AC hp ), corresponding to ohmic-type currents in an Eulerian frame. Streamlines of the magnetic flux density in the air are also included.
Finally, as a further illustration of the coupled physical behaviour, the deformed sphere for a set of different frequencies f = {2000, 2500, 3000, 3500, 4000} Hz is shown in Figure 21.

Test magnet problem with z (longitudinal) gradient coil
The second benchmark problem considered consists of an MRI configuration, where the geometry has been simplified so that it is rotationally symmetric. In addition, the non-rotationally symmetric x and y gradient coils have been removed leaving only the z gradient coil. The configuration of the problem is shown in Figure 22, where the main coils are shown in red, the gradient coils in blue, and the different radiation shields (4K, 77K, and outer vacuum chamber (OVC)) in different green tones. Note that the 3D geometry ( Figure 22B), can be obtained by rotating the r − z plane ( Figure 22A), through 0 ⩽ < 2 rad about the z axis (where (r, , z) are cylindrical coordinates). For further details about this problem, we refer to the work of Bagwell et al 2 where this problem has been previously simulated using an axisymmetric software and the Eulerian approach. In this section, the focus is the 3D simulation of this problem using the Lagrangian approach.
The domain Ω is formed by truncation of Ω c C by a cylinder of radius R 0 = 0.9 m and length L 0 = 2.4 m centred about the magnet. This problem is axisymmetric, but in order to illustrate the performance of our 3D solver, we make an arbitrary choice such that Ω is reduced to one quarter of the original geometry, as shown in Figure 23, and apply the symmetry boundary conditions (29) on the symmetry boundaries. Algorithm 1 is then applied to obtain the computational Lagrangian solution.
A frequency sweep from 1 to 5000 Hz is performed in order to observe the variation of the output power . For this, a mesh of 33 805 tetrahedral elements is generated and different polynomial orders q and p = q + 1 applied. This mesh was generated by locally refining the elements in the conducting shields to accurately resolve the fields in this region and paying attention to avoid issues with the generation of highly stretched tetrahedrons. The case of a system without damping is analysed first, followed by a damped system according to the approach in Section 4.4. In the latter, = 4.5 × 10 −4 corresponds to the first resonant frequency f = 3515 Hz of the system. In a similar way to Section 5.2.1, by applying (26) with = 4.5 × 10 −4 at this resonant peak, we obtain M = 20, which we apply for all frequencies in the sweep. Figure 24 shows the convergence of the dissipated power and kinetic energy under p-refinement for the undamped system. We observe the solution is effectively converged using q = 3 elements, but there are still some differences between the results with q = 3 and q = 4 in the 4K shield. Nevertheless, as it will be shown later, the result achieved with q = 4 is in perfect agreement with the converged result in the case of a system with damping and that of the axisymmetric software using the Lagrangian approach.
The simulations were then repeated for a system with damping leading to Figure 25 for the p-refinement analysis. Notice how the addition of a small amount of damping leads to a faster convergence of the frequency sweeps under p-refinement compared to the undamped case.
Next, in Figure 26, we compare the converged 3D results with those obtained with the axisymmetric code employing both the Eulerian and the Lagrangian approaches as well as with the commercial software NACS 11 (only for the dissipated power), which uses low-order elements. We observe excellent agreement between the 3D and axisymmetric results for the Lagrangian approach. The axisymmetric results using the Eulerian approach are in good agreement with the Lagrangian approach for the OVC shield and for the 77K shield in the case of small frequencies. However, significant differences between the two approaches are observed in the 77K shield for high frequencies and in the 4K shield for almost all frequencies. The reason for these differences can be attributed to the stronger coupling and larger velocities and accelerations present in the shields at these frequencies, which mean that the approximations in our previous Eulerian approach break down. The results obtained with NACS software are in good agreement with the Lagrangian approach in the case of the OVC shield, but differ significantly in the 77K and 4K shields. It can be conjectured that this is due to the inability of low-order elements to accurately resolve the complex coupling mechanisms involved in the problem. Note that the 77K and 4K shields have a bigger conductivity and, therefore, smaller skin depths compared to the OVC shield.
Finally, it can be observed that the resonant peaks in the dissipated power are at similar frequencies to those in the kinetic energy. This is due to the resonant modes of the mechanical problem being driven by the generated eddy currents in the shields. In terms of computational cost, the Lagrangian approach offers considerable advantages over the Eulerian approach since it requires the solution of two small systems as opposed to one larger monolithic system (eg, for p = 4, q = 5 P global = 1 277 370 and Q global = 751 245 compared to an Eulerian, which would involve a P global + Q global sized system for each frequency). Furthermore, as described in Section 4.2.1, the matrices in the Lagrangian approach have nicer properties compared to the indefinite system that would need to be solved for the Eulerian approach, which allow preconditioned iterative solvers to be applied.
In order to illustrate the physical behaviour of the problem, the eddy currents in the deformed 77K shield are shown in Figure 27 together with the streamlines of Re(B B B AC 0,hq ) and contours of its magnitude in the air. These results are shown for two different frequencies, ie, f = 10 and f = 1000 Hz, to show how the magnitude of the eddy currents increases and the skin depth effect becomes visible for the higher frequency. It can also be observed that the deformation increases with the eddy currents.
Finally, in Figure 28, we show a contour plot of |Re(U U U AC hp )| in the deformed 4K shield together with the magnetic flux density streamlines in the surrounding air. The result is shown for a frequency of f = 1000 Hz and the different snapshots correspond to different times. Note that, for a given frequency, the temporal solution can be recovered from (12).

Test magnet problem with x (transversal) gradient coil
A new challenging problem, consisting of the configuration described in Section 5.2.2 with the z gradient coils being replaced by a set of x gradient coils, is now considered. Note that the case of only y gradient coils can be obtained by rotating the current situation through 90 • . This new set of coils produces a magnetic field gradient in the x direction and, therefore, the problem is no longer rotationally symmetric and cannot be simulated using our previous axisymmetric simulation tool. Figure 29 shows the geometry of the problem, including the new set of gradient coils. These coils represent an approximation to the complex fingerprint shape shown in Figure 30 used in real MRI scanners. A similar approximation was also used in the case of the longitudinal coil discussed in Section 5.2.2.
In the case of transversal gradient coils, the description of the current source J J J AC becomes more challenging due to increasing complexity in the geometry. This can lead to issues with producing a description of J J J AC that is both tangential to the coil geometry and that satisfies divJ J J AC = 0. In order to ensure that the current source is divergence free, a mapping of the current source to the space of divergence free functions is performed. To do this, we write J J J AC as the curl of a vector potential and solve a problem analogous to that for A DC . For details, we refer to the work of Seoane 37 where this is discussed in detail.
In the new configuration, the symmetry with respect to the x axis is lost. Therefore, the computational domain is chosen as half of the full domain instead of a quarter, which increases the computational cost significantly. For the discretisation, a mesh of 92 434 tetrahedral elements was used and the polynomial order was increased until convergence was reached, which was found to when q = 4 and p = 3 for the chosen frequency of f = 500 Hz. In order to illustrate how the solution changes when considering transversal instead of longitudinal coils, a contour plot of |Re(J J J o hpq )| in the OVC shield is presented in Figure 31 comparing both cases. The frequency of f = 500 Hz was chosen as this represents already a strong coupling between physics and the numerical computations are not as expensive as in the case of higher frequencies. As expected, it can be observed that the solution in the case of x gradient coil is no longer axisymmetric. The same is repeated for the mechanical field in Figure 32, where we show a contour plot of |Re(U U U AC hp )| on the deformed 77K shield.   Finally, in Figure 33, the streamlines of Re(B B B AC 0,hq ) around the radiation shields are shown, together with contours of |Re(J J J o hpq )| in the shields. The main coils and gradient coils are also included and it can be observed that, in the case of z gradient coils, all the streamlines lie in a plane, whilst in the case of x gradient coils, they follow a more complex helical pattern.
Although transversal gradient coils result in the loss of rotational symmetry of the problem, some symmetries and antisymmetries with respect to the Cartesian planes are still present. In future research, antisymmetries will be exploited to further reduce the cost of the simulations when considering transversal coils.

CONCLUSIONS
In this paper, a new method for the solution of 3D magneto-mechanical problems arising in MRI scanner design has been presented. This is based in a new Lagrangian approach for the description of the coupled physics and, in a forthcoming work, we will also provide a rigorous derivation of this approach. The electromagnetic and mechanical fields are discretised using high-order H(curl) and H 1 conforming finite elements ensuring the continuity requirements for each of the fields are met and the potentially small skin depths are accurately resolved. This greatly improves our previous work, 1,2 which was limited to the case of axisymmetric configurations and small velocities and accelerations. Furthermore, the new Lagrangian approach results in a staggered solver, which presents several computational advantages compared to the fixed point used in our previous Eulerian approach. Single physics problems were used to prove the optimal convergence of the software. Finally, this tool was applied to the solution of industrially relevant MRI models, including the challenging case of transversal gradient coils. Comparisons with the previous axisymmetric software were made for the case of longitudinal gradient coils to validate the new 3D solver. The next steps of our research will focus on the simulation of more complex MRI configurations. where B B B AC 0 ∶= curl A A A AC . For simplicity, we drop the subscript 0 in the following, then, using index notation, this can be written as Thus, assuming is constant in the region of consideration, which will be the case in Ω c C , then (div −1 S) i where, in the final step, the last two terms vanish due to div B DC 0 = 0 and divB B B AC 0 = 0 in R 3 . Next, using the property that in Ω c C ⧵ (supp(J J J AC ) ∪ supp(J DC )).