Design and Characterization of Self-Shielded MRI Gradient Coils With Finite Track Width

Gradient coils are essential for the performance of magnetic resonance imaging systems. Usually, coils are designed assuming thin wire tracks. Here, we design an MR gradient coil set using a more general approach considering the exact track width using the discrete wire approach. The effect of track width on the DC current density distribution and resultant magnetic fields using both loop and Golay coils are first demonstrated. Both, self-shielded X and Z gradient coils of definite width/thickness are designed and optimized. The resistance and inductance of the coils are calculated using the stream functions approach. Track current distribution was used to compute the magnetic fields over the desired volume, and at the cryostat. The linearity of the magnetic field over the volume, the figure of power, and the shielding ratio of the coil are used as parameters in the optimization process. The DC characteristics of the designed coils with definite (small) track width and thickness were compared for verification to that of the corresponding thin wire design where they were found to have approximately similar characteristics. Using our design methodology, the coils’ frequency-dependent resistances and inductances were directly/efficiently calculated. The harmonic and transient eddy current interactions between the longitudinal and transverse gradient coils were computed where track slitting was employed to reduce such interactions. This work stresses the importance of considering coil track width in the design process particularly for wide tracks as well as computing the coil’s figure of merit, harmonic and transient coil characteristics/interactions.

[1], [2], Fig.1. 23 They are commonly designed by one of two methods [3], 24 [4]: the discrete wire method, and the continuous current 25 distribution method. 26 The associate editor coordinating the review of this manuscript and approving it for publication was Jinhua Sheng .
In the discrete wire method [5], [6], [7], [8], [9], the 27 gradient coil turns of specific geometry (e.g., circle, ellipse, 28 semi-elliptic, etc.) and different sizes are positioned on the 29 surface of two cylinders (one for primary and the other for 30 shielding coil). The size and the position of the turns are 31 iteratively adjusted through a computerized algorithm until 32 specific target fields on the Diameter of Spherical Volume 33 (DSV) and the cryostat are obtained. 34 In the continuous current distribution method [10], [11], 35 [13], [14], [15], [16], the target fields on the DSV and the 36 cryostat are predefined where the desired target field (B z ) 37 components on DSV satisfy the gradient linearity constraint 38 where B is the magnetic field, − → R is the vector from the coil 60 segment to the target field point, and R is the length of the 61 vector ( − → R ), − → dl is the element vector along the coil path C,

62
I is the line current that passes in the coil, and µ 0 is the 63 permeability of free space which equals 4π × 10 −7 H m −1 .

64
Commonly, the gradient coil is assumed to have a thin-wire 65 conductor. In reality, coils are not composed of thin wires but 66 they have definite widths as discussed in [20] and [21] and 67 evident from images of gradient coil structures used by some 68 MRI scanner manufacturers whose exact parameters are not 69 necessarily published. In this paper, we investigate the impor-70 tance of considering the coil's track width. We exploited the 71 advantage of the discrete wire method to design gradient coils 72 considering the track width and thickness of the coil's turns. 73 The tracks of the coil are meshed into structured triangular 74 elements and the current density on the coil's tracks is repre-75 sented in terms of the stream functions [22], [23], [24], [25], 76 [26], [27]. Calculating the stream functions is followed by stat cylinder are calculated [28]. Also, cross coil eddy current 83 interaction is calculated where it is shown that track slitting 84 reduces such effects as recommended in [22]. We stress in 85 this work that the X, Y, and Z gradient coil sets affect each 86 other and this is to be considered in coil design in general.

88
The self-shielded gradient coil is represented by two sepa-89 rated concentric cylinders, Fig.1 (a). For the Z gradient coil,  each cylinder is divided transversely into two identical halves, 91 Fig.1 (b). Four groups of turns are arranged on the cylinders: 92 two groups on the inner cylinder for the primary coil and the 93 other two in the outer cylinder for the secondary (shielding) 94 coil. However, For the transverse (X gradient) coil, each 95 cylinder is divided into four identical quadrants, Fig.1 (c). 96 Eight groups of turns, of specific track width and thickness, 97 are arranged on the two cylinders: four groups on the inner 98 cylinder and the other four on the outer cylinder.

99
Each turn's track is meshed into a single layer of structured 100 triangular elements, Fig.2. The nodes in a triangle are locally 101 numbered from 1 to 3 and globally numbered by unique 102 numbers. The direction of local numbering must be consistent 103 (clockwise or counterclockwise) for all triangles in all turns' 104 tracks. According to the continuity equation and using the 105 stream functions [24], [26], [27], [29], the surface current 106 density − → J s on the coil is represented by the stream functions 107 as: (3) 110 VOLUME 10, 2022 [29] and it can be given as: where is a vector that contains the stream functions of  The impedance matrix Z can be given as: where ω is the angular frequency and j is the imaginary unit.

135
If the nodes n and m are shared among the group of 136 triangles N and M , respectively, then the resistance element 137 R nm and inductance element M nm can be given as: The vector r N is pointing to the centroid of the triangle, with 156 an area A N , in the group of triangles N . The vector r M p is 157 pointing to the point p on the middle edge of the triangle, 158 with an area A M in group M . The weight w p is associated 159 with point p which in our work is considered to be equal to 160 1/3. In the case that the nodes n and m belong to the same 161 triangle (N = M ), a closed form of the double integral is 162 used as described in [26].

163
The total current I that passes in the turn track is equal 164 to the difference between the values of stream functions of 165 the boundaries of the turn. As shown in Fig.2, the stream 166 functions of the first boundary nodes depicted by the red color 167 are set to zeroes (for simplicity) while the stream functions of 168 the boundary nodes with green are set to the value of current 169 I . To reverse the current direction, the setting of stream func-170 tions for boundary nodes is simply exchanged. The turn track 171 has two categories of stream functions: the stream function 172 vector b for the boundary nodes (they are known and forced 173 to have certain values as discussed in the above setting) and 174 the stream function vector for the internal nodes i (they 175 are unknown and need to be solved). By breaking down the 176 stream functions vector and impedance matrix Z into a 177 combination of the internal and boundary nodes, the circuit 178 equation (6) can be rewritten as: where Z ii , Z bb , and Z ib are the impedance sub-matrices asso-181 ciated with the internal/internal, boundary/boundary, and the 182 internal/boundary nodes of the turn's track, respectively. The 183 solution to the internal stream functions is given as: The current density can be computed from the stream func-186 tion values using equation (4).

187
For a turn, the calculated i and b are concatenated 188 to create the vector . The per turn resistance R t and 189 inductance L t are calculated using the following formulas 190 as directly inferred from the definition of the electrical and 191 magnetic energies in [27]: The above equations are valid for solid turn track, Fig.3 195 (a). However, it is beneficial to study also the turn with a 196 slitted track, Fig.3 (b). Slitting the turn track into sub-tracks 197 plays important role in reducing the eddy current induced on 198 the turn. In the slitted turn track, Fig.3 (b), the assignment of 199 upper and lower boundary nodes b (green and red) is applied 200 as previously discussed with solid turn. To prevent the current 201 from crossing the boundary of slits, a boundary condition 202 is set where the stream functions of the slits' nodes (cyan 203 and yellow) should be equal to an unknown value that needs 204 to be determined. The stream functions of the slits' nodes 205 matrix H [27] which satisfies the following: The matrix H is a binary matrix with a dimension of N nodes × 211 N ind . N nodes is the number of all nodes on the turn, and N ind 212 is the number of independent nodes.

213
The solution to the stream functions is given as: The vector i is calculated from equations (15,16). Similar to 216 the solid track, i and b for a turn are concatenated to create 217 the vector . And per turn resistance and the inductance are 218 calculated as in equations (13,14).

219
A computational framework from the above equations is  The current density on the turn track, the resistance, and the 245 inductance of the turn at different frequencies are calculated 246 via the framework and also compared with Ansys results of 247 the same turn configuration. For more validation, DC resis-248 tance and inductance of the turn are compared with closed-249 form in [31] and [32].

250
The resultant gradient field for a Golay coil [33] of differ-251 ent track widths (1 to 7 cm) is compared. The radius of the 252 Golay coil is set to 10 cm while the current is 100 A. The 253 Golay's track is meshed similarly to the single loop. For validation, self-shielded Z gradient (longitudinal) and X 256 gradient (transverse) coils of specific track width are designed 257 in this work (comparable to previously published coil designs 258 [5], [28]). The radius and the length of primary and shielding 259 coils of the Z gradient differ from those of the X gradient 260 as will be discussed in the next subsections. In both gradient 261 coils, the DSV is 50 cm. The radius and the length of the 262 cryostat cylinder are 43 cm, and 146 cm, respectively. For all the gradient coil turns, the current densities in the 265 triangular elements on all turns are calculated. The magnetic 266 field from these currents on the DSV (B z component) and the 267 cryostat cylinder (B x , B y , and B z components) are calculated 268 using Biot-Savart law [17], [18], [19] in terms of volume 269 currents as follows: where v is the track elements volume, and − → r is the vector 272 from the triangular elements on the coil tracks to the target 273 points.

274
The following parameters are used as design metrics for 275 coil optimization and performance evaluation [5], [6], [28]: 276 1) The average gradient strength (G m ) over the DSV: where G i is the gradient strength and M is the overall 279 number of points on the DSV. 2) The coil efficiency (η) : where I is the current that passes in the gradient coil. 283 3) The linearity error (LinE): where R is the resistance of the coil.

290
5) The average shielding ratio (SHRa) at the cryostat: extend from 0 to ±z l/2 (z l/2 is the half-length of the cylinder).

301
The coordinates of the turns of the primary or shielding coil 302 as a function of a radius can be expressed as: and N is the number of turns on one cylinder side. A constraint 308 on the distance between any two consecutive turns is given as where its purpose is to minimize the multi-objective function: 324 where x is a vector that represents the locations of the turns.

325
The weighting factors α 1 , α 2 , and α 3 were set, in this work, 326 to 1/3 where equal priority is given to each optimization 327 term. To avoid scaling issues, the objective parameters were 328 normalized according to [34] where The target of the parameters LinE, SHRa, and Gm were 331 set to 5%, 95%, and 45 mT/m, respectively. In both the 332 primary and shielding coils, the locations of the turns are 333 symmetric in z-direction in both halves of the cylinder. The 334 optimization of the turns' location is done only for onehalf 335 of the cylinder. The upper and lower bounds of the turns' 336 location were set between 0-1246/2 mm and 0-1286/2 mm 337 for the primary coil and the shielding coil respectively. Linear 338 inequality constraints were set to the distance between the 339 consecutive turns ( Ax ≤ b in fmincon) 340 After achieving the final Z gradient coil, the whole coil 341 performance is calculated for the solid, slitted tracks as well 342 as for the thin wire coil.

343
In the slitted track Z gradient coil configuration, the track 344 of 10 mm is slitted into three sub-tracks by two slits of 345 1mm width. Each sub-track is divided into three sections of 346 structured triangular elements The turns are presumed to have a quasi-elliptical shape which 349 is mathematically represented as in [5] and [6]. Similar to 350 the approach followed in [28], the turn track width is set to 351 around 5 mm while its thickness is set to 2 mm. Because of the 352 curvature of the coil turn, the turn track width varies between 353 5-5.5 mm. Similar to the Z gradient coil, the X gradient turns' 354 tracks are meshed into a singular layer of structured triangular 355 elements. The coil turn in the circumferential direction is 356 discretized into segments of 20 mm while in the track's 357 width direction is meshed into five sections of structured 358 triangular elements. The radius and the length for the primary 359 cylinder are 320 mm and 1286 mm respectively, while they 360 are 370 mm and 1326 mm for the shielding cylinder similar 361 to [5].

362
Both the primary and the shielding cylinders have four 363 symmetric quadrants. The number of the primary coil turns is 364 searched over 18-23 turns for one quadrant. All the turns have 365 the same center. In the searching process, the primary turns 366 are moved so they occupy 24 possible locations and their 367 center can be located at 50 available locations. The magnetic 368 field over the DSV from the all-possible combinations of 369 primary turns as well as their total resistance are calculated. 370 The elected number for the primary coil has the best FoP with 371 LinE < 10%. Only the combinations containing the elected 372 number of turns with linE<10% are saved as candidate pri-373 mary locations for the next process which involves shielding 374 turns to elect the optimal whole gradient coil.

375
For the shielding coil, the number of turns is searched 376 over 10-16 turns (for one quadrant) which occupy 22 pos-377 sible locations and their center can move over 100 locations. 378 Using brute-force search, the whole gradient coil is searched 379 over the candidate's primary combinations with shielding 380 combinations to get the final whole gradient coil which has 381 the best FoP and is constrained by LinE < 9% and SHRa 382 > 85%. Fig.4 shows the possible locations which can be 383 occupied by shielding turns and their center can move up 384 or down (indicated by the red arrow) over 100 locations. 385 Fig.5 illustrates a flow chart of the algorithm used to optimize 386 the whole self-shielded X gradient coil by brute-force search. 387  Considering the coil with tracks makes it possible to study 420 the interaction between the gradient coils. Via the imple-421 mented framework, the harmonics and transient eddy currents 422 interaction are studied between the designed X and Z gradient 423 coil. In both harmonic and transient interaction analysis, the 424 X gradient coil is activated while the Z gradient coil is set 425 non-active (passive). Two configurations of Z gradient coil 426 are included in the analysis: the solid tracks, and the slitted 427 VOLUME 10, 2022 where dv is the volume element, σ is the conductivity of 442 the coils' material, and * denotes the complex conjugate.

443
The power dissipation (in dBm) due to the eddy currents is where τ is the ramp-up and ramp-down times as depicted in The following differential equation gives the temporal solu-466 tion to the internal nodes stream functions i (t): where M ii , R ii are the inductance and resistance matrices 472 associated with the internal/internal nodes. M ib , R ib are 473 the inductance and the resistance matrices associated with 474 the internal and the boundary nodes, and M bb , R bb are 475 the inductance and resistance matrices associated with the 476 boundary/boundary nodes. By diagonalizing the matrix W , 477 we can write: where U is a matrix contains the eigenvectors of the matrix 481 W, D is a diagonal matrix whose diagonal elements are the 482 eigenvalues (λ 1 , . . . , λ k ) of the matrix W, and k is the number 483 of internal nodes. The general temporal solution to the above 484 first-order differential equation (29) is given as: where t 0 is the initial time, λ is the eigenvalues vector 489 [λ 1 . . . λ k ], and ζ is a dummy variable.

490
The calculated i and b are concatenated to create the 491 vector . The current densities in the triangular elements 492 on the tracks of the active coil (X gradient) and passive 493 coil (Z gradient) are calculated as in equation (4). The net 494 magnetic field from the two coils is calculated on a target 495 point inside the coils. An appropriate boundary condition is 496 set to prevent the excitation current and the induced eddy 497 currents from crossing the edge of the tracks where equation 498 (15) is involved in the above formulations. We preferred to 499 performed once for the system matrices. Since the shape of 506 the pulse is predefined, this allows a direct numerical solution 507 for the transient problem as previously shown in [31]. The computed current densities Am −2 from the framework 511 and Ansys is plotted for the single turn, Fig.7 (a). The Ansys 512 result, Fig.7 (b), partially verifies the core of our computa-513 tional framework where it is obvious that the two current 514 distributions are comparable. The current distribution on the 515 turn is not uniform and it depends on the geometry of the turn.

516
The inner side of the turn track has more current density than 517 the outer side.   are plotted while B x is zero. It is obvious from the plots that 555 when the track width is very small (close to the thin wire), 556 the resultant magnetic fields are identical to those induced 557 by the thin wire. When the track width of the coil increases, 558 a significant difference in the magnetic fields is noticed.

559
Similarly, the current density distribution on the Golay's 560 coil track is computed as shown in Fig.7 (f). The magnetic 561 field B z of the Golay's coil for several track widths as well as 562 the thin wire are displayed on Fig.7 (g). The magnetic fields 563 B x and B y are not included because they are equal to zero. 564 Again, the resultant gradient fields for Golay's coil of wide 565 tracks differ from the thin wire assumption. The optimized Z and X gradient coils are illustrated as 3D 568 plots in Fig.8 (a) and Fig.8 (b), respectively. They are demon-569 strated together in the same 3D plot in Fig.8 (c). The trans-570 verse view and the coronal section view are also illustrated in 571 Fig.8 (d) and Fig.8 (e), respectively. The colors on the plots 572 indicate the direction of the currents on the turns.

573
The DC performances of the Z and X gradient coils are 574 tabulated in TABLE 2 and TABLE 3, respectively. The perfor-575 mances of three different configurations of the Z gradient coil 576  It is also worth noting that the linearity error for 601 DSVs ≤ 50 cm is < 5% as preferable for effective MRI. The 602 linearity of the gradient field (B z ) of the designed gradient 603 coils are illustrated on different planes inside the DSV as 604 shown in Fig.9. Frequency has a noticeable effect on the resistance and a 608 slight effect on the inductance of the coil as shown in TABLE 609 4. As the frequency increases, the resistance of the coil 610 increases noticeably due to the reduction of the effective 611 cross-section of the coil.  The reasons for the slight decrement of inductance with 613 increasing the frequency are discussed in detail in [35].   The temporal magnetic field B z is produced by the active X gradient and passive Z gradient at an arbitrary point P(315,0,300) mm inside the coils with three configurations: Z gradient of wide track, Z gradient of slitted track, and no eddy currents. electromagnetic interaction between the X and Z gradient coil 626 is reduced by the presence of slits in the Z gradient tracks as 627 suggested by [22]. 628 In the transient eddy currents analysis, the magnetic field 629 B z generated by the active X gradient coil and the passive 630 Z gradient coil is calculated at an arbitrary point inside the 631 coils. The magnetic fields B z are plotted versus the time 632 corresponding to Z coil configurations and no eddy currents 633 case as shown in Fig10. The transient eddy current effect is 634 noticeable in the case of the solid track while it is almost 635 not notable in the case of the slitted tracks and the thin wire. 636 The eddy current induced on the solid Z gradient produces 637 a secondary magnetic field that has a distortion effect on the 638 magnetic field created by the X gradient specifically in the 639 transient current locations.