A linearised hp–finite element framework for acousto‐magneto‐mechanical coupling in axisymmetric MRI scanners

We propose a new computational framework for the treatment of acousto‐magneto‐mechanical coupling that arises in low‐frequency electro‐magneto‐mechanical systems such as magnetic resonance imaging scanners. Our transient Newton–Raphson strategy involves the solution of a monolithic system obtained from the linearisation of the coupled system of equations. Moreover, this framework, in the case of excitation from static and harmonic current sources, allows us to propose a simple linearised system and rigorously motivate a single‐step strategy for understanding the response of systems under different frequencies of excitation. Motivated by the need to solve industrial problems rapidly, we restrict ourselves to solving problems consisting of axisymmetric geometries and current sources. Our treatment also discusses in detail the computational requirements for the solution of these coupled problems on unbounded domains and the accurate discretisation of the fields using hp–finite elements. We include a set of academic and industrially relevant examples to benchmark and illustrate our approach. Copyright © 2017 The Authors. International Journal for Numerical Methods in Engineering Published by John Wiley & Sons, Ltd.


INTRODUCTION
Magnetic resonance imaging (MRI) has become a widely used and popular tool in the medical industry capable of diagnosis of many medical ailments, such as tumours, damaged cartilage and internal bleeding as well as its use in neuroimaging. The most popular type of magnet used in these devices are superconducting magnets, consisting of conducting wire contained within a supercooled vessel of liquid helium known as the cryostat, which achieve the field strengths required for high-resolution imaging. Figure 1 shows a typical setup of an MRI scanner, consisting essentially of four main components: the main magnet coils, secondary magnetic coils, the cryostat (comprised of a set of radiation shields that encapsulate the liquid helium immersed magnetic coils) and resistive gradient coils. A set of main magnetic coils produce a strong uniform stationary magnetic field across the radial section of the scanner, required to align the protons of the hydrogen atoms in the patient in the axial direction. The secondary magnetic coils are used to minimise the large stray magnetic fields arising outside of the scanner. The cryostat consists of a set of metallic vessels in vacuum to maintain the supercooled magnet temperatures and shield the coils from radiation. A set of resistive coils inside the imaging volume, known as gradient coils, produce pulsed magnetic field gradients that excite certain regions of the protons to generate an image of the patient.
The simulation of MRI scanners also builds on the expanding literature devoted to magnetomechanics and coupled problems including FEM-BEM coupling [23,24], magneto-mechanical damping machines [25], magneto-mechanical effects on material parameters [26], enhanced basis functions for magneto-mechanical coupling [27] and strongly coupled systems [28]. When the coupling includes acoustic effects, careful treatment of the far field boundary is also required and techniques such as perfectly matched layers (PML), FEM-BEM coupling, infinite elements and absorbing boundary conditions have been developed for this (see, e.g. [29] for a recent review).
Although commercial codes provide an efficient environment for many problems, our interest lies in providing a low-cost dedicated industrial design tool. This tool should operate over a wide range of frequencies more effectively in order to predict the response from an analysis of general gradient coil time signatures, which can be decomposed by Fourier expansion. As such, not only do we require a reliable means of treating the coupled nature of acousto-magneto-mechanical problems, but we also need to be able to accurately resolve the potentially small skin depth in conducting components as well as accurately resolve the propagation of acoustic waves. The aforementioned commercial codes are typically designed with low-order FE discretisations in mind, which require dense meshes for handling the small skin depths and wave propagation at higher frequencies. On the other hand, hp -FE discretisations offer possibilities for high accuracy on locally refined meshes and have been shown to accurately resolve the skin depth in the conducting components [30], handle the complex coupling [16][17][18] and to resolve the propagation of the acoustic waves [31] and we, therefore, choose to adopt them here. Furthermore, in the computation of unbounded domains using a PML, they have also been shown to offer superior performance [32].
Building on the established hp -FEM methodology, the main novelty of our work is to provide a new rigorous theoretical framework for the simulation of acousto-magneto-mechanical effects in MRI scanners, which forms the basis of our design tool. We undertake the consistent linearisation of the transient equations and arrive at a simplified monolithic single-step strategy in the case of harmonic gradient coil excitations. It greatly improves on our previous work [2], which required nonphysically motivated simplifying assumptions and resulted in a fixed point strategy with a growth of iterations for increasing frequency. We also extend our framework to incorporate the effects of the acoustic field and propose a rigorous set of interface conditions, which couple the various physics together. The entire framework is suitable for three-dimensional geometries discretised by hp -FEM, but such simulations would be prohibitively expensive for the industrial design cycle. As our interest lies in the development of a rapid design tool, we focus on the simulation of problems on axisymmetric geometries.
The presentation of the material proceeds as follows: In Section 2, we outline the governing equations, the coupling between the fields and present a fully coupled non-linear transient transmission problem. Section 3 presents the consistent linearisation of the weak form of the transmission problem, and we derive a simplified monolithic strategy in the case of harmonic gradient coil excitations. Then, in Section 4, we briefly discuss the reduction to axisymmetric geometry, the computational far field treatment and the hp-FE discretisation. Section 5 presents a selection of numerical case studies to validate each physical field independently and highlights the computational challenges of small skin depth and high-frequency wave propagation. Then we present a series of results for complex coupled case studies of industrial relevance before closing with concluding remarks.

COUPLING APPROACH
In the following, we describe the governing equations and coupling methodologies that link the electromagnetic, mechanical and acoustic behaviour of an MRI scanner. We begin, in Section 2.1, with the transient eddy current model to describe the electromagnetic response from a conducting region Ω c when illuminated by a low-frequency background magnetic field. This field arises from a current source J s with support in an unbounded region of free space R 3 ∖Ω c , as illustrated in Figure 2. Then, in Section 2.2, we present the corresponding transient mechanical response of Ω c resulting from electromagnetic stresses generated in this region. Finally, in Section 2.3, we present the transient  acoustic response resulting from the vibration of Ω c in R 3 ∖Ω c . The complete coupled transmission problem is stated in 2.4.

Electromagnetic description
For our chosen application, the transient eddy current approximation of Maxwell's equations is valid [2], where the displacement current terms are neglected because of the high conductivity of the conducting components and the low frequency of the exciting currents. A rigorous justification involves the topology of the conducting region [33]. Defining E, H, D, B as the electric, magnetic, electric flux and magnetic flux intensity vectors, respectively, and introducing a vector potential A such that B = ∇ × A, this model can be described by ‡ where we assume x to be measured from the centre of Ω c . The previous equation assumes the regions to be homogenous and isotropic, such that B = H and J o = E where denotes the electrical conductivity and the magnetic permeability. The solenoidal external current sources J s are assumed to lie in free space, R 3 ∖Ω c , where = 0 and = 0 = 4 × 10 −7 H/m. The term J l = u∕ t × ∇ × A denotes the Lorentz currents where u is the mechanical displacement in the conducting region Ω c . The vector potential A satisfies the transmission conditions on the conductor boundary Ω c where [·] Ω c denotes the jump on this interface and n is a unit outward normal vector to Ω c .

Mechanical description
The conducting region Ω c is assumed to behave elastically and the mechanical displacements u to satisfy ‡ The temporal gauge has been applied in Ω c and the Coulomb gauge in R 3 ∖Ω c [30]. Note we set A = O ( |x| −1 ) as x → ∞ according to the mathematical model described by Ammari, Buffa and Nédélec; [33] here, the big O notation implies that the rate is at least as fast as |x| −1 and can be faster in practice; for details, see the aforementioned paper.

Coupled transmission problem
Combining the statements from the previous Sections 2.1, 2.2 and 2.3, we arrive at the following transmission problem for describing our coupled acousto-magneto-mechanical system in a time where we have chosen to set the initial conditions for the fields to be zero, corresponding to a system at rest at t = 0. An illustration of the fully coupled system is shown in Figure 4. For our desired application, the system (7) is excited through the current source J s (t). In practice, the application allows for the decomposition J s (t) = J DC +J AC (t), where J DC corresponds to the static current source of the main magnetic coils and J AC (t) the transient current source of the gradient coils [2]. This decomposition, illustrated in Figure 5, allows us to introduce the following static problem: where we have assumed a similar decomposition of the Dirichlet displacement condition u D = u D DC + u D AC (t).

LINEARISATION
With developing a FE framework for the approximate solution of (7) and (8) in mind, we linearise weighted residual statements of the transmission problems. To this end, it is convenient to introduce the following which will be used to describe the weak solutions to the dynamic and static transmission problems, where H(curl, R 3 ) and H 1 (R 3 ) have their usual definitions (e.g. [34]). We start with the treatment of the simpler static problem (8) and then continue to our approach for the transient system (7).

Linearisation of the static problem
Consider possible weak solutions (A DC , u DC ,p DC ) ∈ X × Y(u D DC ) × Z and the associated residuals for all (A , u ,p ) ∈ X × Y( ) × Z. The directional derivatives of these residuals are where n = n + = n − and we have introduced the linearised electromagnetic stress tensor At a continuous level, equations (9,10) can be used to introduce the Newton-Raphson iteration: Find The initial guesses are such that To permit the computational solution of (11), a spatial FE discretisation is required, which we will discuss in Section 4. However, at this stage, it is already useful to note that, because of the specific nature of the equations, (11a) can be solved independently, followed by (11c) and then (11b) without iteration. Moreover, if the system is solved monolithically, the solution will converge to ( In our previous approach [2], we neglected the effects of the static displacement, driven by the static magnetic field, as we were primarily interested in computing the output power of the system that depends only upon the transient displacements, which were assumed to be harmonic. However, in this new formulation, we also include the static effects of the displacements to allow for a rigorous treatment of the linearised transient scheme and maintain consistency of the physical fields. for all (A , u ,p ) ∈ (X × Y( ) × Z). The directional derivatives of these residuals are

Consider possible transient weak solutions (A, u,p)(t) ∈ (X × Y(u D AC ) × Z)[0, T] and the associated residuals
In the preceding equations, we have found it convenient to use the alternative form of ∇ · e introduced in (5) when linearising Rp.
One strategy for solving the temporal system, after spatial discretisation, would be to adopt a discrete time integration scheme and then apply the Newton-Raphson algorithm at each time-step to solve the non-linear equations using the directional derivatives computed previously. However, such an approach is computationally expensive and, in the interests of developing a fast computational technique, we adopt a different strategy.
Rather than integrating the equations in time and solving the non-linear equations at each time step, we choose instead to linearise the full time dependant equations about the static solution. This linearisation in the context of MRI scanners is motivated by the knowledge that the static DC current source is several orders of magnitude stronger than the weaker AC time varying source, leading to a strong DC field and a weaker time varying field. Similar techniques, involving the additive split of a non-linear problem to a series of linear problems, have been successfully applied to the field of computational mechanics such as analysis of structural membranes [35][36][37], high-order mesh generation [38] and in biomedical applications [39]. In this case, the residuals of the dynamic problem becomẽ and the associated directional derivatives take the form Finally, recalling that A DC , u DC ,p DC are all time invariant, we see that the residuals and the directional derivatives in (14) and (15), respectively, are linear in the time dependent terms A , u ,p and J AC . This motivates the time harmonic representations where denotes the angular frequency of the driving current in the gradient coils in the case of a harmonic excitation. In reality, the gradient coils are driven using non-harmonic excitations, but their time signals can be decomposed into its different frequency modes and the same approach applied. The solution of the linear harmonic problem becomes as follows: find We may then decompose the full temporal solution into its static and time varying components, which in the case of a single frequency excitation are given by The directional derivatives in Equation (16) explicitly become In the absence of pressure, the formulation previously proposed in [2] for the dynamic case is identical to that obtained in the preceding equations; however, previously, a number of assumptions had to be made in order to drop terms to arrive at the considered system, which is now no longer the case. Moreover, in the preceding formulation, the equations are rigorously established through a linearisation of the dynamic system and the solutions can be obtained by a single monolithic solve rather than an iterative fixed point scheme, previously proposed.

COMPUTATIONAL TREATMENT
In this section, we discuss the key components for the computational solution to Equations (16) and (11). This includes the reduction for rotationally symmetric geometries in Section 4.1, the far field treatment in Section 4.2 and the hp -FE discretisation in Section 4.3.

Reduction for rotationally symmetric geometries
The formulation proposed in Section 3 is valid for general three-dimensional domains involving a conducting region surrounded by an unbounded region of free space containing the current sources.
To a first approximation, the geometry of an MRI scanner is close to cylindrical and the strong DC current source J DC has only the component J DC in cylindrical coordinates (r, , z). However, of the three sets of AC gradient coils, it is only the z-gradient coil that exhibits rotational symmetry, and in our quest for a rapid industrial design tool, we must neglect the x and y gradient coils. Under these assumptions, our simplified MRI scanner is rotationally symmetric with respect to the azimuth, the problem reduces to solving for A = A (r, z)e , u = u r (r, z)e r + u z (r, z)e z ,p =p(r, z) [2]. Here e r , e and e z denote the standard basis of the cylindrical coordinate system. The reduction of the full threedimensional problem to the axisymmetric meridian (r, z) plane Ω m is shown in Figure 6. However, by projection of this plane, full three-dimensional results are still achieved.
When transformed to the axisymmetric domain, the spaces in which the weak solutions are sought in the variational statements (16) and (11) must also be adapted. In general, this leads to necessity to seek for solutions in weighted spaces to ensure the fields are well behaved at the radial axis [40]. To avoid the complexity of weighted spaces, we transform the fields A = rÂ and u r = r̂r according to that proposed in [2], and note that Â ∈ H 1 (Ω m ),ũ ∶= [̂r, u z ] ∈ H 1 (Ω m ). The acoustic pressure presents no difficulty becausep(r, z) ∈ H 1 (Ω m ). The treatment of the bilinear and linear forms associated with the terms in (11) and (16) follow similar steps to that presented in [2], and hence, we only give as an example where Ω m c is the projection of Ω c on to the meridian plane, ∇ m ∶= [ ∕ r, ∕ z] T and the factor of 2 results from the integration of the azimuthal direction. This factor similarly appears in the treatment of all of the other terms in the linear system and therefore cancels.

Far field treatment
The strong forms of the dynamic and static problems include radiation and decay conditions, which describe the behaviour of A, A DC ,p andp DC as |x| → ∞. To allow the computational treatment of the problem, the unbounded free space region R 3 ∖Ω c is truncated at a finite distance from Ω c and the region Ω n is created, which contains all the current sources. The three-dimensional computational domain Ω ∶= Ω n ∪ Ω c ⊂ R 3 becomes Ω m ∶= Ω m n ∪ Ω m c ⊂ R 2 for axisymmetric problems. This means that integrals over R 3 ∖Ω c in (11) and (16) become integrals over Ω n and then, in turn, Ω m n . On Ω m , the static decay of A, A DC ,p DC is approximated by fixing Ω m to be located sufficiently far from the region of interest and setting Naturally, the quality of the approximation improves as the size of Ω n (Ω m n ) is increased. However, the harmonic pressure field, and in particular the treatment of DR AĈ p (p )[p], cannot merely be approximated by truncating the domain and fixingp = 0 at Ω m , because this would result in reflections that would pollute the computational domain because of its wave behaviour. The radiation condition (4) in the continuous problem describes the correct decay of this field, which must also be approximated computationally. Therefore, we add a PML [32] Ω pml (Ω m pml ) to the exterior of Ω n (Ω m n ) so that the computational domain now becomes Ω = Ω c ∪ Ω n ∪ Ω pml (Ω m = Ω m c ∪ Ω m n ∪ Ω m pml ). For terms other than DR AĈ p (p ) [p] in (16), Ω pml can be merely thought of as a free space extension of Ω n . However, the aforementioned term (18) is treated differently as where 1 , Λ 2 are both complex functions of position in the layer and reduce to identity on Ω m n ∩ Ω m pml . The coefficients of these functions can be established through a complex coordinate stretching of the domain Ω m n following the approach in [32]. For the axisymmetric case, this complex stretching is equivalent to introducing the complex position-dependent functions In the preceding equation, the complex coordinate transform z s is described as a power law in terms of the distance to the layer d s and thickness of the layer t s , as illustrated in Figure 7b. The prime indicates differentiation with respect to the argument, and explicitly, we choose where s = [r, z]. The choice of a power law of degree 5 and user-defined thickness t s is somewhat arbitrary, provided that the resulting complex field behaviour in Ω m pml is properly resolved. If this is accomplished, the acoustic pressure field is absorbed without reflection and we can setp = 0 on Ω m .

hp-FE discretisation
Building on the previous success in [2], we choose to adopt the hp-FE H 1 conforming basis functions, proposed by Zaglmayr and Schöberl [41,42], for accurately discretising the fields Â ,ũ to account for the domain truncation introduced in Section 4.2. The corresponding discrete linear harmonic problem for ( A hp e ,ũ hp ,p hp ) = (r Â hp e , r̂r hp e r + u z hp e z ,p hp ) ∈X( ) ∩ rX hp e × Y(u D AC ) ∩ (rX hp e r + X hp e z ) ×Z( ) ∩ X hp can be developed analogously. We choose not to give the explicit formulae for the terms because they can be found by following the approach for axisymmetric problems, which we have previously described in [2]. The complete algorithm is summarised in Algorithm 1.

NUMERICAL RESULTS
In this section, we present a series of academic and industrial numerical examples in order to demonstrate the capabilities of the presented framework. Firstly, we include examples with analytical solutions to demonstrate the independent validation of the electromagnetic, mechanical and acoustic fields for uncoupled problems in Sections 5.1. Secondly, we include validation of coupled physics problems and application to challenging industrial benchmarks in Section 5.2.

Validation of single-physics problems
5.1.1. Conducting sphere in uniform alternating magnetic field. A closely related problem to the solution of (1), with J s = J l = 0, is that of a conducting object located in free space excited by a uniform harmonic background magnetic field of amplitude H 0 and frequency . In this case, −1 0 ∇× A → H 0 as |x| → ∞. For a spherical conductor Ω c = {x ∶ |x| 2 ⩽ R 2 } of radius R, permeability s and conductivity s , as illustrated in Figure 8, an analytical axisymmetric solution is presented in [43] for A .
We A hp | Ω n = A | Ω n and we set A hp ∶= A hp . We generate a coarse mesh of 578 unstructured quasiuniform triangular elements of maximum size h = 0.5m and use here, and subsequently, a blending function approach to represent the exact geometry of the sphere's surface [44]. This function avoids any geometrical error in the solution because of coarse approximation of the boundary. In light of the smooth nature of the solution, we consider uniform polynomial enrichment corresponding to p = 1, 2 … , 10 on this mesh and plot in Figure 9 the relative error measures ||A − A hp || L 2 ∕||A || L 2 , ||A − A hp || H 1 ∕||A || H 1 against the number of degrees of freedom for varying frequencies of the alternating magnetic field. Figure 9a shows the convergence of the error against the number of degrees of freedom on a logarithmic scale, where each point represents a polynomial refinement and the different curves correspond to different frequencies and error measures. Each line indicates a downward sloping   curve suggesting that the convergence is exponential. This is confirmed by plotting the error on a logarithmic scale against the number of degrees of freedom raised to the power 1∕2 on an algebraic scale in Figure 9b. After a pre-asymptotic region, each curve becomes a straight line indicating that convergence of the numerical to the analytical solution is exponential with respect to the square root number of degrees of freedom for p-refinement of this problem. This corresponds to the expected rate for smooth solutions, as reported by Babuška and Guo [45] ‖ . As the frequency of the alternating magnetic field increases, the gradient of the lines in Figure 9b reduces, indicating that, although still exponential, the rate of convergence is lower. Physically, this is due to the smaller skin depths s = √ 2∕( ), [43] which characterises the depth to which the eddy currents J o = E ≈ i A hp e decay to 1∕e of their surface value, associated with higher . It is possible to improve the gradient in the plots in Figure 9b by using a graded mesh towards the sphere's surface. To illustrate the different fields and skin depth effects for different frequencies, Figure 10 shows the contours of |B| ≈ |∇×(A hp e )| for the various frequencies of the converged solutions. This figure illustrates the smaller s for higher and the need to use higher fidelity discretisations to capture the solution with the same degree of accuracy.

Mechanical shell.
To verify the treatment of the elasticity system, we consider the solution of equation (2) for e = 0 and static displacements u ≠ u(t) for the case of a spherical mechanical shell Ω = Ω c = {x ∶ r 2 i ⩽ |x| 2 ⩽ r 2 o } of inner and outer radii r i , r o , respectively. The inner and outer surfaces of the shell are subject to traction conditions −p i n and −p o n resulting from internal and external pressuresp i ,p o respectively, on different parts of Ω N as illustrated in Figure 11. This problem is axisymmetric and has the analytical solution [46], which can be expressed in terms of the cylindrical displacement components (u r , u z ).
Specifically, we solve the problem corresponding to r i = 0.5m, r o = 1m, E = 210 × 10 9 Pa, = 0.49,p i =p o = 10 4 Pa so that the shell is nearly incompressible. As described in Section 2.2, we must fix part of the boundary of the shell to avoid it floating away. We solve a suitably simplified Figure 11. Mechanical shell subject to interior and exterior pressure: problem setup. version of (20) for a single iteration where we choose to fix the displacements DC[0] uhp = u r e r +u z e z on a small boundary segment Ω D according to the analytical solution. The region of computation corresponds to Ω m = Ω m c = {(r, z) ∶ r 2 i ⩽ (r 2 + z 2 ) ⩽ r 2 o }, which we discretise by a quasi-uniform mesh of 68 unstructured triangular elements of maximum size h = 0.5m. We perform the same p-refinement study that was described in Section 5.1.1 and now measure convergence using ||u − u hp || L 2 ∕||u|| L 2 , ||u − u hp || SNS ∕||u|| SNS where u hp ∶=ũ hp . The results shown in Figure 12a and 12b illustrate similar trends to those shown previously indicating that exponential convergence with respect to the number of degrees of freedom raised to the power 1∕2 is also achieved through p-refinement for the mechanical problem. In particular, p-refinement serves as a method for overcoming volumetric locking that is known to be associated with the displacement formulation of elasticity for nearly incompressible material [47,48] and leads to exponential rates of convergence of the error measured in the SNS norm. Although we do see stagnation of convergence when the norms of the error reach 10 −13 , which coincides with the numerical precision of our computation. The displacements in the radial and axial directions of the shell, obtained using p = 10, are illustrated in Figure 13.  [4,4]m andp hp = 0 on Ω m . We fix a quasi-uniform mesh of 408 unstructured triangular elements of maximum size h = 0.5m and perform the same refinement study of p = 1, 2 … , 10 as previously and measure convergence using ||p sc −p schp || L 2 ∕||p sc || L 2 , ||p sc −p schp || H 1 ∕||p sc || H 1 , wherep schp ∶=p hp . Figure 15a and 15b illustrates similar downward sloping trends, as shown in the previous examples. However, in this case the pre-asymptotic region is now affected by the increase in wave number. For higher wave numbers there is an initial stage of increase in error, which results from wave dispersion effects, discussed in [31] and [50]. This effect is overcome by further increasing p and eventually   results in the same expected exponential rates of convergence as before, confirmed in Figure 15b. This again indicates that exponential convergence with respect to the number of degrees of freedom raised to the power 1∕2 is also achieved through p-refinement for the acoustic problem, provided sufficiently high refinement is used to eliminate numerical dispersion. For the case of k = 4 ∕3 m -1 for p > 7, the convergence behaviour is suboptimal because of the effect of the PML, which is an approximate absorbing boundary condition, but nevertheless, accurate solutions are still obtained. The finest solution, using p = 10, of the scattered pressure field arising from the incident pressure field for wave numbers k = [4 ∕3, 10, 30] m -1 are illustrated in Figure 16.

Coupled multi-physics problems
5.2.1. Acoustic wave scattering of thin elastic shell. We consider a coupled acousto-mechanical problem consisting of a thin elastic shell of thickness t, mid surface radius R and material parameters s , and E. The shell is placed in a background medium described by f and c f and is illuminated by harmonic incident pressure fieldp in . The configuration is illustrated in Figure 17. This problem requires the solution of Equation (7) in absence of electromagnetic coupling and naturally lends itself to a harmonic treatment. For thin shells, the solution to this problem can be approximated by  the Kirchhoff shell theory [49,51] and the total pressure exterior becomesp =p in +p sc +p r . For the incident fieldp in =p 0 e ikz , the componentp sc corresponds to the hard scattering by the sphere and p r to the radiated pressure [49,51].
We  ])∖{(r, z) ∶ r 2 + z 2 ⩽ (R − t∕2) 2 }m 2 with the same PML settings as in Section 5.1.3. As in Section 5.1.2, to avoid the shell from floating away, we fix a small boundary segment Ω D to have displacements u = 0. The problem is driven by the incident pressure field in the form of a Neumann condition set on the external boundary of the shell {(r, z) ∶ r 2 + z 2 = (R + t∕2) 2 } as n · ∇p hp = −n · ∇p in and the coupling according to the interface conditions in equation (7). The inside boundary of the shell is left free and the acoustic effects inside of the shell are ignored.
Given that the convergence rates of the mechanical and acoustic fields have already been verified in Sections 5.1.2 and 5.1.3, we now instead directly compare the effectiveness of both h-refinement and 1344 S. BAGWELL ET AL.
p-refinement with the analytical solution for k = 4 ∕3m −1 . In Figure 18, we plot various hp-enriched solutions for the line segment 1.025 m ⩽ r ⩽ 5.6 m, z = 0, taken from the outer surface of the shell to the truncated boundary. For both h-refinement and p-refinement, the computed solution tends to the analytical for r ⩽ 4 m. However, for h-refinement, a mesh of h ≈ 0.1m with 15 752 elements, with 7779 unknowns, is required to obtain good agreement with the analytical solution, with a level of accuracy of O(10 −2 ). On the other hand, using p = 2 on a mesh with 530 elements requires only 1021 unknowns for comparable accuracy in the solution. If we further refine p to p = 4, then the number of unknowns increases to 4163, but with improvement in the relative accuracy of two orders of magnitude. The PML is defined by the grey area, in which the computed solution is non-physical and absorbed. Figure 19 shows the comparisons in computed and analytical solutions for higher wave numbers of k = [10, 30]m −1 . For both cases, the order of p = 4, used to obtain the finest solution in Figure 18, offers reasonable agreement with the analytical solution and is able to capture the higher frequencies of the waves. However, for increasing wave numbers, our computed solution requires even further prefinement in order to accurately capture the solution in these regions because of the aforementioned dispersion effects. The solution case for p = 4 and a suitably refined solution of each case is plotted Figure 18. Elastic shell subject to an incident acoustic pressure field: Effects of h-refinement and p-refinement on the acoustic pressure field. Figure 19. Elastic shell subject to an incident acoustic pressure field: effects of p-refinement for high wave number k on the acoustic pressure field.
against the analytical solution in Figure 19a and b. The mesh density required, for p = 1 elements, to capture the high-frequency wave effects at k = 30m −1 becomes impractical for such geometries. The interaction between the pressure field and the displacements of the shell is illustrated in Figure 20, where the computed deformed shape of the shell is plotted in the surrounding acoustic field.
The success of high-order hp discretisations for the preceding case studies motivates our strategy for the following industrial example, which does not have an analytical solution.

Siemens benchmark problem.
We now consider an industrially relevant benchmark problem, proposed by Siemens Magnet Technology, in which a simplified quarter-size representation of an MRI scanner was modelled and previously presented in [2]. The problem setup comprises of the same main components illustrated in Figure 1, with a reduced complexity in the coil configuration. The setup comprises 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 , which make up Ω c and each with different material parameters ( , , , E, ). A pair of main coils, with static current source J DC , are located on the outside of the three shields and a pair of gradient coils, with alternating current source J AC (t), are located within the imaging bore. Both are assumed as Biot-Savart coils and are illustrated in Figure 21.
We treat this problem computationally for two cases; in the first case, we apply a suitably simplified version of Algorithm 1 in which we neglect the acoustic effects and focus on the purely magneto-mechanical coupling mechanisms, as in [2]. In the second case, we consider the fully coupled acousto-magneto-mechanical systems in Algorithm 1. We truncate the non-conducting Figure 20. Elastic shell subject to an incident acoustic pressure field: deformed shell interacting with surrounding acoustic pressure field.  . Similar mode shapes are also obtained for the other shields [54]. For each spike in the kinetic energy, or resonant frequency, a corresponding three-dimensional axisymmetric mode shape of the shield is included. Higher resonant frequencies excite higher order sinusoidal modes as the figures illustrate. The resonant frequencies of the kinetic energy in the shields coincide with the frequencies experienced in the output power, which suggests that the primary source of excitation in the conductors is that of the eddy currents dissipated in the shields.
Finally, we illustrate the complex behaviour of the magnetic flux lines, acoustic contour lines and deformed structure in Figure 26 for a range of frequencies in both the eddy current dominant (low frequency) and the resonance (high frequency) regions. This figure illustrates that in the low-frequency region (120Hz ⩽ f ⩽ 1500Hz), a patient in the bore can become exposed to sufficiently higher noise levels than that of the exterior region of the scanner, with the sound radiating and decaying outwards in space. This is due to the dominance of the harmonic magnetic field arising from the gradient coils, located inside the imaging volume, which gives rise to the source term in the acoustic Helmholtz system in (3). As the frequency increases (1500Hz ⩽ f ⩽ 3700Hz), the effect of the mechanical resonance begins to dominate and the acoustic field is further excited by the displacement of the shields and the sound intensity outside of the scanner increases. The case of f = 4075Hz illustrates a higher frequency mode shape and the effect of the displacement on the acoustic field. In this case, the magnetic field is further perturbed because of the increase in Lorentz currents resulting from the acoustic excitation. Notably, in each case, the greatest sound intensity is that inside the bore tube, suggesting that the highest noise levels are experienced by the patient.

CONCLUSION
In this paper, we have provided a new rigorous theoretical framework for the simulation of acoustomagneto-mechanical effects in MRI scanners, which provide the basis of our design tool. We have provided a consistent linearisation of the transient equations and have arrived at a simplified monolithic single-step strategy in the case of harmonic gradient coil excitations. This greatly improves our previous work, [2] which required non-physically motivated simplifying assumptions and resulted in a fixed point strategy with a growth of iterations for increasing frequency. We have further extended our approach to include acoustic effects and discretised the resulting framework by hp-FEM to ensure a robust tool that provides the accurate resolution of small skin depths and wave propagation effects. The optimal convergence of the physical fields is demonstrated for a series of challenging single and multi-physics test cases in axisymmetric coordinates. The predictive capability has been illustrated by applying our approach to rapid evaluation of the performance of a MRI scanner model whose acousto-magneto-mechanical response has been analysed in detail. The next steps of our research involves developing a full three-dimensional simulation tool for MRI scanners based on hp-FEM.