An Analytical Approach to the Design of Quiet
Cylindrical Asymmetric Gradient Coils in MRI
LAWRENCE K. FORBES,1 MICHAEL A. BRIDESON,1 STUART CROZIER,2 PETER T. WHILE1
1
School of Mathematics and Physics, University of Tasmania, Hobart 7001, Tasmania
School of Information Technology and Electrical Engineering, University of Queensland, St. Lucia,
QLD 4067, Australia
2
ABSTRACT: We consider the design of asymmetric gradient coils in a conventional
cylindrical bore magnetic resonance imaging system. The gradient coils are switched on
and off during the scan, and are therefore subject to Lorentz forces that cause them to
buckle during the switching process. This in turn creates a pressure wave within the coil,
and that gives rise to acoustic noise. We present a simplified linearized model for the
deflection of the coil due to electromagnetic forces, which is amenable to solution using
analytical methods. Closed-form solutions for the coil deflection and the pressure pulse
and noise level within the coil are obtained. These are used to design new coil winding
patterns so as to reduce the acoustic noise. Sample results are shown both for
unshielded and shielded gradient coils. Extensions of this model are indicated, although
it is suggested that the advantages of the present closed-form solutions might then not
be available, and fully numerical solutions may be needed instead.
Ó 2007 Wiley Periodicals, Inc.
Concepts Magn Reson Part B (Magn Reson Engineering) 31B: 218–236, 2007
KEY WORDS: gradient coils; switching functions; Lorentz forces; elastic deformation;
acoustic noise reduction
INTRODUCTION
The design of electromagnetic coils to produce a
desired magnetic field is an important problem in
magnetic resonance imaging (MRI) technology, and
has been the subject of a great deal of investigation.
Conventional MRI equipment has cylindrical geometry, and the permanent magnet, gradient coils, and radio-frequency probes are arranged around the surface
of the cylinder. The patient then lies within the cylinReceived 18 April 2007; revised 8 June 2007;
accepted 15 June 2007
Correspondence to: Lawrence K. Forbes; E-mail: larry.forbes@utas.
edu.au
Concepts in Magnetic Resonance Part B (Magnetic Resonance
Engineering), Vol. 31B(4) 218–236 (2007)
Published online in Wiley InterScience (www.interscience.wiley.
com). DOI 10.1002/cmr.b.20095
Ó 2007 Wiley Periodicals, Inc.
drical bore of the device. For imaging purposes, it is
typically required to generate a magnetic field with
(axial) z-component that varies linearly in the three
orthogonal x-, y- and z-directions, over a certain
region of interest within the coil. This is accomplished by the gradient coils, which are switched on
and off repeatedly during the scan in many imaging
applications. A description of MRI technology is
given in the book by Jin (1).
The design problem with gradient coils is therefore to arrange windings on the surface of the cylinder in such a way as to produce the desired linear
gradient field with the region of interest, namely the
DSV (diameter of the sensitive volume). It is now
known that, in general, this is an ill-conditioned
design problem, in the sense that two very different
winding patterns on the cylinder can nevertheless
give similar magnetic fields within the DSV; consequently, the inverse problem, in which the desired
magnetic field is specified in advance, does not have
218
DESIGN OF ASYMMETRIC GRADIENT COILS IN MRI
a clear unique solution for the required winding
pattern.
One method of overcoming this difficulty is the
target-field approach suggested by Turner (2, 3). This
method considers the MRI coil to be notionally infinite in length. The advantage of this approach is that
Fourier transforms are then immediately available to
solve for the required current density on the cylindrical coil, since the Fourier transform and its inverse are
unique, and this avoids ill-conditioning. Of course,
coils in practice are finite in length, but Turner’s technique can be modified to account for this, to some
extent, by requiring the current density to vanish
at points outside the actual coil (4). This is usually
combined with the use of a damping function in the
Fourier-transformed space, so as to guarantee the
convergence of the transforms, and Jin (1) refers to
this as ‘‘apodization.’’
Alternative solution methods are also available for
this design problem, and likewise overcome the illconditioned nature of the task. Crozier and Doddrell
(5) used a simulated annealing optimization algorithm to find the location of discrete coils of wire and
the required current in each. This technique is
extremely robust, although it takes much iteration to
converge. Those authors applied it to the design of
zonal coils, for which the geometry of the windings
is relatively simple, but the method is expected to be
more difficult to use in situations that demand more
elaborate winding patterns around the cylinder.
Forbes and Crozier (6–8) developed an inverse
method for designing winding patterns on cylindrical
coils, in which the finite length was accounted for explicitly in the formulation. Their method is a type of
target-field approach, similar in some respects to that
of Turner (2), and it treats the Biot-Savart law as a
first-kind integral equation for the current density on
the coil, once the desired magnetic field has been
specified. It makes use of Tikhonov regularization to
treat the ill-conditioned nature of the problem [see,
for example, Delves and Mohamed (9) or Neittaanmäki et al. (10, p 62)] and can cope with designs for
zonal or tesseral fields of a quite general nature, with
coils that are either shielded or unshielded. Recently,
a similar technique has been applied to the design of
shielded biplanar coils (11).
For many patients, an MRI scan can be a claustrophobic experience (12). Discomfort is no doubt exacerbated by the fact that the interior of the cylindrical
coil is a noisy environment, caused by the fact that
the gradient coils are subject to significant Lorentz
forces [see Jackson (13)] that result from the interaction of the magnetic field in the coil with the currents
on its surface. As the gradient current is pulsed in
219
time, the coil deforms accordingly, giving an acoustic pressure wave that is experienced as noise by the
patient.
A method for reducing the noise was suggested by
Chapman and Mansfield (14). They relied on the concept of balancing the Lorentz forces on every conductor in the coil, by matching each individual wire
with a segment carrying an equal but opposite current. Their model suggests that, for a sinusoidal current at 500 Hz, a reduction of 10 dB in noise output
could be expected. A similar reduction occurs at
higher frequencies also, but is regarded by Mansfield
et al. (15) as less than satisfactory, since noise levels
can be about 130 dB, which represents a dangerous
situation. These authors have suggested active acoustic control using force-balanced gradient coil layers,
in which extra coils are added with the sole purpose
of opposing the deflections caused by the gradient
set, although these coils will themselves affect the
desired magnetic field.
Mechefske et al. (16) have used a finite-element
method to model the deformation of the gradient coil
in response to Lorentz forces, and have shown that
their predictions are in good agreement with experimental measurements. At frequencies of about 2 kHz
they observe noise levels in the 120–140 dB range.
More recently, Shao and Mechefske (17, 18) have
used analytical methods in the analysis of noise
within cylindrical ducts, assuming a monochromatic
sinusoidal behavior with time. Their methods are
similar in some respects to those presented here, and
as they suggest (18), analytical methods are computationally faster and allow greater understanding of the
solution. In the present article, we show that analytical methods have the added advantage that they may
also be incorporated directly into optimization techniques in a reasonably straightforward way, and this is
a great benefit in design. Nevertheless, this requires a
physical model of sufficient simplicity to permit the
use of these techniques, and further discussion of this
point is given in concluding Section 5.
THE MATHEMATICAL MODEL
We consider a cylindrical duct of length 2L, with the
z-axis of a Cartesian coordinate system lying along
its axis. The coil is thus located over the interval L
< z < L as shown in Fig. 1. The inner radius of the
coil is a, and it is supposed that the primary windings
are located at this radius. The coil is made of material
of thickness h, Young’s modulus E, and Poisson’s
ratio v. Its outer radius is denoted b ¼ a þ h, as indi-
Concepts in Magnetic Resonance Part B (Magnetic Resonance Engineering) DOI 10.1002/cmr.b
220
FORBES ET AL.
field is specified on the surface of this cylinder, and
also on a second cylinder with radius c2 (<c1) within
the region of interest, for greater accuracy. These
two inner target cylinders are shown in Fig. 1. When
a shield winding is present on the outer surface r ¼ b
of the coil, then a third target cylinder of radius c3
(>b) is also employed, and the intention is to make
the magnetic field drop to zero on that surface. Given
the cylindrical geometry of the coil and target zones
shown in Fig. 1, it will be convenient to use cylindrical polar coordinates (r, y, z) according to the usual
relations x ¼ rcos y and y ¼ rsin y.
Magnetic Field Calculation
Figure 1 A sketch of the cylindrical coil and target
zones. The primary is located at inner radius r ¼ a and
the shield is at r ¼ b ¼ a þ h. The coil mid-radius is rM
¼ a þ h/2. The inner target cylinders of radii c1 and c2
are located asymmetrically over the interval pL < z < qL.
The field is minimized on outer radius c3. [Color figure
can be viewed in the online issue, which is available at
www.interscience.wiley.com.]
From the Biot-Savart law, the magnetic induction
field B at a field point r in the space surrounding the
coil is found from
m
BðrÞ ¼ 0
2p
ZL Z2p
L
0
ZL Z2p
L
In this expression, the functions jPy and jSy refer to
the (azimuthal) y-component of the current density
vectors on the primary and shield windings, respectively. Under steady-state conditions, they are related
to the (axial) z-components jPZ and jSZ of the current
density vectors by the continuity equations divjP ¼ 0
m0
2p
krb 0 rk3
dA0
½1
In this expression, the source points ra0 and rb0 lie on
the inner and outer surfaces r ¼ a and r ¼ b of the
coil, respectively. The two vectors jP and jS (A/m)
are, respectively, the primary and shield current densities on conductors placed on the inner and outer
surfaces of the coil. (It is assumed in the derivation
of Eq. [1] that current flows on both sides of the conducting elements).
The induction field B in Eq. [1] is related to the
magnetic field H (A/m) through the constitutive relation B ¼ m0 H. We are primarily concerned with the
axial component of this field, which may be obtained
from Eq. [1] in the form
½a2 þ r 2 2ar cosðy0 yÞ þ ðz0 zÞ2 3=2
b
2p
dA0
kra 0 rk3
ZZ
ðr rb 0 Þ jS ðrb 0 Þ
r¼b
½r cosðy0 yÞ a jPy ðy0 ; z0 Þ
ðr ra 0 Þ jP ðra 0 Þ
r¼a
cated in the Figure, and the shield windings are
located on this outer surface.
The region of interest within the coil is taken to
be a cylinder of radius c1 positioned co-axially with
the coil. In this investigation, we allow this region to
be located either symmetrically or asymmetrically
with respect to the coil length, to account for situations in which the patient is positioned toward one
end of the coil to overcome claustrophobia, for
example. This is achieved by locating the region of
interest over some interval pL < z < qL in which the
dimensionless constants p and q have the relationship
1 < p < q < 1. When p ¼ q the region of interest
is positioned symmetrically; otherwise an asymmetric gradient coil will be designed. The desired target
a
HZ ðr; y; zÞ ¼
2p
ZZ
0
dy0 dz0
½r cosðy0 yÞ b jSy ðy0 ; z0 Þ
½b2
þ
r2
0
2br cosðy yÞ þ
ðz0
2 3=2
zÞ
dy0 dz0 : ½2
on r ¼ a and divjS ¼ 0 on r ¼ b. It follows from a
vector identity that the two components of each current density vector can be obtained from streamfunctions cP(y,z) and cS(y,z), on the primary and
shield locations, respectively, by means of the
expressions
Concepts in Magnetic Resonance Part B (Magnetic Resonance Engineering) DOI 10.1002/cmr.b
DESIGN OF ASYMMETRIC GRADIENT COILS IN MRI
jPy ¼
qcP
;
qz
jPZ ¼
1 qcP
a qy
on r ¼ a
[3a]
from which it follows that the azimuthal component
must be
and
jSy ðy; zÞ ¼
jSy
qc
¼ S;
qz
jSZ
1 qcS
¼
b qy
[3b]
þ
Equally spaced contours of each streamfunction
correspond to actual winding patterns on each coil,
as shown in the tutorial article by Brideson et al.
(19).
As the axial current density components jPZ and
jSZ are required to be zero at each end of the coil z ¼
6L, it follows that this component on the primary
surface r ¼ a may be taken to have the form
jPZ ðy; zÞ ¼
N
M X
X
2mL
m¼1 n¼1
npa
þ
N
M X
X
½PPmn cos my
m¼1 n¼1
þ
QPmn
and
cP ðy; zÞ ¼
N
X
2L
n¼1
þ
np
8
9
npðz þ LÞ>
:
;
sin my cos>
2L
PP0n
np
m¼1 n¼1
½4b
½PPmn cos my
8
9
npðz þ LÞ>
:
;
þ QPmn sin my sin>
2L
½4c
respectively. Similar arguments based on Eq. [3b]
may be made for the shield windings on r ¼ b. The
axial current density component may be taken to be
jSZ ðy; zÞ ¼
M X
N
X
2mL
m¼1 n¼1
npb
QSmn
½PSmn sin my
8
9
npðz þ LÞ>
>
:
;;
cos my sin
2L
½5b
cS ðy; zÞ ¼
N
X
2L
8
9
npðz þ LÞ>
:
;
PS0n cos>
np
2L
n¼1
M X
N
X
2L
np
m¼1 n¼1
þ
QSmn
½PSmn cos my
8
9
npðz þ LÞ>
>
:
;: ½5c
sin my sin
2L
In these expressions [4] and [5], the integers M and
N are to be chosen suitably large, and the primary
coefficients PPmn , QPmn and shield coefficients PSmn ,
QSmn are as yet unknown. The algorithm for determining them will be outlined in Section 3.
Deformation of the Cylinder
8
9
npðz þ LÞ>
>
:
;
cos
2L
M X
N
X
2L
N
M X
X
½PSmn cos my
and the streamfunction for the shield windings must
take the form
½4a
8
9
npðz þ LÞ>
:
;
PP0n sin>
2L
n¼1
8
9
npðz þ LÞ>
>
:
;
sin
2L
8
9
npðz þ LÞ>
:
;
þ QSmn sin my cos>
2L
þ
From Eq. [3a], the azimuthal current density component and the streamfunction on the primary then
become
N
X
PS0n
m¼1 n¼1
½PPmn sin my
8
9
npðz þ LÞ>
P
>
:
;:
Qmn cos my sin
2L
jPy ðy; zÞ ¼
N
X
n¼1
on r ¼ b:
221
½5a
It is necessary now to consider the equation of motion
for the deformation of the coil in response to Lorentz
forces. In general this is a difficult task, as the governing Cauchy momentum equation is highly nonlinear,
even for Hookean materials. To make progress with
analytical methods, it is necessary to approximate
under the assumption that the deformations are small
relative to the radii of curvature of the coil itself, and
that the material is elastically isotropic with density
rC (kg/m3). For linearized elastodynamic problems,
the momentum equation then becomes (20, p 102)
rC
3
X
qsji
dvi
¼ rC fi þ
;
qxj
dt
j¼1
i ¼ 1; 2; 3:
[6]
This vector equation has been written in component
form, and the terms vi (m/s) are the components of
the velocity vector in the three orthogonal space
directions. If quantities ui (m) represent the components of the displacement vector, then vi ¼ dui/dt.
This time derivative, and that in Eq. [6], are material
Concepts in Magnetic Resonance Part B (Magnetic Resonance Engineering) DOI 10.1002/cmr.b
222
FORBES ET AL.
derivatives, but in the linearized analysis to be presented here, the inertial contributions to these terms
will be ignored.
The symbols sij in Eq. [6] represent the nine components of the stress tensor in the coil. These are
related to the components of the strain tensor eij
through the generalized Hookean law (21, p 262)
sij ¼ dij
3
X
ekk þ 2Geij ;
i; j ¼ 1; 2; 3
[7]
k¼1
in which the two constants L and G are Lamé coefficients that are related to the Young’s modulus E and
the Poisson ratio v for the coil material through the
expressions
G¼
E
vE
and ¼
:
2ð1 þ vÞ
ð1 þ vÞð1 2vÞ
i; j ¼ 1; 2; 3
[9]
and is thus symmetric, since it is invariant if indices i
and j are interchanged.
The body force vector f per mass, with components fi in Eq. [6], is obtained from the Lorentz force
law in the form
f¼
1
ðj BÞ:
rC h
q2 u 1
¼ ðj BÞ þ ð þ GÞrðr uÞ þ Gr2 u:
qt2
h
[11]
The two Lamé coefficients L and G in this relation are given in Eq. [8].
In view of the cylindrical geometry of this
problem, it is appropriate to express the mechani-
[12]
where uR is the radial component of displacement in
the outwardly directed radial direction given by the
unit vector er and similarly for the other two components in the two other orthogonal directions. The
three components of Eq. [11] may now be written in
cylindrical polar coordinates by means of the relations given in Batchelor (22, p 602), for example.
This gives
rC
rC
q2 uR 2
q
¼ BZ0 jy þ ð þ GÞ ðr uÞ
h
qr
qt2
8
9
u
2 quy >
:r2 uR R
;
þ G>
r 2 r 2 qy
q2 uy 2
1 q
ðr uÞ
¼ BR jZ þ ð þ GÞ
2
h
r qy
qt
8
9
u
2 quR >
: r 2 uy y þ
;
þ G>
r 2 r 2 qy
rC
[13a]
½13b
q2 uz
2
q
¼ BR jy þ ð þ GÞ ðr uÞ þ Gr2 uz ;
2
h
qz
qt
[13c]
in which the divergence in polar coordinates is
[10]
(Other body forces, including gravity, are ignored in
this analysis, since they are small in comparison to
the Lorentz force). Equations [6]–[9] are now combined with this Lorentz term (Eq. [10]) and quadratic terms in the displacement are ignored. This
results in the linearized deflection equation for the
coil, which may be expressed in vector form as
rC
u ¼ uR er þ uy ey þ uZ eZ ;
[8]
The symbol dij is the Kronecker delta term, and it
has the value 1 when i ¼ j and 0 otherwise. According to linearized elasticity theory, the strain tensor is
given in terms of gradients of the displacement components according to the formula
9
8
1>
qui quj >
>
eij ¼ : þ >
;;
2 qxj qxi
cal deformation Eq. [11] in cylindrical polar coordinates (r,y,z), consistently with the approach
taken for the magnetic Eqs. [2]–[5]. In these coordinates, the displacement vector u has the representation
ru¼
1 qðruR Þ 1 quy quz
þ
;
þ
r qr
r qy
qz
and
r 2 uz ¼
8
9
1 q > quz >
1 q2 uz q2 uz
:r
;þ
þ 2
r qr
r 2 qy2
qr
qz
is the usual scalar Laplacian operator. It has been
assumed in Eq. [13] that the current density flows
on both sides of a conducting surface, giving the
factor of two in the first terms on the right-hand
sides of each equation; in addition, the (nonlinear)
self-induced radial components of j B are in opposite directions either side of the current sheet,
and so cancel, leaving only the contribution from
the large axial magnetic field component BZ ¼ BZ0
produced by the permanent magnets in the MRI
equipment.
Concepts in Magnetic Resonance Part B (Magnetic Resonance Engineering) DOI 10.1002/cmr.b
DESIGN OF ASYMMETRIC GRADIENT COILS IN MRI
Equation [11] and its component forms [13] represent a substantially simplified model of the deflection
of the coil. However, these equations nevertheless
still constitute a complicated vector system of equations, and further simplifications are desirable. In particular, we are primarily concerned with the radial
component [13a] of this momentum equation, since
this radial displacement component perpendicular to
the surface of the cylindrical coil will be mostly responsible for the pressure pulses in the air inside the
coil, and hence the noise experienced by the patient.
The azimuthal and axial deflection components [13b]
and [13c] are therefore ignored in this approximation.
Furthermore, it is assumed that the radial component
uR is approximately independent of radial coordinate
r, and that the contribution from the term quy/qy in
Eq. [13a] is small relative to the other terms. Finally,
we make a quasi-equilibrium argument, in which it is
assumed that elastic waves within the coil occur on a
faster time scale than the intervals over which the
current in the gradient coils changes. Under this
approximation, the term q2uR/qt2 in Eq. [13a] is
ignored, so that the coil is essentially regarded as
adjusting rapidly to the changing conditions brought
about by altering the current in the coil windings.
With these additional simplifications, Eq. [13a]
now takes the further approximate form
8 2
9
2
1 q uR q2 uR uR >
>
>
0 ¼ BZ0 jy þ G: 2 2 þ 2 2 >
;:
h
r qy
qz
r
[14]
This expression [14] is taken to hold along the coil
mid-surface rM ¼ a þ h/2 indicated in Fig. 1. Making use of Eq. [8] and allowing for current density
contributions from both the primary and the shield
then gives
1 q2 uR q 2 uR uR
þ 2 2 ¼ D0 ðjPy þ jSy Þ
2
qz
rM
rM
qy2
[15a]
in which it is convenient to define the constant
D0 ¼
4ð1 þ vÞBZ0
:
Eh
[15b]
The two current density components on the right-hand
side of Eq. [15a] are as given in the expressions [4b] and
[5b]. Consequently, Eq. [15] may be used to determine
the radial displacement component uR of the coil.
Pressure Waves
223
movement of the inner coil wall. If rA and pA represent the density and pressure, respectively, of the air
in the coil and vA is the air velocity vector, then the
governing equations are conservation of mass
qrA
þ r ðrA vA Þ ¼ 0;
qt
[16]
conservation of momentum
qvA
1
þ ðvA rÞvA þ
rpA ¼ 0
rA
qt
[17]
and the isentropic gas relation
pA ¼ KrgA :
[18]
The gravitational body force has been ignored in
Eq. [17], since its effect is insignificant over the time
scales at which acoustic effects occur inside the coil.
The quantity K in Eq. [18] is a constant, and the
exponent g ¼ cp/cv ¼ 1.4 is the ratio of specific heats
for air. These well-known nonlinear equations of
motion may be found in the book by Batchelor (22),
for example.
To make analytical progress with modeling the
acoustic effects in the coil, it is again necessary to
linearize Eqs. [16]–[18]. Thus it is assumed that
rA ¼ rA0 þ rA1
vA ¼ 0 þ vA1
pA ¼ pA0 þ pA1 :
[19]
The terms with subscript 0 in Eq. [19] represent the
ambient conditions, and in this paper we have used
the values rA0 ¼ 1.226 kg/m3 and pA0 ¼ 1.013
105 N/m2. The perturbation quantities with subscript
1 are supposed to be small terms, of which quadratic
and higher powers may be ignored. When Eqs. [19]
are substituted into Eq. [16]–[18], the linearized mass
conservation equation
qrA1
þ rA0 r ðvA1 Þ ¼ 0
qt
[20]
is obtained, along with the approximate momentum
equation
qvA1
1
þ
rpA1 ¼ 0
rA0
qt
[21]
and the linearized isentropic relation
To calculate the noise levels within the coil, it is necessary now to consider the response of the air to the
pA1 ¼ ðgpA0 =rA0 ÞrA1 :
[22]
Concepts in Magnetic Resonance Part B (Magnetic Resonance Engineering) DOI 10.1002/cmr.b
224
FORBES ET AL.
These linearized Eqs. [20]–[22] are essentially the
standard acoustic equations assuming small disturbances (23), and may be combined to give a single
wave equation
q2 pA1
¼ c2A0 r2 pA1
qt2
8
9
1 q > qpA1 > 1 q2 pA1 q2 pA1
2
:r
;þ
þ
¼ cA0
r qr
r 2 qy2
qr
qz2
½23
for the pressure perturbation pA1. The constant cA0 is
the isentropic sound speed, and may be calculated
from
c2A0 ¼ g pA0 =rA0 ¼ g RTA0
[24]
after use has been made of the ideal gas law in the
form pA0 ¼ rA0RTA0. The ambient temperature is TA0
and R is the universal gas constant. The Helmholtz
equation used by Shao and Mechefske (24) for modeling acoustic effects in the coil may be derived from
the wave Eq. [23] under the assumption that a purely
sinusoidal time dependence of the perturbation pressure field occurs. In this article, however, more general switching sequences in time will be considered.
Boundary Conditions
Finally, it is necessary to consider the boundary conditions for the air on the inner wall of the coil cylinder. As the coil starts from rest,
pA1 ¼ 0 on r ¼ a at t ¼ 0:
[25a]
Similarly, the air within the coil is initially at rest,
so that Eq. [20] and [22] give
qpA1 =qt ¼ 0 at t ¼ 0:
[25b]
We suppose in addition that a steady-state current
in the coils is established in some switching time tS.
The radial component of the displacement vector is
assumed to follow the time-dependant behavior
uR ðy; z; tÞ ¼ fS ðtÞUR ðy; zÞ;
[26]
in which the quantity fS is a (dimensionless) switching function that starts at zero at time t ¼ 0 and rises
to the value 1 at the switching time tS. For analytical
modeling, it is desirable that fS should have a continuous second derivative for t > 0, and in this article
we have therefore chosen it to be the cubic-spline
function
Figure 2 A sketch of the cubic-spline switching function
fS(t) in Eq. [27], showing the switching time tS.
8
t<0
< 0;
fS ðtÞ ¼ ðt=tS Þ3 3ðt=tS Þ2 þ 3ðt=tS Þ; 0 < t < tS :
:
1;
t > tS
[27]
A sketch of this dimensionless function is given in
Fig. 2. By taking the radial component of the linearized momentum Eq. [21] and combining it with
boundary condition [26], it follows that the perturbed
air pressure in the coil must satisfy the second
boundary condition
qpA1
¼ rA0 fS 00 ðtÞUR ðy; zÞ on r ¼ a:
qr
[28]
Thus, the wave Eq. [23] is to be solved subject to
conditions [25] and [28].
Noise Level
The noise level within the coil may be estimated
from the formula
9
8
jpA1 j>
>
>
>
[29]
SPL ¼ 20 log10 :
;
pref
(see, for example Ref. 23, p 125), in which the reference pressure pref ¼ 2 105 N/m2 is assumed. The
sound pressure level in Eq. [29] is measured in decibels. This quantity (Eq. [29]) is computed over the
interior of the coil, but for ease of understanding, we
show the sound pressure level (Eq. [29]) at the midpoint z ¼ (p þ q)L/2 of the larger target radius r ¼ c1
within the primary coil (in Fig. 1) at y ¼ 0, for the
results in Section 4, to enable simple comparison of
acoustic effects in different coils.
Concepts in Magnetic Resonance Part B (Magnetic Resonance Engineering) DOI 10.1002/cmr.b
DESIGN OF ASYMMETRIC GRADIENT COILS IN MRI
THE SOLUTION PROCESS
The current density and streamfunction on the primary surface r ¼ a and shield surface r ¼ b are given
in Eqs. [4] and [5], respectively. The radial component uR of the coil displacement may now be calculated by solving Eq. [15a]. It is convenient to express
the solution in the form
8
9
N
X
npðz þ LÞ>
>
:
;
UR ðy; zÞ ¼
AM
sin
0n
2L
n¼1
þ
N
M X
X
As yet, however, the coefficients in the representations [4] and [5] are unknown.
At first, it would appear to be possible to determine these coefficients PPmn , and so on, by specifying
a desired target field HZ in Eq. [2] and solving the
first-kind integral equations to obtain the current densities. Indeed, when Eqs. [4] and [5] are substituted
into Eq. [2], after some algebra there results an equation of the form
HZ ðr; y; zÞ ¼
½30
þ
in which the Fourier coefficients may be found from
the relations
D0 ðPP0n þ PS0n Þ
2
½1=rM
AM
mn ¼
BM
mn ¼
2
þ ðnp=2LÞ
D0 ðQPmn þ QSmn Þ
2 þ ðnp=2LÞ2
½ð1 þ m2 Þ=rM
ZL Zp
Umn ðr; z; aÞ ¼
a
p
0
ZL Zp
L
0
M X
N
X
½PSmn cos my
þ QSmn sin myUmn ðr; z; bÞ: ½32
;
:
½31
In this expression, it has been convenient to define
the intermediate functions
½r cos y0 a sinðnpðz0 þ LÞ=ð2LÞÞ
½a2 þ r 2 2ar cos y0 þ ðz0 zÞ2 3=2
dy0 dz0
½r cos y0 a cos my0 cosðnpðz0 þ LÞ=ð2LÞÞ
½a2 þ r 2 2ar cos y0 þ ðz0 zÞ2 3=2
and these integrals may be evaluated numerically
using the composite trapezoidal rule.
The system of Eqs. [32] evaluated on an appropriate mesh of field points nevertheless represents a
highly ill-conditioned system of linear equations for
the Fourier coefficients. As has been observed elsewhere (6–8), it cannot simply be inverted with any
degree of accuracy. Instead, we impose target fields
HTZ (c1; y, z) and HTZ (c2; y, z) on the two inner cylinders shown in Fig. 1, and the additional target field
HTZ (c3; y, z) ¼ 0 on the outer radius r ¼ c3. The aim
is now to minimize the two functionals
Ej ðPPon ; PPmn ; QPmn ; PS0n ; PSmn ; QSmn Þ ¼
ZqLZp
½HTZ ðcj ; y; zÞHZ ðcj ; y; zÞ2 dydz; j ¼ 1; 2 ½34a
pL p
m¼1 n¼1
QPmn sin myUmn ðr; z; aÞ
N
X
PS0n U0n ðr; z; bÞ
n¼1
m¼1 n¼1
2 þ ðnp=2LÞ2
½ð1 þ m2 Þ=rM
L
;
D0 ðPPmn þ PSmn Þ
a
U0n ðr; z; aÞ ¼
p
PP0n U0n ðr; z; aÞ
N
M X
X
½PPmn cos my
m¼1 n¼1
AM
0n ¼
N
X
n¼1
½AM
mn cos my
8
9
npðz þ LÞ>
>
:
;;
þ BM
mm sin my cos
2L
225
dy0 dz0 ;
m 1;
½33
on the inner target zones, along with the additional
functional
E3 ðPPon ; PPmn ; QPmn ; PS0n ; PSmn ; QSmn Þ ¼
Z L Zp
HZ2 ðc3 ; y; zÞdydz
½34b
L p
on the outer target zone, subject to the additional
requirement that certain other penalty functions
should be minimized simultaneously. The axial magnetic field component HZ in these expressions [34] is
calculated from Eq. [32].
It is necessary to comment briefly on the condition
[34b] used to reduce the strength of the magnetic
field in the region r > c3 beyond the shielding coil.
Here, only the axial component of the field is taken
Concepts in Magnetic Resonance Part B (Magnetic Resonance Engineering) DOI 10.1002/cmr.b
226
FORBES ET AL.
into account explicitly, although it is also necessary
for the shields to reduce the transverse components
in order for shielding to be as effective as possible. It
is possible to include all three field components in
the integrand in Eq. [34b], although at the expense of
significant additional calculation. We have found,
however, that the axial components of the field are
also reduced substantially by minimizing the quantity
[34b] in comparison with the corresponding
unshielded coil. This is discussed in more detail in
(8), where it was shown that reduction of the total
field could be done very effectively using shields that
are substantially longer than the primary, although of
course with the disadvantage of adding significantly
to the effects of claustrophobia.
In the attempted minimization of acoustic noise
within the coil, an obvious choice for a penalty function is
F1 ðPP0n ; PPmn ; QPmn ; PS0n ; PSmn ; QSmn Þ ¼
ZL Zp
u2R ðy; zÞdy dz;
½35
L p
in which the coil deflection function is as given in
Eq. [30], with coefficients [31]. Following Forbes
and Crozier (8), we also allow the option of additional streamline curvature functionals
FP ðPP0n ; PPmn ; QPmn Þ
¼
ZL Z p
¼
qR
¼ 0;
qPPjk
qR
¼0
qQPjk
qR
¼ 0;
qPS0k
qR
¼ 0;
qPSjk
qR
¼ 0;
qQSjk
j ¼ 1; . . . ; M:
½39
2
jr cP ðy; zÞj dy dz [36]
on the primary coil r ¼ a, and
ZL Zp
qR
¼ 0;
qPP0k
k ¼ 1; . . . ; N;
2
L p
FS ðPS0n ; PSmn ; QSmn Þ
that here their value is not known in advance, and
must be determined empirically. By minimizing the
residual R in Eq. [38], then each of the positive-definite terms in Eq. [34]–[37] is also minimized, and the
constants l1, lP, and lS should be chosen large
enough to make the problem acceptably well conditioned, but still small enough to achieve the original
aim of minimizing the squared-error terms in Eqs.
[34]. This process is known as Tikhonov regularization (9, 10). It is usually sufficient to choose the
smoothness coefficients lP and lS to be of order
1012. For materials appropriate to the construction
of the coil former, the Young’s modulus is approximately E ¼ 1.3 1010 N/m2 and the Poisson ratio is
of order v ¼ 0.2, so that D0 in Eq. [15b] has a small
numerical value; consequently, constant l1 can be
taken to be large, so that the terms in Eq. [38] are all
in balance. We have obtained good designs with this
parameter as large as l1 ¼ 1017, as will be discussed
in Section 4.
The residual error (Eq. [38]) is minimized with
respect to the unknown Fourier coefficients by
requiring that
jr2 cS ðy; zÞj2 dy dz [37]
L p
on the shield r ¼ b. Minimizing these expressions
[36] and [37] is equivalent to smoothing the winding
pattern designs on the primary and shield. The two
streamfunctions in these expressions are given in
Eqs. [4c] and [5c].
To create a well-conditioned system that can be
solved for the Fourier coefficients, we combine Eqs.
[34]–[37] to create the total residual error function
R ¼ E1 þ E 2 þ E 3 þ l 1 F1 þ l P F P þ l S F S :
[38]
The three constants l1, lP, and lS in the expression (38) play a role similar to that of Lagrange multipliers in constrained minimization problems, except
Each of these represents a linear system of equations
involving the Fourier coefficients (because the original functionals [34]–[37] are quadratic in these quantities). Thus a matrix equation can be developed, from
which the coefficients PPjk , and so on, can be found. In
fact, there are very significant extra efficiencies available by exploiting the structure of these Eq. [39]. It
turns out that the equations are strongly de-coupled,
so that it is possible to compute the 2N coefficients
PPjk , PSjk k ¼ 1, . . . , N independently at each value of
index j. Similarly, the 2N coefficients QPjk , QSjk k ¼
1, . . . , N can likewise be determined separately for
each j. The details are similar to those presented by
Forbes and Crozier (8) and so will not be repeated
here. The de-coupled structure of the system (Eq.
[39]) means that the present algorithm for obtaining
these Fourier coefficients has optimal efficiency.
Once these coefficients have been found, it is then
straightforward to evaluate expressions [26] and [30]
for the deflection of the coil in response to Lorentz
forces. Finally, it is necessary to obtain the acoustic
pressure perturbation pA1 within the coil, by solving
the wave Eq. [23], subject to boundary conditions
Concepts in Magnetic Resonance Part B (Magnetic Resonance Engineering) DOI 10.1002/cmr.b
DESIGN OF ASYMMETRIC GRADIENT COILS IN MRI
[25] and [28] on the inner wall of the coil. This is
somewhat difficult, but can be done with the use of
Laplace transforms in time. We suppose that the
pressure perturbation has the representation
þ
tively the modified Bessel functions of the first and
second kind, of order zero (25, p 374). Here, it is
convenient to define
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
8np92 8 s 92
> ;
> :
x ¼ : ; þ:
2L
cA0
8
9
N
X
npðz þ LÞ>
>
:
;
pA1 ðr; y; z; tÞ ¼
F0n ðr; tÞ sin
2L
n¼1
N
M X
X
½Emn ðr; tÞ cos my
m¼1 n¼1
8
9
npðz þ LÞ>
:
;: ½40
þ Fmn ðr; tÞ sin my cos>
2L
in 0 , r , a ½41a
subject to
on
r ¼ a:
Zt
fs 00 ðtÞK0n ðr; t tÞdt;
[47]
0
[42]
The sound speed in Eq. [41a] is given by the
expression [24], and the coefficient AM
0n in Eq. [42] is
defined in the relation [31]. The Laplace transform of
the desired function is defined in the usual way as
F~0n ðr; sÞ ¼
Z1
est F0n ðr; tÞdt
[43]
0
and we denote the transform of the second derivative
of the switching function in Eq. [27] as
est fs 00 ðtÞdt
0
¼
rA0 AM
0n
in which the kernel function K0n is given by the
inverse Laplace transform (26, p 406) as
qF0n
¼ rA0 fS00 ðtÞAM
0n
qr
Z1
[46]
in which D(s) is the transform defined in [44].
According to the convolution theorem for Laplace
transforms (26, p 393), the original function F0n(r,t)
may now be recovered from Eq. [46], and becomes
[41b]
and the condition
DðsÞ ¼
DðsÞI0 ðrxÞ
F~0n ðr; sÞ ¼ rA0 AM
;
0n
xI1 ðaxÞ
F0n ðr; tÞ ¼
F0n ¼ 0 and qF0n =qt ¼ 0 at t ¼ 0
[45]
The second-kind function K0, however, must be
deleted from the solution since it becomes
unbounded as r ? 0. The boundary condition [42]
then yields the transform [43] of the solution in the
form
This form [40] is substituted into the wave Eq.
[23], and yields partial differential equations for each
coefficient function.
The first coefficient in Eq. [40], for example, satisfies
8
9
8np92
1 q2 F0n
1 q > qF0n >
:
;
:
;
r
F0n þ
¼
2L
r qr
qr
c2A0 qt2
227
8
9
6 >1
1
t
: expðsts Þ s >
;
t3s s2 s2
s
½44
When the Laplace transform of Eq. [41a] is taken,
using the initial conditions [41b], it becomes evident
that the function F~0n (r; s) in Eq. [43] involves the
two functions I0 (rx) and K0 (rx), which are respec-
1
K0n ðr; tÞ ¼
2pi
cþi1
Z
est
ci1
I0 ðrxÞ
ds
xI1 ðaxÞ
[48]
and x is the function of the transform variable s
defined in Eq. [45]. The path of integration in Eq.
[48] consist of a vertical line in the complex s-plane,
as sketched in Fig. 3, on some path Re{s} ¼ c that lies
to the right of all the singularities of the integrand.
Although the integral in Eq. [48] notionally involves the function x, which is defined as a fractional
power involving s, it may nevertheless be shown that
the integrand does not possess branch-type singularities. There are only simple pole singularities present
in the integrand of Eq. [48], and these occur at the infiðnÞ
ðnÞ
nite number of locations s ¼ 6ib0;1 and s ¼ 6ib0;k ,
k ¼ 2, 3, . . . , in which we have written
ðnÞ
np
and
2Lffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
s
8np92 8j 92
: 1;k >
; ;
¼ cA0 : ; þ>
2L
a
b0;1 ¼ cA0
ðnÞ
b0;k
k ¼ 2; 3; . . .
½49
Concepts in Magnetic Resonance Part B (Magnetic Resonance Engineering) DOI 10.1002/cmr.b
228
FORBES ET AL.
Equation [47] now gives an expression for the
original function F0n(r,t).
A similar solution process yields the remaining
Fourier coefficient functions in Eq. [40]. The convolution theorem of Laplace transforms again gives
Rt 00
Emn ðr; tÞ ¼ rA0 AM
mn fS ðtÞKmn ðr; t tÞdt
0
Fmn ðr; tÞ ¼ rA0 BM
mn
Rt
fS00 ðtÞKmn ðr; t tÞdt ½51
0
with kernel function
2c2A0
a
1
X
Jm ðr=aÞj0 m;k
ðnÞ
h
i
sinðbm;k tÞ: ½52
ðnÞ
2
0
0
k¼2 bm;k 1 ðm=j m;k Þ Jm ðj m;k Þ
Kmn ðr; tÞ ¼
Figure 3 Contour in the complex transform s-plane, for
evaluating the inverse Laplace transform. Simple poles
occur along the imaginary axis.
The symbol j1,k in Eq. [49] refers to the k-th zero
of the Bessel function J1(z), following the notation in
Abramowitz and Stegun (25, p 370). Thus the poles
occur on the imaginary axis in the s-plane, as
sketched in Fig. 3. The integral in Eq. [48] is therefore evaluated by summing the residues at each of
the poles, to give the infinite series
K0n ðr; tÞ ¼
"
#
ðnÞ
1
2c2A0 sinðb0;1 tÞ X
J0 ððr=aÞj1;k Þ
ðnÞ
þ
sinðb0;k tÞ : ½50
ðnÞ
ðnÞ
a
b
k¼2 b J0 ðj1;k Þ
0;1
pA1 ðr; y; z; tÞ ¼ 2rA0
0;k
This result is again obtained using a contour similar to that shown in Fig. 3. In this expression, the
quantities
ðnÞ
bm;k
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
8np92 8j0 92
> m;k ;
>
¼ cA0 : ; þ:
2L
a
[53]
have been defined, similarly to Eqs. [49]. The nota0
in Eq. [52] and [53] refers to the k-th zero
tion jm,k
of the derivative Jm0 (z) of the first-kind Bessel
function of order m, again following the notation
of (25).
The final form of the perturbation pressure pA1 is
obtained by combining Eqs. [40], [47], and [50]–[52]
to yield
8
9
N
c2A0 X
AM
npðz þ LÞ>
ðnÞ
0n
>
:
;
Sðt;
b
;
t
Þ
sin
S
0;1
2L
a n¼1 bðnÞ
0;1
ðnÞ
8
9
1 J ðr=aÞj
N
X
c2A0 X
npðz þ LÞ>
0
1;k Sðt; b0;k ; tS Þ
M
>
:
;
2rA0
A0n
sin
ðnÞ
2L
a n¼1
b0;k J0 ðj1;k Þ
k¼2
8
9
N
M X
c2A0 X
npðz þ LÞ>
M
M
>
:
;
2rA0
ðA cos my þ Bmn sin myÞ cos
2L
a m¼1 n¼1 mn
ðnÞ
1
X
Jm ððr=aÞj0 m;k ÞSðt; bm;k ; tS Þ
h
i
: ½54
ðnÞ
2
k¼1 bm;k 1 ðm=j0 m;k Þ Jm ðj0 m;k Þ
Concepts in Magnetic Resonance Part B (Magnetic Resonance Engineering) DOI 10.1002/cmr.b
DESIGN OF ASYMMETRIC GRADIENT COILS IN MRI
In this expression, it has proven convenient to define
the quantity
Sðt; b; tS Þ ¼
Zt
fS00 ðtÞ sinðbðt tÞÞdt:
0
This integral can be evaluated easily, and takes the
value
229
x-Gradient Coils
To begin, we consider an x-gradient coil, so that in
polar coordinates, the two inner target fields are
taken to be
HTZ ðc1 ; y; zÞ ¼ Hmax cos y
HTZ ðc2 ; y; zÞ ¼ Hmax ðc2 =c1 Þ cos y;
[55]
for t > tS :
in which Hmax ¼ Bmax/m0, which here has the value
1.59 104 A/m.
Winding patterns for both the primary and the
shield are shown in Figs. 4(a,b), respectively, for a
coil designed with the choice of Tikhonov regularization parameter l1 ¼ 1 in Eq. [38]. The dashed lines
in each picture correspond to windings in which the
current is negative. Because of the highly asymmetri-
In this Section, results are presented for gradient
coils of half-length L ¼ 0.5 m, inner radius a ¼ 0.3
m, and outer radius b ¼ 0.4 m; the thickness of the
coil material is therefore h ¼ 0.1 m. Its Young’s
modulus is taken to be E ¼ 1.3 1010 N/m2 and
Poisson’s ratio for the coil is v ¼ 0.2. The background magnetic field (in which the gradient coils
are immersed) is assumed here to produce the constant axial field BZ0 ¼ 2 T and the gradient coil is
taken to generate a maximum induction field of
Bmax ¼ 0.02 T. These parameters have been chosen
to correspond to actual coil geometries. The inner
target region (illustrated in Fig. 1) is assumed to be
positioned very asymmetrically, with parameters
p ¼ 0.7 and q ¼ 0.1, and the two inner target radii
are c1 ¼ 0.2 m and c2 ¼ 0.1 m. The outer target radius is generally taken to be c3 ¼ (3/2)b, where b is
the outer (shielding) radius of the coil, although this
value is sometimes varied to give better quality
fields as will be seen presently, particularly when zgradient coils are to be designed. The minimization
algorithm [39] was implemented in the FORTRAN
language, and took about 25 min run time to produce a shielded design, using an AMD Athlon
1600þ CPU computer with 256 MB RAM memory.
A significant advantage, in terms of computer run
time, is gained by caching the functions Umn in Eq.
[33] so that they are computed only once and then
stored, and this has been done in our routines. For
all the results shown in this Section, the two
smoothness coefficients lP and lS in the residual
Eq. [38] were chosen to be 1012.
Figure 4 Winding patterns for an x-gradient coil with
inner and outer radii a ¼ 0.3 m and b ¼ 0.4 m, for (a) the
primary coil and (b) the shield coil. This system has been
designed with regularizing parameter l1 ¼ 1. Dashed lines
indicate reversed windings.
6 t
sin bt
1 þ cos bt
Sðt; b; tS Þ ¼ 2
btS
btS tS
Sðt; b; tS Þ ¼
6
bt2S
for t < tS
sinðbðt tS ÞÞ sin bt
cos bt þ
btS
PRESENTATION OF RESULTS
Concepts in Magnetic Resonance Part B (Magnetic Resonance Engineering) DOI 10.1002/cmr.b
230
FORBES ET AL.
cally located inner target regions, the windings are
consequently asymmetrical in the axial z coordinate,
and in fact, there is a small reverse winding present
at the opposite end to the target regions, in the case
of the primary coil shown in Fig. 4(a). The horizontal
axis in each diagram corresponds to the azimuthal
angle y around the circumference of the coil, so that
each winding pattern in Fig. 4 is to be thought of as
wrapped around the cylinder at the appropriate radius. Thus each winding pattern resembles a Golay
coil [see Jin (1)], in which the direction of the current
in the shield opposes that in the primary. The relatively small value of the parameter l1 ¼ 1 used in
designing this coil represents a situation in which
there is essentially no attempt to minimize the overall
coil deflections that lead to acoustic noise. The winding patterns in Fig. 4 therefore correspond to coils in
which noise effects are ignored, and resemble quite
closely the designs obtained for gradient coils in the
algorithm of Forbes and Crozier (8).
By contrast, Figs. 5(a,b) show the primary and
shield windings for an x-gradient coil designed with
parameter l1 ¼ 1017. This value was chosen after careful experimentation, and represents about the optimum
trade-off between being able to reproduce the target
fields Eq. [55] faithfully on the one hand, and reducing
the acoustic noise on the other. The winding patterns
for this new coil are very different to those in Fig. 4,
and it is clear that, at least in this instance, the ability
to minimize acoustic noise within the coil is achieved
by a large number of reverse windings on both the primary and the shield, particularly at the top of each picture, in the portion of the coil far from the target zone
(which is located over the region 0.35 < z < 0.05
m). It is possible, however, that such a coil with a large
number of reverse windings in close proximity may be
inefficient in terms of its power use, and may also pose
problems from the point of view of local Lorentz
forces between different sections of the same coil.
For each coil design, we have evaluated the perturbation pressure pA1 from the solution (Eq. [54]),
using the MATLAB programming environment
(which has library functions for evaluating Bessel
functions). The coefficients PPmn , and so on, produced
by the FORTRAN implementation of algorithm [39]
are read into the MATLAB routine, and used to create the coefficients AM
mn , and so on, in Eq. [31]. These
are then used to evaluate expression [54]. The major
practical difficulty to be overcome in implementing
Eq. [54] is to generate a table of values of the zeros
0
j1,k and jm,k
of the first-kind Bessel functions and
their derivatives. This task has received significant
attention in the literature, and includes techniques
such as the use of a high-order Newton’s method
Figure 5 Winding patterns for an x-gradient coil with
inner and outer radii a ¼ 0.3 m and b ¼ 0.4 m, for (a)
the primary coil and (b) the shield coil. This system has
been designed with regularizing parameter l1 ¼ 1017.
Dashed lines indicate reversed windings.
(27), a tri-diagonal matrix approach (28) and global
Newton-like methods (29). We have investigated a
number of these techniques, but it is interesting and
perhaps surprising that the most reliable results were
obtained simply by using the bisection algorithm [see
Atkinson (30, p 43)]. Good estimates for the zeros of
Bessel functions and their derivatives are available
from Abramowitz and Stegun (25, p 371) and from
these it is reasonably straightforward to construct
intervals that enclose each zero, thus guaranteeing
that the simple bisection method will indeed converge
to the desired root in each case. Once the perturbation
pressure pA1 has thus been evaluated, the sound pressure level is then calculated using Eq. [29].
Figure 6 gives a comparison between the noise levels produced by the two different coils designs shown
in Figs. 4 and 5, obtained with regularization parame-
Concepts in Magnetic Resonance Part B (Magnetic Resonance Engineering) DOI 10.1002/cmr.b
DESIGN OF ASYMMETRIC GRADIENT COILS IN MRI
Figure 6 Sound pressure level (in decibels) as a function
of time, for the un-optimized coil in Fig. 4 (dashed line)
and the optimized coil in Fig. 5 (solid line). The switching time is 0.25 ms.
ters l1 ¼ 1 and l1 ¼ 1017. For convenience of comparison, sound pressure levels are shown at the point (r, y,
z) ¼ (c1, 0, (p þ q) L/2), which is located at the midpoint on the outer surface of the inner target regions
shown in Fig. 1. In each case, the switching time in Eq.
[27] is tS ¼ 2.5 104 s (0.25 ms). Time t is shown on
the horizontal axis, and the vertical scale gives the
sound pressure level in decibels. For the coil design
shown in Fig. 4, in which acoustic effects are essentially not taken into account, the noise history is drawn
with a dashed line and indicates a peak level of about
124 dB. The solid line corresponds to the acoustic
noise produced by the optimized coil in Fig. 5, for
which l1 ¼ 1017, and shows peak noise levels in the
first pulse of about 116 dB. Thus the optimized coil of
Fig. 5, while giving very differently shaped winding
patterns to those in Fig. 4, nevertheless only generates
a reduction of about 8 dB in overall noise level.
The essential requirement for an x-gradient coil is,
of course, that the field it generates should match Eq.
[55] to a high degree of accuracy in the inner target
zones. For a shielded coil, the field should also die
away as rapidly as possible beyond the coil. To check
that these key requirements are satisfied, the axial
field components HZ generated by the coils in Figs. 4
and 5 are shown in Figs. 7(a,b) respectively, on the
mid-plane z ¼ (p þ q)L/2 of the inner target regions,
along the line y ¼ 0. In addition, the target field (Eq.
[55]) is also drawn on each diagram, and appears as
the diagonal dashed line appropriately labeled. The
inner target region c1 < x < c1 is indicated in Figs.
7(a,b) with vertical dashed lines.
For the shielded coil of Fig. 4, for which noise
generation is effectively not penalized, Fig. 7(a)
231
shows that the target field (Eq. [55]) is indeed
matched extremely closely over the entire target
region and significantly beyond it. Likewise, the field
exterior to the coil falls away to zero very rapidly.
The portion of the field shown in Fig. 7(b), corresponding to the deflection-minimized winding pattern
of Fig. 5, also indicates that this coil is capable of
matching the target field (Eq. [55]) in the interior and
is able to produce only small fields externally. However, it is clearly not as precise as in Fig. 7(a), and
thus it is evident that the minimization of acoustic
noise using winding patterns alone comes at the cost
of slightly poorer quality magnetic fields.
Although Figs. 7(a,b) have confirmed that coils
designed here are indeed capable of replicating the target fields accurately in the required regions, Fig. 6
Figure 7 The axial field computed on a transverse line
y ¼ 0, z ¼ (p þ q)L/2 through the centre of the inner zones,
for (a) the regular shielded coil of Fig. 4 and (b) the deflection-minimized coil of Fig. 5. The computed field is shown
with a solid line, and the target field with a dashed line. The
edges x ¼ 6c1 of the target region are also indicated.
Concepts in Magnetic Resonance Part B (Magnetic Resonance Engineering) DOI 10.1002/cmr.b
232
FORBES ET AL.
Figure 8 The coil shape in response to Lorentz forces,
for the x-gradient coil in Figure 5. Displacements have
been magnified by a factor of 104 for ease of viewing.
[Color figure can be viewed in the online issue, which is
available at www.interscience.wiley.com.]
showed reductions in noise levels inside the coil of
about only 8 dB. We have experimented in detail with
coil designs at higher values of the regularizing parameter l1 and we have observed that noise levels can, in
fact, be reduced as much as desired, simply by increasing l1 as much as needed. However, the ability of the
coil to produce a field that matches the target field (Eq.
[55]) is seriously compromised by this process. Consequently, the value l1 ¼ 1017 used in Figs. 5–7 is about
optimum for the x-gradient coil investigated here.
Following Mechefske et al. (16) and Yao et al. (31),
we show in Fig. 8 the deflection caused to the coil by
Lorentz forces, for the design in Figs. 5(a,b). The radial deflection function uR has been computed from
Eq. [30], and has been amplified by a factor of 104 so
as to make it visible. The shape drawn in Fig. 8 is thus
the surface r ¼ a þ 104uR, for some time t > tS at
which the function fS in Eq. [27] has the value one.
Clearly the greatest deflections occur in the asymmetrically located target region (at the left end of the coil),
as a result of which the coil undergoes an overall lateral displacement towards the negative x-direction.
Smaller deflections also occur at the end of the coil
further from the target region, but these appear to be
more circularly symmetric. The small undular deflections on the coil surface are indications of the patterns
of reverse windings that feature in Fig. 5.
The results for these shielded coils have also been
compared against unshielded designs. Consistently
with previous work (6–8), we have observed that the
performance of the coil, in terms of its ability to match
the target field faithfully, is actually enhanced slightly
by the presence of the shields. In addition, shields con-
tribute to a slight reduction in the noise level within
the coil. For example, the optimized (shielded) coil
with l1 ¼ 1017 in Fig. 6 generates a peak noise level of
about 116 dB, but when shields are removed, this
increases to 119 dB. The overall shape of the noise history profile is very similar to that of the optimized coil
in Fig. 6, and so is not shown here.
The winding pattern of the unshielded coil is displayed in Fig. 9, at the same values of the coil design
parameters as for Fig. 5(a). The two patterns are
quite similar over the target region at the bottom of
each diagram, but the shielded coil in Fig. 5(a) possesses a region of more strongly alternating current
directions at the top end of the diagram, far from the
target zone. These are necessary in order for the field
external to the shielded coil to drop more rapidly to
zero, and are additionally responsible for the reduced
overall deflection and consequent slight reduction in
noise level. They may, however, constitute a difficulty in terms of practical manufacture of the coil.
z-Gradient Coils
We conclude this presentation of results with a brief
investigation of the behavior of the z-gradient coil.
This device is designed to produce a linear field in
the axial direction, so that the target fields on the two
inner radii have the forms
HTZ ðc1 ; y; zÞ ¼ HTZ ðc2 ; y; zÞ
2Hmax h z p þ qi
¼
: ½56
2
ðq pÞ L
Figure 9 Winding pattern for an unshielded x-gradient
coil with inner and outer radii a ¼ 0.3 m and b ¼ 0.4 m.
This system has been designed with regularizing parameter l1 ¼ 1017. Dashed lines indicate reversed windings.
Concepts in Magnetic Resonance Part B (Magnetic Resonance Engineering) DOI 10.1002/cmr.b
DESIGN OF ASYMMETRIC GRADIENT COILS IN MRI
As indicated by Forbes and Crozier (8), these zonal
coils require less computational effort to design using
this strategy, since they only require calculation of
the coefficients PP0n and PS0n in Eqs. [4] and [5], and
all the higher-order (tesseral) coefficients are zero.
This means, in particular, that the axial current-density components jPZ and jSZ in Eqs. [4a] and [5a] are
identically zero.
The axial field HZ computed along the coil axis
x ¼ y ¼ 0 is contrasted for an unshielded coil in Fig.
10(a) and a shielded coil in Fig. 10(b). In both cases,
the regularizing parameter for the deflection was
again taken to be l1 ¼ 1017. Unlike the case of the xgradient coil shown in Fig. 7, it turns out that the
type of shielding imposed here does degrade slightly
the ability of the z-gradient coil to match the target
Figure 10 A comparison of the axial field on the centreline of a z-gradient coil for (a) an unshielded coil and (b)
a shielded coil. Both results have been obtained with regularizing parameter l1 ¼ 1017. The vertical dashed lines
in each picture indicate the positions of the coil ends and
the asymmetrically located target zone.
233
Figure 11 Sound pressure level (in decibels) as a function of time, for the unshielded z-gradient coil (dashed
line), the pure shielded coil (dashed-dot line, l1 ¼ 1) and
the deflection-minimized shielded coil (solid line, l1 ¼
1017). The switching time is 0.25 ms.
field (Eq. [56]) accurately, as is evident from Fig. 10.
Nevertheless, the unshielded coil in Fig. 10(a) clearly
generates much larger fields at the ends of the coil,
which is an undesirable feature in that case. The
shielding condition of minimizing exterior fields was
imposed on the cylindrical surface c3 ¼ 1.3b in Fig.
10(b), since it is found that moving this external surface closer to the shield windings improves the quality of the interior magnetic fields.
Figure 11 compares the sound pressure levels produced inside the shielded coil, the deflection-minimized shielded coil obtained with l1 ¼ 1017 and the
unshielded coil also obtained with l1 ¼ 1017. The
optimized coil is able to reduce the noise level by
about 5 dB, when compared with the purely shielded
coil (l1 ¼ 1), although an examination of the winding
patterns indicates it does this by the use of several coil
sections with reversed currents. Interestingly, the
unshielded coil creates significantly less noise within
the coil, as is evident from Fig. 11, although of course
it produces a field which decays only slowly exterior
to the coil. The winding patterns corresponding to the
three curves in Fig. 11 are not shown here, as they
consist simply of circular strips around the cylindrical
former, for these z–gradient coils.
To conclude this presentation of results, the
deflections corresponding to each of the three curves
in Figure 11 are shown in Figs. 12(a–c). The displacements are again magnified by a factor of 104 for
ease of viewing, so that the surfaces shown represent
graphs of r ¼ a þ 104uR. Each coil retains circular
symmetry, as the windings are simply arranged circularly around the cylinder. The unshielded coil is
Concepts in Magnetic Resonance Part B (Magnetic Resonance Engineering) DOI 10.1002/cmr.b
234
FORBES ET AL.
shown in Fig. 12(a), and it is characterized by fairly
small displacements over most of the coil, except at
the end furthest from the target zone, where the
deflections are moderately large. This is consistent
with the result shown in Fig. 10(a), for which the
field becomes large at the far end z ¼ 0.5 m. This
coil produces the lowest noise levels in Fig. 11, but
of course, it cannot prevent field leakage beyond the
coil, or equivalently, it cannot protect the interior
field from stray influences external to the coil.
The case of a purely shielded coil that is uncorrected for overall deflection is shown in Fig. 12(b). It
is essentially the type of coil designed with the algorithm of Forbes and Crozier (8). It is evident that
such a coil allows large deflections at the end furthest
from the target zone, and thus is expected to be the
noisiest of the three designs. This is confirmed by the
results of Fig. 11. It must be recalled, however, that
the radial deflections have been multiplied by a factor of 104, so that Fig. 12(b) represents a severe exaggeration of the actual shape.
The deflection-minimized coil, obtained with regularization parameter l1 ¼ 1017, is displayed in Fig.
12(c). It has the smallest overall deflection of the
three coils shown, and thus produces the smallest
amount of noise, as shown in Fig. 11. Nevertheless,
it possesses more undulations in its surface, corresponding to the regions of reversed current in the coil
windings.
DISCUSSION AND CONCLUSIONS
Figure 12 The coil shape in response to Lorentz forces,
for the z-gradient coils in Fig. 11. The diagrams show (a)
the unshielded coil, (b) the pure shielded coil (l1 ¼ 1) and
(c) the deflection-minimized shielded coil (l1 ¼ 1017). Displacements have been magnified by a factor of 104 for ease
of viewing. [Color figure can be viewed in the online issue,
which is available at www.interscience.wiley.com.]
A method has been presented in this article for computing the deflection of a cylindrical gradient coil
because of Lorentz forces, and for incorporating this
into an optimization strategy for coil design. The
technique uses Tikhonov regularization, as in previous studies (6–8), to overcome the ill-conditioned
mathematical structure of the problem. In addition, a
closed-form solution for the acoustic noise within the
coil has been presented, and by designing coils to
minimize overall deflection, the noise level is thereby
reduced.
An approximate Eq. [15] has been derived for calculating the deflection of the coil. It has the advantage of being a relatively simple linear partial differential equation that can be solved in a reasonably
straightforward manner, once a form for the current
density in the coil has been chosen. Equation [15]
approximates the coil deflection as being mainly radial, although the results, such as in Fig. 8, suggest
that other displacement modes might also usefully be
considered in a more complete model of the deforma-
Concepts in Magnetic Resonance Part B (Magnetic Resonance Engineering) DOI 10.1002/cmr.b
DESIGN OF ASYMMETRIC GRADIENT COILS IN MRI
tion. This could be done, in the present quasi-steady
context, by solving an approximation to the full linearized steady-state Eqs. [13a]–[13c] with the time
derivatives removed, although significantly greater
complexity in the mathematics would occur. This
approach is presently being investigated, however. A
fuller theory of this type, in which all three components uR, uy, and uZ of the displacement vector in Eq.
[12] are accounted for, may also allow the possibility
of considering the precise manner in which the coil
has been tethered in the MRI equipment; this detail
has been overlooked in the present treatment based
on the simplified Eq. [15].
It has been assumed in the present article that the
coil is unloaded, so that the presence of a patient
within the primary has not been considered. This
could likewise be accounted for in a more complete
theory, although at the expense of considerably more
difficult mathematics. Indeed, it is likely that the
present advantage of a closed-form mathematical solution that can be incorporated into optimization
strategies, would be lost, and finite difference or finite element methods would be needed instead. It is
to be expected, however, that the presence of a
patient will modify the sound pressure level within
the coil.
The results presented here show that the noise
within the coil may be reduced by shielding, and
may be reduced further by redesigning the winding
patterns using the technique in this article. There is,
however, a trade-off between acoustic noise reduction and the accurate matching of the target field, so
that a practical limit exists to the amount by which
noise can be reduced within conventional coil geometry. Indeed, the present article possibly indicates the
limits of noise reduction that can be achieved by
optimizing coil windings alone. Further gains will at
least involve a more complete analysis of the full
deflection possibilities for the coil and the details of
how it is tethered, but may ultimately require nonconventional coil geometries.
In this analysis, it has been assumed that the coil
moves rapidly to a quasi steady-state configuration,
and that the noise within it is therefore caused by the
primary pressure pulse associated with its deformation. It is possible to account for the unsteady nature
of the coil deflection more fully, by including the
time-dependent terms in the Eqs. [13a]–[13c], and a
more realistic (periodic) switching sequence than that
shown in Fig. 2. This, however, involves considerably more complicated mathematics than used here,
and may in fact lose the advantages of the analytical
approach adopted in this article, and as such is
beyond the scope of the present work.
235
ACKNOWLEDGMENTS
This work is supported in part by Australian
Research Council grant DP0343350.
REFERENCES
1. Jin J. 1999. Electromagnetic analysis and design in
Magnetic Resonance Imaging (Biomedical Engineering Series). Boca Raton: CRC Press.
2. Turner R. 1986. A target field approach to optimal
coil design. J Phys D: Appl Phys 19:147–151.
3. Turner R. 1993. Gradient coil design: a review of
methods. Magn Reson Imaging 11:903–920.
4. Turner R. 1994. Electrical coils. U.S. Patent 5, 289, 151.
5. Crozier S, Doddrell DM. 1993. Gradient-coil design
by simulated annealing. J Magn Reson Ser A
103:354–357.
6. Forbes LK, Crozier S. 2001. A novel target-field
method for finite-length magnetic resonance shim
coils, Part 1: Zonal shims. J Phys D Appl Phys 34:
3447–3455.
7. Forbes LK, Crozier S. 2002. A novel target-field method
for finite-length magnetic resonance shim coils, Part 2:
Tesseral shims. J Phys D: Appl Phys 35:839–849.
8. Forbes LK, Crozier S. 2003. A novel target-field
method for magnetic resonance shim coils, Part 3:
shielded zonal and tesseral coils. J Phys D: Appl Phys
36:68–80.
9. Delves LM, Mohamed JL. 1985. Computational
Methods for Integral Equations. Cambridge: Cambridge University Press.
10. Neittaanmäki P, Rudnicki M, Savini A. 1996. Inverse
Problems and Optimal Design in Electricity and
Magnetism (Monographs in Electrical and Electronic
Engineering, Vol. 35). Oxford: Clarendon Press.
11. Forbes LK, Crozier S. 2004. Novel target-field method
for designing shielded bi-planar shim and gradient
coils. IEEE Trans Magn 40:1929–1938.
12. Fishbain D, Golberg M, Labbe E. 1988. Long-term
claustrophobia following MRI. Am J Psychol 145:
1038–1039.
13. Jackson JD. 1999. Classical Electrodynamics, 3rd ed.
New York: Wiley.
14. Chapman BLW, Mansfield P. 1995. A quiet gradientcoil set employing optimized, force-shielded, distributed coil designs. J Magn Reson Ser B 107:152–157.
15. Mansfield P, Haywood B, Coxon R. 2001. Active
acoustic control in gradient coils for MRI. Magn
Reson Med 46:807–818.
16. Mechefske CK, Yao G, Li W, Gazdzinski C, Rutt
BK. 2004. Modal analysis and acoustic noise characterization of a 4T MRI gradient coil insert. Concepts
Magn Reson Part B. 22B:37–49.
17. Shao W, Mechefske CK. 2005. Acoustic analysis of a
gradient coil winding in an MRI scanner. Concepts
Magn Reson B Magn Reson Eng 24:15–27.
Concepts in Magnetic Resonance Part B (Magnetic Resonance Engineering) DOI 10.1002/cmr.b
236
FORBES ET AL.
18. Shao W, Mechefske CK. 2005. Analysis of the sound
field in finite length infinite baffled cylindrical ducts
with vibrating walls of finite impedance. J Acoust
Soc Am 117:1728–1736.
19. Brideson MA, Forbes LK, Crozier S. 2002. Determining complicated winding patterns for shim coils using
stream functions and the target-field method. Concepts Magn Reson 14:9–18.
20. Aris R. 1962. Vectors, Tensors and the Basic Equations of Fluid Mechanics. New York: Dover.
21. Boresi AP, Chong KP. 1987. Elasticity in Engineering Mechanics. New York: Elsevier.
22. Batchelor GK. 1977. An Introduction to Fluid Dynamics. Cambridge: Cambridge University Press.
23. Kinsler LE, Frey AR. 1962. Fundamentals of Acoustics, 2nd ed. New York: Wiley.
24. Shao W, Mechefske CK. 2005. Analyses of radiation
impedances of finite cylindrical ducts. J Sound Vib
286:363–381.
25. Abramowitz M, Stegun IA, eds. 1972. Handbook of
Mathematical Functions. New York: Dover.
26. Marsden JE. 1973. Basic Complex Analysis. San
Francisco: W.H. Freeman.
27. Temme NM. 1979. An algorithm with ALGOL 60
program for the computation of the zeros of ordinary
Bessel functions and those of their derivatives.
J Comput Phys 32:270–279.
28. Ikebe Y, Kikuchi Y, Fujishiro I. 1991. Computing zeros and orders of Bessel functions. J Comput Appl
Math 38:169–184.
29. Segura J. 2000. Bounds on differences of adjacent zeros of Bessel functions and iterative relations between
consecutive zeros. Math Comput 70:1205–1220.
30. Atkinson KE. 1978. An Introduction to Numerical
Analysis. New York: Wiley.
31. Yao GZ, Mechefske CK, Rutt BK. 2005. Acoustic
noise simulation and measurement of a gradient insert
in a 4T MRI. Appl Acoustics 66:957–973.
Concepts in Magnetic Resonance Part B (Magnetic Resonance Engineering) DOI 10.1002/cmr.b