Transient solutions to nonlinear acousto‐magneto‐mechanical coupling for axisymmetric MRI scanner design

In this work, we simulate the coupled physics describing a magnetic resonance imaging (MRI) scanner by using a higher‐order finite element discretisation and a Newton‐Raphson algorithm. To apply the latter, a linearisation of the nonlinear system of equations is necessary, and we consider two alternative approaches. In the first approach, ie, the nonlinear approach, there is no approximation from a physical standpoint, and the linearisation is performed about the current solution. In the second approach, ie, the linearised approach, we realise that the MRI problem can be described by small dynamic fluctuations about a dominant static solution and linearise about the latter. The linearised approach permits solutions in the frequency domain and provides a computationally efficient way to solve this challenging problem, as it allows the tangent stiffness matrix to be inverted independently of time or frequency. We focus on transient solutions to the coupled system of equations and address the following two important questions: (i) how good is the agreement between the computationally efficient linearised approach compared with the intensive nonlinear approach and (ii) over what range of MRI operating conditions can the linearised approach be expected to provide acceptable results for outputs of interest in an industrial context for MRI scanner design? We include a set of academic and industrially relevant examples to benchmark and illustrate our approach.


INTRODUCTION
In recent years, magnetic resonance imaging (MRI) scanners, illustrated in Figure 1, have become an indispensable tool for use in medical imaging. Their capability to nonintrusively diagnose a range of medical ailments, such as tumours, 1 damaged cartilage, fractures, 2 internal bleeding, and even detection of multiple sclerosis 3 make them a desirable imaging acoustic effects in the human head. 30 More recently, a few attempts to analyse the magneto-mechanical coupling in MRI systems have been considered; with the modelling of axisymmetric superconducting solenoids in self magnetic fields 31 and efficient low-order FE solvers for magneto-mechanical coupling in the work of Rausch et al, 32 which utilised calculation of body forces to generate a weakly coupled algorithm. This work was later extended to include acoustic effects in the aforementioned authors' other work. 33 Our previous work 34 focused on the solution of coupled magneto-mechanical behaviour of geometrically idealised (axisymmetric) MRI scanners, with arbitrarily high-order finite elements, 35,36 through a stress tensor approach to avoid direct calculation of the electromagnetic body forces. Given the relatively small-scale phenomena that emerge in such magnetic and acoustic applications, such as small skin depths in conducting components (known as skin effect) 37 and high frequency acoustic waves, high-order elements are necessary to accurately resolve such effects and ensure accurate solutions, as shown in our other work. 13 Furthermore, by choosing to linearise the transient nonlinear system of coupled equations about the static fields, the formulation simplifies greatly, resulting in a transient linear system of equations dependent on the transient magnetic field excitation, as shown in our other work. 13 Similar techniques, involving the additive split of a nonlinear problem to a series of linear problems, have been successfully applied to the field of computational mechanics, such as analysis of structural membranes, [38][39][40] high-order mesh generation, 41 and in biomedical applications. 42 This formulation permits a simple linearised approach that can be solved in the frequency domain and, given the relatively small frequency ranges of MRI applications, provides rapid solutions to industrial problems. However, the applicability of this linearisation will depend on the strength of the coupling between the acoustic, mechanical, and electromagnetic fields. In the context of MRI scanners, it is thus imperative to answer two important questions: i) how good is the agreement between the computationally efficient linearised approach compared with the intensive treatment of the fully nonlinear system, denoted the nonlinear approach and (ii) over what range of MRI operating conditions can the linearised approach be expected to provide acceptable results for outputs of interest in an industrial context for MRI scanner design? The nonlinear approach is trusted as there is no approximation from a physical standpoint.
The aim of this paper is to address these two questions. This is achieved through the following novelties.
1. We revisit the linearisation of the nonlinear equations presented in our other work 13 and provide an alternative formulation, suitable for transient simulation of the full nonlinear problem, using a Newton-Raphson and alpha-type time integration schemes. We refer to this as the nonlinear approach. 2. We demonstrate numerically that the linearised approach is in excellent agreement with the nonlinear approach and derive estimates of the energy present in the nonlinear approach, which is not captured by the linearised approach. 3. We undertake a rigourous comparison of the linearised and nonlinear approaches in terms of outputs of interest for a set of challenging industrially motivated examples.
In this paper, we begin by presenting the coupled transient transmission problem and its associated initial conditions in Section 2. Then, in Section 3, we present the continuous weak treatment of the nonlinear and linearised approaches, derive estimates of the energy not captured by the linearised approach and present a transient Newton-Raphson procedure for the nonlinear iteration of the continuous fields in the two approaches. In Section 4, we briefly recall the spatial discretisation of the system by means of high-order finite elements. We include the temporal discretisation, where we have implemented a second order generalised -scheme, to allow for inclusion of numerical dissipation into the model. The system is then solved through an iterative monolithic Newton-Raphson solution procedure, which is summarised in an algorithm. We then conclude with numerical results in Section 5, which have been obtained by using the in-house software written by the group for this work.

COUPLED SYSTEM
Previously, in our other work, 13 we presented the fully coupled system of nonlinear equations describing the magnetic, mechanic, and acoustic behaviours of an MRI scanner. The work included the transmission conditions present at the interface between the conducting and nonconducting regions, as well as initial and far-field conditions of the system. From this coupled transmission problem, a simpler linearised scheme was introduced by a suitable additive split of the exciting current source J s (t) as J s (t) = J DC +J AC (t). In this section, we briefly recall the fully coupled transmission problem derived in our earlier paper, describe in detail the physical motivation for the additive split of J s (t), and provide a set of realistic initial conditions.

Transient nonlinear system
The governing coupled transmission problem, illustrated in Figure 2, describing the acousto-magneto-mechanical behaviour of MRI scanners is given by the following.
subject to an appropriate set of initial conditions. In (1), the magnetic flux density B is related to the magnetic field H and the magnetic vector potential A by B = H = ∇ × A, where is the magnetic permeability of the medium, is the acoustic pressure field in free-space region (of air) R 3 ∖Ω c , and u is the displacement field of the conducting bodies Ω c . The eddy currents J e (A) = A t and Lorentz currents J l (A, u) = u t × (∇ × A) are expressed in terms of the field variables, respectively. The notation describing the system of equations in (1) is the same as in our other work, 13 where the magnetic vector potential A, mechanical displacements u, and acoustic pressurêhave been introduced above. Here, we use u t ∶= u∕ t and u tt ∶= 2 u∕ t 2 and note that J s are the magnetic source currents, x is the position vector, is the electrical conductivity, and m (u) ∶= tr ( (u)) I + 2G (u), is the Cauchy stress tensor. In the above, and G denote the Lamé parameters, the material density, (u) ∶= (∇u + ∇u T )∕2 the linear strain tensor, I the identity tensor, T the transpose, n the unit normal vector, and the magnetic component of the Maxwell stress tensor.
The current source J s (t) in (1) can be decomposed through the additive split J s (t) = J DC + J AC (t), where J DC is a static current source provided by the main magnet and J AC (t) is a time-varying current source through a series of conducting (gradient) coils. Thus, the supports of J DC and J AC (t) are different and their (maximum) magnitudes and orientations are defined during the MRI design process. In all current MRI designs, they are also chosen so that |J DC | ≫ |J AC | by several orders of magnitude to provide a strong background static magnetic field, which is perturbed by weaker time-varying magnetic fields generated from the gradient coils and their interaction with the conducting components. Furthermore, J DC is always present and J AC (t) is only nonzero during imaging sequences. 43 In absence of J AC (t), the system in (1) simplifies considerably as there is no time variation and it is instead described by static solutions A DC , u DC ,̂D C . 13 If = 0 = 4 × 10 −7 H/m in R 3 , then u DC = 0 and̂D C = 0; however, in general, this is not the case.
Taking advantage of the fact that MRI scanners are typically maintained at full field strength, whilst in clinical use, as they are maintained at superconducting temperatures, and that the gradient (time-varying) fields are then only applied during imaging sequences, the initial conditions for (1) are set as the solution of the static field components

WEAK TREATMENT OF THE CONTINUOUS PROBLEM
The transient problem described by (1) represents a nonlinear system of coupled equations including electromagnetic, mechanic and acoustic effects in an MRI environment. We first establish the continuous weak form of this coupled nonlinear approach and present the continuous weak form of the linearised approach, previously derived in our other work. 13 Based on estimates of the energy present in the nonlinear approach, which is not captured by the linearised approach, presented in the Appendix, we summarise conditions under which we expect the linearised approach to provide an accurate representation of the solution of the nonlinear approach. A Newton-Raphson iteration is then presented to resolve the nonlinearity in the continuous weak forms.

Weak forms of the nonlinear and linearised approaches 3.1.1 Nonlinear approach
In our nonlinear approach, we assume that there exists a fixed point weak continuous solution for all weighting functions (A , u ,̂) ∈ (X × Y (0) × Z), denoted by a superscript , wherẽ and H(curl, R 3 ) and H 1 (R 3 ) have their usual definitions, eg, the work of Monk. 44

Linearised approach
The continuous version of the linearised approach, previously derived in our other work, 13 corresponds to seeking weak solutions (A DC + A AC (t)), (u DC + u AC (t)), (̂D C +̂A C (t)), which are obtained from a direct current DC and an alternating current AC stage.

DC-stage
The static solutions

AC-stage
The transient weak solutions (A AC , u AC ,̂A C )(t) ∈ (X(0) ×Ỹ (0, 0) ×Z(0))(0, T] are the solution of the linearised problem In the above, the notation DR A (A ; A DC , u DC )[A AC ] denotes the directional derivative of R A about the static solutions (A DC , u DC ) in the direction A AC with similar meanings for the other terms. 45 By deriving the required terms, this can be explicitly written as: Find (A AC , u AC ,̂A C )(t) ∈ (X(0) ×Ỹ (0, 0) ×Z(0))(0, T] such that for all (A , u ,̂) ∈ (X × Y (0) × Z). In the above, the linearised electromagnetic stress tensor, introduced previously in our other work, 13 is

Energy not captured by the linearised approach
In the Appendix, we estimate the energy present in the nonlinear approach, which is not captured by the linearised approach. We find that this will be small, at a continuous level, provided the following criteria are met for all time t: Provided that these conditions hold, we conjecture that the linearised approach will be in good agreement with the nonlinear approach. We argue that these conditions will be met on physical grounds for MRI scanners as follows: 1) requires that the velocities are small in comparison to and ||B DC || 2 L ∞ (Ω c ) , which is typically the case for MRI scanners; 2) requires that the displacements and strains (displacement gradients) are small in comparison to ||B AC || 2 L ∞ (Ω c ) and ||B AC || 2 , which is again typically the case for MRI scanners; 3) requires the pressure gradients in the coils are small, which is typically the case for MRI scanners under the Biot-Savart coil assumption; and 4) requires the B DC field obtained from the strong main magnet results in a ||B DC || L ∞ (Ω c ) that is orders of magnitude larger than ||B AC || L ∞ (Ω c ) , obtained from the weaker AC coils and smaller field perturbations caused by eddy and Lorentz currents in the conductors. In Section 5, we will demonstrate numerically that the 2 approaches are in excellent agreement for such problems. We include a simple model below to illustrate the applicability of these criteria.

Simple model relating field and current strengths
Whilst the true field strength B = H can only be found after solving (1) by means of the solution of the continuous weak forms for the nonlinear or linearised approaches, described in Sections 3.1.1 and 3.1.2, there is merit in also considering a simple model to relate B to the applied current density in the coils. This is because manufacturer's data, available for MRI scanners, is typically only quoted in terms of the maximum capable main and gradient field strengths at the central axis of the imaging bore. [6][7][8] Therefore, to determine the operating ranges of the current densities on in-use scanners, the relationships given by this simplified model offer useful insight.
In the simplified model, as shown in Figure 3, all transient, eddy current, and coupling effects are neglected, as is the mutual inductance between coils. Coil i has a constant cross-sectional area a i , is circular, and hence, rotationally symmetric and so is best expressed in terms of cylindrical coordinates (r, , z) and carries a uniform current density where B(0, , z) = B z (r = 0, z)e z along this axis, R i is the radius, and Z i the axial position of the ith current source relative to the scanner's central axes. Provided the aforementioned assumptions are enforced, this relationship can be applied to obtain the current densities that are associated only with the main magnet (J s = J DC ) or with the gradient coil (J s = J AC ) from their produced field strengths. It also provides a guide as to the ranges of ratios of |J AC |∕|J DC | and ||B AC || L ∞ (Ω) ∕||B DC || L ∞ (Ω) over which the linearised approach is applicable.

Nonlinear iteration
We now present a nonlinear iteration to resolve the nonlinear terms present in (4) and (5), resulting in linearised problems (not to be confused with linearised approach) being solved at each iteration. The unknowns in the linearised problems are continuous update fields, which we relate to the continuous solution of (4) and (5). We make explicit those terms, which are time dependent.

Newton-Raphson for the nonlinear approach
To resolve the nonlinearity in (4) a nonlinear iteration, which requires the solution of a linearised problem at each iteration, is performed. This is expressed in the form of the following iterative Newton-Raphson procedure: Find As standard, this linearisation is established under the assumption that the updates 45 Given , and performing iterations of this kind, leads to a sequence of continuous iterates ( [1] , u [1] ,̂ [ 1] ) (t), (A [2] , u [2] , with the goal of converging to the continuous solution is chosen to be sufficiently close to (A, u,̂)(t) this is expected to be the case. The requirements for the convergence of Newton-Raphson schemes are well documented and are consequently not discussed further, see e.g. 45 It is instructive to introduce the notation q(t) ∶= (q A , q u , q̂)(t) = (A, u ,̂)(t) and to rewrite (9) in the form of the Newton-Raphson iteration: [1] (t), q [2] (t), … ∈W(q DC , u D ) are expected to converge to q(t) ∈W(q DC , u D ). This formulation has the benefit that the forms of M, C, and K, are associated with the mass, damping, and stiffness contributions to the system, respectively, and can be found to be ‡ The system residual R q easily follows from (4) and is defined as

Newton-Raphson for the linearised approach DC-Stage
A nonlinear iteration, with a linearised problem being solved at each iteration, is also required to resolve the nonlinearity in (7). For this, we employ the Newton-Raphson procedure:  (5). This problem can also be formulated in terms of q DC

AC-stage
For completeness, we also recast (7) to: for all q ∈ W(0), which does not require iteration as it is linear in the unknown fields; hence, the reason we call it the linearised approach. The bilinear formsM,C andK are associated with the mass, damping and stiffness contributions in the linearised approach, and can be expressed in terms of the definitions in (13) as In this case the solution of (17) will contain only the transient alternating current (AC) component of the fields, and thus the complete description is q(t) = q DC + q (t). The linear nature of this system also allows for it to be transformed to the frequency domain by the application of a Fourier transform, which was the approach followed in our other work. 13

DISCRETISATION OF THE SYSTEM
In this section, we describe the spatial and temporal discretisation of (11) for the nonlinear approach and (16,17) for the linearised approach. The eventual aim of this formulation is to handle fully coupled nonlinear problems, which, for the nonlinear approach, will employ an iterative procedure at each time step to fully resolve the solution while, for the § The directional derivatives in (15) have been previously stated in our other work 13 and, in the following, we relate them to the explicit expressions in (13) to avoid repetition. linearised approach, will involve a static iterative solve followed by a linear time integration without iteration. The presentation includes a brief description of the spatial discretisation using hp-finite elements, which result in a semidiscrete system of coupled second order ODEs. We then present the temporal discretisation of the system of equations to allow for solutions in the time domain and close the section by presenting a solution algorithm for the nonlinear and linearised approaches, which summarises the steps required once the problems are fully discretised.

Spatial discretisation
The spatial discretisation of the linearised system of coupled acousto-magneto-mechanical equations in (16) and (17) is extensively covered in our other work 13 and the work of Bagwell 47 and the spatial discretisation of (11) is analogous. Hence, we only briefly summarise the approach. We limit ourselves to the treatment of idealised MRI scanners that are assumed to be rotationally symmetric in the azimuthal direction so that the problem is axisymmetric. Our domain then becomes the meridian plane Ω m = {(r, z) ∶ 0 < r < ∞, −∞ < z < ∞}, the current sources have the form J s = J s (r, z)e , and the unknown fields are A = A (r, z)e , u = u r (r, z)e r + u z (r, z)e z and̂=̂(r, z), where e r , e , and e z are the standard bases for the cylindrical coordinate system (r, z, ). However, by projection of Ω m the full 3-dimensional results are still achieved.
When transformed to the axisymmetric domain, the spaces in which the weak solutions are sought in the variational statements (11,16,17) must be adapted following the approach described in our other work 13 and the work of Bagwell. 47 To overcome the need for the solution in weighted spaces, to ensure the fields are well behaved at the radial axis, we transform the fields as A = rÂ , u r = rû r and note thatÂ ∈ H 1 (Ω m ) andũ ∶= [û r , u z ] ∈ (H 1 (Ω m )) 2 and recall from our other work 13 that the acoustic pressurê(r, z) ∈ H 1 (Ω m ) presents no difficulty. The continuous update fields and test functions in the Newton-Raphson iterations are also treated in an analogous manner.
To achieve a finite computational domain, Ω m is truncated by the introduction of infinite elements 48 for the far-field treatment ofÂ and a perfectly matched layer approach 49 for the far-field treatment of̂. The far-field decay of the update fields is also treated analogously. For the details, we refer again to our other work 13 and the work of Bagwell. 47 Then, a nonoverlapping partition (, , ) of Ω m is achieved through the generation of a hybrid unstructured triangular-quadrilateral mesh, where  denotes the set of vertices,  the set of edges, and  the set of cells. In each element, we employ the H 1 conforming hierarchic finite element basis functions proposed by Schöberl and Zaglmayr, 50,51 which allow for arbitrary increases in element order p and local refinement of the mesh spacing h, for the spatial discretisation of the fields A ,ũ, and̂and the associated continuous update fields.
For a discretisation consisting of uniform order p elements, and a mesh consisting of N  vertices, N  edges, N  T triangular cells and N  Q quadrilateral cells, there will be N = N  + ( − 1)N  + 1 2 ( − 2)( − 1)N  T + ( − 1) 2 N  Q degrees of freedom for (each component of) each field.

Nonlinear approach
Following a spatial discretisation using a uniform order approximation, the semidiscrete vector of updates and unknowns for problem (11) for each time t and the corresponding semidiscrete solution to the nonlinear approach reduces. Finding q [k] (t) at each time t ∈ (0, T] satisfying the following initial value problem: q | to 0 will be quadratic. 45 In the above, M, C [k] , K [k] , and R [k] q are the discrete linearised mass, damping and stiffness matrices, and residual vector, respectively, which are obtained by discretising the terms in (13) and depend, with the exception of the mass matrix, on the solution at iteration [k]. Explicit expressions for the matrices can be found in the work of Bagwell. 47 The vector q DC ∈ R 4N is the converged solution of (21) below.

Linearised approach DC-stage
In a similar manner, the linearised approach takes the form of a time-independent Newton-Raphson iteration. Find with the iteration performed until |R DC[k] q | < TOL, which follows from the discretisation of (16). The convergence of |R DC[k] q | to 0 will again be quadratic.

AC-stage
The linearised approach also requires the solution of the semidiscrete initial value problem. Find q which follows from the discretisation of (17), whereM,C,K, andR q are the discrete mass, damping and stiffness matrices, and the residual vector, respectively, obtained by discretising the terms in (18). Explicit expressions for the matrices can be found in the work of Bagwell. 47 The semidiscrete solution for the linearised approach is q(t) = q DC + q (t).

Temporal discretisation
There are a considerable number of alternative time integration techniques that could be employed for the temporal integration of (20) and (22), such as Euler schemes, 52 trapezoidal schemes, 53 and the Newmark method. 54 In this paper, we choose to adopt a second-order generalised-scheme (see the work of Chung and Hulbert 55 ) to discretise our system of equations. This scheme is designed to allow for tailor-made numerical dissipation in the solution and is of sufficient degree, given the second-order nature of (20) and (22). Given the complexity of the coupled system, it is difficult to know the true initial conditions for the full transient problem, when applying a forced excitation. The numerical dissipation of the generalised-scheme is therefore beneficial as it can be used to damp out any artificial frequencies that pollute the solution. There exists also a number of high-order integration schemes, such as backward differentiation formulas, 56 explicit Runge-Kutta and leap-frog type methods, 57,58 that allow for higher-order accuracies in time and larger time steps, although with increased temporal accuracy comes a significant increase in computational cost. However, for our problem, the required time step size is dictated by the need to capture the excitation frequency and resonant frequencies of the system (typically in the region of 5000 Hz). For this reason, the need for higher-order accuracy is significantly outweighed by the requirement of smaller time step sizes anyway, and thus, the generalised-scheme appears a sensible choice.

Generalised-time integration scheme (second order)
In this section, we focus on the development of the temporal discretisation of (20), which can similarly be developed for (22). The implicit form of this scheme evaluates the system of nonlinear equations in (20) at an intermediate time step where C| [ The generalised-scheme, described in the aforementioned work, 55 n+1 , known as an acceleration-based formulation. This can be manipulated into a displacement based formulation, whereq [k] n+1 and q [k] n+1 are expressed in terms of q [k] n+1 , as where , Δt is the time step size and ∞ denotes the user-specified value of the spectral radius in the high-frequency limit. ¶

Predictor-multi-corrector step
A predictor-corrector approach has been followed from the implementation standpoint. 59 The prediction step [k = 0], based on a initial guess of the displacement field q [0] n+1 , is defined as where (26a) and (26b) are consistent with (25a) and (25b), respectively. In other words both the velocity and accelerations predictors preserve second-order accuracy, as discussed in the work of Jansen et al. 60 If we return our attention to (25b) and (25a) and substitute the fields at iteration step [k+1], we obtain that the update variables of the first-and second-order fields arëq

Fully discrete nonlinear approach
Using the relations between the updates, in (27), we may rewrite the discrete system, in (23), in terms of the update in the zeroth-order solution vector q n+1 at time t n+1 as the following Newton-Raphson procedure. Find q with the predictors (26) used to provide the initial guessesq [0] n+1 andq [0] n+1 , which are used to start the iteration over [k] until |R [k] q | < TOL. ¶ Traditionally with the generalised method, is used to refer to the first stability parameter; however, given our formulation reserves for the material conductivity, we use instead the alternative symbol .

Fully discrete linearised approach AC-stage
The fully discrete version of (23) for the linearised approach is similarly expressed. Find q n+1 ∈ R 4N for n = 0, 1, 2, … , N Δt = T∕Δt such that where the system matricesM,C, andK are independent of the previous temporal solutions, andR q depends only on the evaluation of J AC at time level t n+1− and so can be solved in a single step, which is a simplified version of (28).

Solution algorithm
A general algorithm for computing the transient variation in the fields for both the nonlinear and linearised approaches, proposed above, under a generalised-scheme is summarised in Algorithm 1. In the nonlinear approach, the algorithm involves 2 nested loops, the first over n performs the integration of the coupled system over time and, the second, over [k] resolves the nonlinearity through the application of a Newton-Raphson iteration. In the linearised approach, this simplifies to a single loop over n to perform the time integration. Furthermore, Figure 4 summarises the steps required to construct the discretised Newton-Rapshon iterative scheme in (28) from the transient nonlinear transmission problem (1).

NUMERICAL EXAMPLES
In this section, we discuss two benchmark MRI test problems to justify that the conditions presented in Section 3.2 do hold in practice. We study the linearity of the fields for a range of different current density values in the coils to determine the validity of the linearised approach.

Test magnet problem
First, we consider solutions to the industrially relevant test magnet problem, which was previously presented in our other work 13 as the "Siemens benchmark problem." The conducting region Ω c of the test magnet consists of three metallic shields, known as the outer vacuum chamber (OVC) Ω OVC c , 77 • K radiation shield Ω 77K c , and 4 • K helium vessel Ω 4K c , each with different material parameters ( , , , G, ). The exact geometries and material parameters of the conducting components are commercially sensitive and as such are not displayed in this paper. By choosing to study only the Z-gradient coils and noting that the currents in the main coils and geometry of the conducting components are rotationally symmetric, the problem may be reduced to an axisymmetric description and solved in cylindrical coordinates (r, , z) on the meridian plane Ω m . The full 3-dimensional representation of this simplified MRI scanner is depicted in Figure 5.
A pair of main magnet coils, each with a static current source J DC = J DC (r, z)e , are located on the outside of the three shields and a pair of Z-gradient magnet coils, each with alternating current source J AC (t) = J AC (t, r, z)e , are located within the imaging bore, both of which are assumed as Biot-Savart conductors and are located in free space. Realistic excitations of the gradient coil are nonsinusoidal in nature; however, for the purposes of comparison between the linearised and nonlinear approaches, we restrict ourselves to sinusoidal excitations, where the current density is described by J AC (t) = Re(|J AC |e i t ). We consider several frequencies of excitation of = 2 [1000, 1500, 2000]rad/s, which lie outside of the resonance region of ≥ 2 [3500], predicted in our other works. 13,34 The magnitude of the static |J DC | and gradient current sources |J AC | to be considered are obtained from manufacturer's data. [6][7][8] This data quotes the maximum capable flux density on the central axis of the imaging bore (r = 0) and to obtain the corresponding current densities we have used the model described in Section 3.2.1, the results of which are summarised, for key clinical field strengths, in Table 1.
Given criteria 4) in Section 3.2, and hence the implication of the current density ratios, the greatest nonlinearity should appear for a magnet with the weakest static field of 1.5 T and the strongest gradient field of 100 × 10 −3 T/m. This would result in a ratio between the static and gradient current density values of |J AC |∕|J DC | ≈ 7.2% for this problem. We therefore study the two approaches across a range of current density ratios of |J AC |∕|J DC | = [5,10,15,20]% to provide a rigourous test of the linearised approach for applications of higher levels of nonlinearity than current MRI scanners are capable of. The problem is subject to the following boundary conditions: we fix the Dirichlet boundaries of the conductors to u D = 0 to fix the conductors in space, as in Section 2, and we set the value of the magnetic vector potential on the outer boundary to A = 0 due to the eddy current decay. The initial conditions of the problem are defined by (3). We treat this problem computationally for both the linearised and nonlinear approaches. We truncate the nonconducting free-space region, comprised of air, and create the domain Ω m , which is the same as in our other work. 13 In terms of spatial discretisation, we analyse the solution for a single unstructured mesh of 2842 triangles of maximum size h = 0.25 m, but with substantial refinement in the conductors Ω m c , resulting in 570 elements in Ω m c . The mesh parameters used are the same as in our other work 13 ; however, due to improvements in the mesh generator used, 61 resulted in a better distribution of elements and hence fewer elements. We apply a single layer of 18 infinite (quadrilateral) elements on the outer boundary of Ω m to resolve the static decay of the magnetic field, so that the boundary condition A = 0 is effectively imposed at infinity. We consider p-refinement for elements of order p = [1,2,3,4,5]. The temporal discretisation used to resolve these waves are studied for a time step size of Δt = 2 ∕( N Δt ), where we vary the number of points per wavelength of the excitation frequency N Δt = [10,20,30,40]. The spectral radius of the -scheme time integrator ∞ allows for damping of certain frequency regimes. For ∞ = 1 the amplitude of the wave is fully preserved and no damping of any high frequencies is introduced. For ∞ = 0, the scheme is fully dissipative and higher frequency waves are completely damped; however, the damping of physical modes also occurs. Given that we wish to recover only physical modes in our problem, we therefore chose the value of ∞ = 0.8, to allow for numerical damping of any nonphysical high frequencies induced through the forcing, whilst still preserving the physical lower frequency waves.
Solutions to this problem are obtained by applying Algorithm 1 for both the nonlinear and linearised approaches across the range of current densities and excitation frequencies discussed.

Convergence study in outputs of interest
We consider quantities of industrial interest for both the mechanical and electromagnetic fields. For the electromagnetic field, we are concerned with the output power dissipation in the conducting components P o which reduce, in a time harmonic representation, to where A( ) and u( ) are the complex amplitudes of their respective fields as A(t) = Re (A( ) e i t ) and u(t) = Re (u( ) e i t ), as described in our other work. 13 To compute these quantities for the test magnet problem, using either the transient solutions obtained by the linearised or nonlinear approaches, Algorithm 1 is run at specific N Δt , Δt, ∞ for a sufficiently long time until the field responses reach steady state periodic solution. Equation (30) is approximated by Gaussian quadrature for the spatial integration (with sufficient quadrature points corresponding to the degree of the integrand) and a Trapezoidal rule for the temporal integration so that no further approximations are made. In order to determine suitable spatial resolution in the solution, we perform a p-refinement study. We apply Algorithm 1 for N Δt = 30, ∞ = 0.8 using the linearised approach to compute the output power and kinetic energy across the frequency  Figure 6. For p ≤ 3 the computed curves, illustrated in Figure 6, do not match one another and so the results have not yet reached convergence. The curves for p = 4 and p = 5, however, are practically indistinguishable for both the output power and kinetic energy and suggests that for elements of order p = 4 the results are sufficiently converged. We hence decide to adopt elements of p = 4 for all subsequent computations of this problem on the mesh specified previously.
A similar study to determine the required time step size was also carried out, where a range of number of time steps per wavelength N Δt = [10,20,30,40] for p = 4 elements were studied. The results suggest that N Δt = 30 offers sufficient temporal resolution to capture the amplitudes of the dominant frequency as well as frequencies twice the dominant frequency 2 . The importance of this frequency doubling will be explained later in Section 5.  47 Using these parameters, in Section 5.1.4, we compute the quantities in (30) by applying Algorithm 1 for the linearised and nonlinear approaches in the time domain and compare them with our previous frequency domain solver. 13 However, first we consider the results for the linearised and nonlinear approaches for transient electromagnetic and mechanical fields.

Electromagnetic field
To perform comparisons between the linearised and nonlinear approaches, we compare the transient response of the magnetic field for both approaches. Given that the output power of the conductors P o Ω c , described above, is driven by the temporal derivative of the magnetic vector potential A t , we measure and compare the response of this field. Figure 7 summarises the transient results of both approaches for the = 2 [2000]rad/s sinusoidal excitation. The graphs in the left-hand column plot the time signal obtained from both the linearised approach (in red) and nonlinear approach (in black). The right hand column plots the corresponding frequency spectrum of the signal, obtained by performing fast Fourier transforms (FFTs). We can see from this figure that the linearised approach provides an accurate approximation of the magnetic vector potential across the full range of current density ratios. The relative norm of the difference in the time signals is 1.829 × 10 −4 for |J AC |∕|J DC | = 5% and 2.1888 × 10 −4 for |J AC |∕|J DC | = 20%. The magnitude of the frequencies across the spectrum are almost identical and the linearised approach is even capable of capturing all the fundamental frequencies, around = 2 [3500 − 4000]rad/s. Figure 8 plots the transient response and corresponding frequency spectrum for a current density ratio of |J AC |∕|J DC | = 20% and a range of excitation frequencies. These plots again illustrate the agreement between the linearised and nonlinear approaches for different excitation frequencies.

Mechanical field
We now compare the response of the mechanical field for both the linearised and nonlinear approaches. Given that the kinetic energy of the conductors E k Ω c , described in Section 5.1.1, is driven by their mechanical velocity u t , we measure and compare the response of this field. Figure 9 summarises the transient results of both approaches for the = 2 [1000]rad/s sinusoidal excitation. The graphs in the left-hand column plot the time signal obtained from both the linearised approach (in red) and nonlinear approach (in black). From the plots it appears as though the two approaches offer good agreement, especially when looking at the time signals. The relative norm of the differences between signals being 4.039 × 10 −3 for |J AC |∕|J DC | = 5% and 1.465 × 10 −2 for |J AC |∕|J DC | = 20%. However, in the frequency spectrum there appears to be an extra frequency at = 2 [2000]rad/s that is picked up in the nonlinear approach, but not in the linearised approach. This term appears due to the nonlinearity in the Maxwell stress tensor (2), which we can observe is quadratic in the magnetic field. From a decomposition of the magnetic field into static and dynamic components, as shown in Section 3.1.2, it can be shown that this nonlinear term comprises of a product of the dynamic component of the field with itself, which disappears in the linearised approach. This term causes a frequency doubling effect, in other words it results in a component of excitation of the mechanical field that is double the frequency of the AC currents. So for a = 2 [1000]rad/s wave, as in Figure 10, an excitation at = 2 [2000]rad/s would also appear, which matches exactly with the results obtained. However, from Figure 10, it is clear that the magnitude of this term is far smaller than that of the amplitude associated with the exciting frequency and thus has little effect on the solution. As the ratio of the current densities increases so too does the magnitude of this term. However, even for a ratio of 20% the magnitude is still several orders smaller than the main excitation and smaller also than the most dominant resonant frequencies. Thus the magnitude of this doubled frequency component provides a useful measure in determining the nonlinearity of the problem. Figure 10 plots the transient response and corresponding frequency spectrum for a current density ratio of |J AC |∕|J DC | = 20% and a range of excitation frequencies. We see from these plots that for the AC current frequencies of 1500 Hz, the agreement between the two approaches results in almost indistinguishable time signals. The relative norm of the differences between signals being 0.0349 for = 2 [2000]rad/s AC currents' differences in the time signal become more visible. This is because the doubled frequency excitation component of 4000 Hz lies within the resonance region. When exciting close to the resonance region, the problem results in matrices of high condition numbers, which are close to singular. Consequently, the tangent stiffness matrix inversion becomes more challenging and less reliable and, as a result, can lead to differences in the amplitudes across the frequency spectrum. Despite this effect however, the differences in the time signal are still very small and the characteristics of the system and prediction of the resonance region remain well captured by the linearised approach.

Comparison of linearised and nonlinear approaches
To benchmark the accuracy of the solution from the linearised approach with the nonlinear approach, we compare the computation of the outputs of interest for the two approaches, presented in Section 5.1.1, across the range of frequencies ≤ 2 [5000]rad/s. Figure 11 illustrates the outputs of interest computed by the linearised approach in both frequency and time domain and the nonlinear approach in time domain, using the definitions in (30) and (31) for a current density ratio of |J AC |∕|J DC | = 10%.
Using the frequency domain approach, described in our other work, 13 47 From the plots, the curves produced by the linearised and nonlinear approaches are almost indistinguishable across the frequency spectrum which suggests that the linearised approach provides a very accurate approximation to the full nonlinear approach across the full spectrum for |J AC |∕|J DC | = 10% for the two outputs of interest. In fact, given that the individual fields, analysed in Sections 5.1.2 and 5.1.3, also show very good agreement for |J AC |∕|J DC | = 20%, we can hypothesise this to be the case also for the outputs of interest as they are directly related.  We now compare the displacements of the mechanical shields by plotting the displacements of the OVC, in three dimensions, at interesting instances of the time signal for the radial velocity (u r ) t computed for the case of |J AC |∕|J DC | = 20% and = 2 (2000)rad/s in Figure 12. The chosen time instances across the time signal are plotted, where the difference in the mechanical velocities between the linearised and nonlinear approaches is noticeable. The contour plots of the displacement fields are all scaled such that the colour maps between the linearised and nonlinear approaches are the same. The displacements in Figure 12 are scaled by several orders of magnitude to show visually the displacement shapes of the shield. We can see that for all snapshots the differences between the two approaches are almost indistinguishable. This suggests that even for J AC |∕|J DC | = 20% and an excitation frequency of = 2 (2000)rad/s, where the doubled frequency component of = 2 (4000)rad/s resides in the resonance region, see Section 5.1.3, that the linearised approach provides accurate and comparable results to the full nonlinear approach. For the results of the other shields, see our other work. 47 The average computational timings, per time step, for the two approaches across a range of different element orders p are summarised in Table 2. The comparison between the computational timings of the two approaches suggests that the linearised approach is orders of magnitude more efficient in terms of computational cost than the nonlinear approach. For low-order p = 1 elements the linearised approach requires 383 times less computational effort than the nonlinear approach. Whereas, for higher-order p = 5 elements the linearised approach requires 164 times less computational effort than the nonlinear approach. This speed-up factor appears to offer an inverse exponential behaviour with p, which is due to the higher requirement on the solver for higher-order elements. Nevertheless, the computational timings displayed in the table suggest that the linearised approach offers orders of magnitude increase in computational efficiency over the nonlinear approach.

Realistic magnet problem
We now consider a more realistic problem that represents, very accurately, the sorts of MRI scanner designs currently used in clinical operation. The geometry is illustrated in Figure 13A. This problem consists of a similar construction to the previous problem, where the conducting region is comprised of the three radiation shields each with different material parameters ( , , , G, ). The geometry of the radiation shields, however, is more complex and their topology represents that of closed cylindrical shells of trapezoidal cross section, with curved face end sections. However, despite the increased complexity in topology, the geometry is still cylindrical and can be treated as axisymmetric. Again, the exact geometries and material parameters of the conducting components are commercially sensitive and as such are not displayed in this paper. The configuration of the static main coil consists of the same block cross section as the test magnet problem, but contains more sets of coils including also a set of secondary coils, which act to minimise the magnetic stray field by reversing the polarity of the magnetic field. The gradient coils of this problem represent a far more realistic Z-gradient coil structure, which also contains a set of primary and secondary coils for shielding. The coils are sourced in the same way as the test magnet problem. The cross section of this problem, projected onto the positive half axisymmetric meridian domain Ω m (r, z ≥ 0), is illustrated in Figure 13B. Given criteria 4) in Section 3.2, and hence the implication of the current density ratios, the greatest nonlinearity is expected to appear for a magnet with weakest static field of 1.5 T and strongest gradient field of 100 × 10 −3 T/m. This would result in a ratio between the static and gradient current density values of |J AC |∕|J DC | ≈ 12.3% for this problem. We therefore study the two approaches for an extreme case of a current density ratio of |J AC |∕|J DC | = 15% to rigourously test the linearised approach. The form of the excitation of J AC is as described in Section 5.1.
In terms of boundary conditions, we truncate the nonconducting free-space region, comprised of air, and create the domain Ω m and apply those similar to the test magnet problem, whereby we set the Dirichlet boundaries of the conductors to u D = 0 to fix the conductors in space. We also set the magnetic vector potential to be A = 0 on the outer boundary.
In terms of spatial discretisation, we analyse the solution for an unstructured mesh of 19 218 triangles of maximum size h = 0.25 m, but with substantial refinement in the conductors Ω m c , resulting in 4085 elements in Ω m c . We then apply a single layer of 40 (quadrilateral) infinite elements on the outer boundary of Ω m to resolve the static decay of the magnetic field. Having carried out similar convergence studies to those presented for the test magnet problem, we have chosen to use order p = 4 elements, N Δt = 30 time steps per excitation frequency , and set ∞ = 0.8 for the solutions to this problem.
We now compare the displacements of the mechanical shields by plotting the velocity magnitude of the system across the time signal using both approaches for |J AC |∕|J DC | = 15% and = 2 (1500)rad/s in Figure 14. Given that the differences in both the time signal and the displaced shapes of the OVC between the two approaches are almost indistinguishable, we highlight instead snapshots of the mechanical displacement in the inner OVC shell at various time instances across the time signal for the linearised approach. For the full results of the other shields (see our other work 47 ). We plot the displaced OVC and static magnetic field for the static problem, in Figure 15, and the corresponding gradient fields, at various time instances, in Figure 16. The displacements in Figures 14, 15, and 16 are scaled by several orders of magnitude to show visually the displacement shapes of the OVC.

CONCLUSION
In this paper, we have described a framework for both the linearised and nonlinear approaches to obtain transient solutions to the acousto-magneto-mechanical coupling in MRI scanners. We have demonstrated numerically that the linearised approach is in excellent agreement with the nonlinear approach and have estimated the energy present in the nonlinear approach, which is not captured by the linearised approach. In terms of the two questions we posed at the start of the paper, through the numerical examples presented in Section 5, we have shown that: (i) there is accurate agreement between the two approaches and that, not only does the linearised approach provide accurate approximations to the transient response of the fields, but also accurately predicts the quantities of interest and the resonance behaviour of MRI scanners, and (ii) the linearised approach provides accurate results across the full region of interest in the frequency spectrum ≤ 2 [5000]rad/s for a range of current density ratios J AC ∕J DC ≤ 20% for a range of outputs of interest. In terms of current clinical MRI scanner applications, this ratio is restricted to around 4% to 12% and thus our analysis validates the use of the linearised approach in providing accurate solutions in current and future MRI scanner design, with orders of magnitude saving in the computational cost over the nonlinear approach.

R u term
Following similar steps to above we find that where C depends on the size of Ω c and Ω N c .