Analytical Solutions of Eddy-Current Problems in a Finite Length Cylinder

Magnetic field and eddy currents in a cylinder of finite length are calculated by separation of variables. The magnetic field outside the cylinder or inside the bore of the hollow cylinder and shell is expressed in terms of Bessel functions. Both axial and transverse applied fields are considered for the solid and hollow cylinders. The equations for the vector potential components are transformed in one-dimensional equations along the radial coordinate with the consequent integration by the method of variation of parameters. The equation for the scalar electric potential when required is also integrated analytically. Expressions for the magnetic moment and loss are derived. An alternative analytical solution in terms of scalar magnetic potential is derived for the finite length thin shells. All formulas are validated by the comparison with the solutions by finite–element and finite-difference methods.


Introduction
The paper presents an analytical method of calculation of steady-state magnetic fields and eddy-currents in the cylinder of a finite length placed in the external axial or transverse magnetic field.It is known that for the infinitely long cylinder the closed form solutions were known in different forms [1][2][3][4].However as it was pointed out by many authors [1,2,5,6] that for the finite length cylinder the general analytical solution had not been available.For axial symmetry the distribution of eddy currents induced in a conducting rod of finite length by a coaxial coil is given in [7].In general, the problem of a finite length cylinder earlier was solved by different numerical methods, and the solution was reported by several authors [6,[8][9][10][11].The conductive cylinder in transverse field was also included in so-called FELIX test problems in 80-ties [12][13] to verify the FEA codes.The analysis of eddy-currents in the cylinder is of interest for many practical applications such as electromagnetic shielding, designs of MRI /NMR components, conductive components of accelerator magnets, fusion reactors, induction heaters, turbogenerators and other electrical machines.The main aid of this paper is to derive a general analytical solution in terms of Bessel functions for eddy currents induced in a conductive cylinder by the quasi-static electromagnetic field.

Formulation of the problem
We consider the conductive non-magnetic cylinder with the constant electrical conductivity  and magnetic permeability  0 placed in the applied magnetic field oscillating with the constant frequency .The vector magnetic potential A satisfies Amperes equation where  0 is the magnetic permeability of free space , j is the current density for the area r=[0, r 0 ] , =[0, 2], z=[-z 0 ,z 0 ] and zero outside this region.The boundary condition at infinity 0 where A 0 is the vector potential of applied magnetic field.From Ohm's law we can express the current density in terms of the vector and scalar V potentials of electromagnetic field as: where  is the electrical conductivity of cylinder.The induced currents form the closed loops, and therefore current density satisfies the condition for the solenoidal field The vector potential allows us to apply an additional gauging condition so that the field can be determined uniquely from the Maxwell equations.If we apply Coulomb gauge ( ) equations ( 1), (4) become   V (5) with the coupled boundary conditions Equation ( 6) provides a zero component of eddy-currents normal to the surface Γ of the conducive cylinder.The scalar potential is defined to be zero outside the conductor.After substitution (3) in (1a) we have a system of four differential equations that are coupled.
In orthogonal system of coordinates (1a) or (1b) can be solved by the method of separation of variables when the boundary Γ coincides with the pieces of coordinates surfaces.A conductive cylinder falls in this category of geometries.However we have not found in literature the complete analytical solution for the eddy-current problem when the field is applied orthogonal to the cylinder axis.The major obstacle probably is in finding the analytical solution of the field in the whole space outside and inside the cylinder without subdividing the space in sub-regions.However in case of linear conductors there is no a rapid change of the flux density or field on its boundaries at least for moderate frequencies.Therefore the vector potential of magnetic field can be sought as a continuous smooth function in the whole space having piecewise conductive properties.The only problem seems to appear for the electric scalar potential that has to be coupled accurately with the vector magnetic potential on the boundaries of the conductor.

Cylinder in axial field
When cylinder is placed in the axial magnetic field (Fig. 1) the problem can be significantly simplified.Only one component of vector potential is enough to formulate the problem.The scalar electric potential is not needed in this case.The complex azimuthal component of vector potential A  exp(it) inside and outside the cylinder satisfies the equation and the far field boundary condition where B 0 is the amplitude of applied field .The vector potential of uniform field in form (2a) can be found by integrating the equation  [5,14].The current density in (7) can be expressed in terms of the time derivative of vector potential Since the magnetic permeability is the same through the whole space we can express the solution in the form of Fourier series accounting for the field symmetry as In accordance with Grinberg method [15] after substituting (8) in (7) and integrating over the interval [0,z  ] with the weight (2/z  )cos m z we have where Using a method of variable parameters we can resolve (9) in the sum of a general solution of the homogeneous equation and a particular solution of ( 9) where J 1 and H 1 1 are Bessel and Hankel functions, respectively, C m , D m are the constants.To keep the potential finite at r=0, we have to eliminate the functions that have singularities on the cylinder axis.This is achieved by zeroing D m = 0 in (12).To satisfy the boundary condition (2a), we have to assume C m = 0 as well.Then, after replacing the functions of complex argument in (12) by the modified Bessel functions we rewrite equations for harmonics of vector potential as: After substitution j m using (10) and approximation of the integrals by quadratures on a regular grid r n =r n-1+ h rn , n=1,N in the interval [0,r 0 ] one can derive a system of linear equations for the potential harmonics, where , h r is the grid size over the radial coordinate,  kk' is the Kronecker delta, w n are the quadrature weights.For the linear approximation of potential along the radial coordinate w n =1/2 when n=1,N and w n =1 otherwise.A better accuracy is achieved when integrating the Bessel functions around each node of the grid analytically.For example, for internal nodes where  and  K are the combinations of modified Bessel and Struve functions defined in [16].
Formulas for harmonics of vector potential ( 14) can be easily modified for the hollow cylinder with the inner radius r i .The only change required in (14) (that is now valid for r≥r i ) is the replacement of zero by r i as a lower limit in all integrals over the radial coordinate.For hollow cylinders, it is of practical interest to calculate the field on the axis.The vector potential inside the bore 0 r r The radial component of flux density outside the cylinder The field inside has a similar structure but the modified Bessel functions of the first and second kind in ( 18) are swapped.The axial component of flux density can be derived as well.In general, it will include K 1 ( m r) and its derivative over the radial coordinate outside the cylinder and I 1 ( m r) and its derivative inside the cylinder if hollow.For the thin cylindrical shells the number of steps over the radial coordinates can be significantly reduced.For just one step the set of equations ( 15) degrades to the system of equations for the potential harmonics at the same radial coordinate.
The total current induced in the cylinder is The moment of the induced currents is calculated as The power dissipated in the cylinder is Figure 6: Relative amplitude of flux density on the axis of the thin cylindrical shell of r 0 =1 cm, z 0 =1 cm, thickness h r =0.35mm at different frequencies : solid line= analytics,  = FEA..
Formulas for the field components, moment and power loss have been validated by the comparison with the results from FEA code [17].The solid and hollow copper cylinders of radius 1 cm and various lengths have been considered.The results of comparison are presented in Fig. 2-Fig.6.The number of harmonics M in these numerical tests typically is in the range between 20 and 30, and number of steps in radial grid, N is between 20 and 30 as well.The step  rn depends on the frequency and it is selected based on the skin depth.The far field boundary ( 2) is applied at z inf that typically is selected in the range of 5z 0 or higher.To reduce the oscillations in the solution, Lanczos smoothing is applied in all series over the axial coordinate.

Solid cylinder
Let's assume that the external field is applied in x-direction (Fig. 7).The induced current density in the cylinder has radial, peripheral and axial components, so the problem becomes a 3D one.We denote the components of vector magnetic potential as A r ,A  ,A z , respectively.The potential satisfies (1a) in the whole space when Coulomb gauge is selected.The current density is non-zero for the area r=[0, r 0 ] , =[0, 2], z=[-z 0 ,z 0 ], and is zero outside this region.
The boundary condition at infinity [1] Equation (1a) in the component notation can be written as The scalar electric potential is also required to satisfy (5).
Because the magnetic field has a tangential symmetry over xz-and xy-planes and normal symmetry over zy-plane the potential components can be presented as where Since the cylinder has a finite length of L =2z 0 the current density and electric scalar potential are zero beyond its volume while the magnetic field should be extended in zdirection far enough from the cylinder top and base faces.We assume that the magnetic field is equal to the applied field at z  z z  .Thus the current density and electric potential are expanded in the Fourier series over the interval z=[-z 0 ,z 0 ] while the vector potential is fit in the interval z=[z  ,z  ].For the current density and electric scalar potential in the cylinder, we also account for the model symmetry where Equation (1e) contains only one component of vector potential and can be solved similar to (7).After substitution series (24) in (1e) and integration over [-z  ,z  ] the equation for the potential harmonics is where Figure 7: Conductive cylinder of radius r 0 and axial length 2z 0 in uniform transverse field directed in x-axis Similarly to the solution of ( 9) we can resolve (29) in the sum of a general solution of the homogeneous equation and a particular solution.Because the current density harmonics (30) is dependent on the radial coordinate we use the method of variable parameters to express the harmonics A zm in the form of integrals Equations ( 1c) and (1d) for the radial and angular components of vector potential contain the crossing terms.To avoid coupling between the equations, few manipulations on (1c) and (1d) should be performed.First, let's exclude the polar angle from the consideration by substituting A r , A  into (1c) and (1d) in accordance with ( 22)-( 23), respectively r r r r r where Next, we combine (33) and (34), namely we add and subtract them from each other with obtaining the equations for the new variables Note that the variables  +r and  -r in (35) and (36) are decoupled now, and the solutions can be expressed in terms of Bessel functions of integer index.After multiplying ( 35) and (36) by 1/z  sin m z and integrating over [-z  ,z  ] the equations for the harmonics are The solution of one-dimensional equations ( 37) and (38) can be written in terms of Bessel functions similar to solution of (32) The current density harmonics j rm , j m , j zm in (32),(39), and (40) can be determined from the condition for the divergence-free current density (4).Let's impose the Coulomb gauge for the vector potential.Then the scalar electric potential satisfies the Laplace equation (5).After substituting (28) in (5), multiplying it by 1/z 0 sin k z and integrating over [-z 0 ,z 0 ] the equations for the harmonics of scalar potential are where the right part in (41) appears because the normal derivative of scalar potential on the cylinder surface at z= z 0 is different from zero.It can be further expressed through the time derivative of vector potential in accordance with The boundary conditions for the harmonics of scalar potential are Note that the boundary condition for the derivative of scalar potential (44) is a result of zero radial component of current density at r =r 0 .In general, (41) can be resolved through the integrals similarly to (14 where the coefficients G k are determined from the boundary condition (43) as: are the first derivatives of Bessel functions.The scalar potential derivative at z= z 0 can be expressed through  - 1  A div   .Indeed (4) can be rewritten as: The divergence in the right part of (47) is not zero anymore if we take the ends of cylinder into consideration.For the Coulomb gauge A div   is zero everywhere inside the cylinder but the surface where it has a jump of that is called as a surface divergence [18].Because the conductivity is discontinuous at z= z 0 and the integration totally includes the jump of  we have Thus the right part in (41) can be replaced by the integral of A div   after accounting for the angular dependence of vector potential components ( 22)-( 24) After replacing the vector potential components by harmonics over the axial coordinate and subsequent integration we obtain where

The coefficients zk
A ~are the result of fitting the axial component  z to the cosine harmonics along the cylinder length.We use values of potential up to the ends of the interval [-z 0 ,z 0 ] where -i z is zeroed because the jump in the conductivity.This allows us to account for the surface divergence at z=z 0 .The computation of q k using ( 50) is more complex than using a straight expression (42).However including A div   in the equation for the scalar potential enforces the divergence-free condition (4).To combine the advantages of two methods for computation of q k , one can add A div  (=0) to the jump of potential in expression ( 42) The latter expression for q k accurately includes the surface divergence and enforces (4) at the same time.
The harmonics of current density can be expressed through the harmonics of potentials as follows: The solutions (32),( 39),( 40) and (45) contain the unknown coefficients of harmonics A r+,m A -r , m , A zm ,V k .After approximating the integrals by a high accuracy quadratures on the grid r n =r n-1+ h rn , n=1,N built for the interval [0,r 0 ], these expressions form a system of linear equations of order of 4NM where M is the number of harmonics in series ( 22)-( 24),(28).The structure of the system is similar to (15) but it is more complex since the current density components in (32), ( 39), (40) should be expressed through the vector and scalar potential harmonics using (52).To simplify the solution, the system of equations can be solved also using the iterative approach.In this case the current density (52) is split so that the parts associated with the time derivatives of vector potential are included into the matrix coefficients while the terms associated with the electric potential are kept in the right part of system.
The analytical method has been applied to the copper cylinder of radius, r 0 of 1cm and half length, z 0 of 1cm.The number of harmonics M has been selected in the range ≥30, and the far field boundary conditions have been applied at z≥3z 0 .The distribution of current density is presented in Fig. 8-10 in comparison with the results from FEA [17] and FDM [19].A pure hexahedral mesh in all directions has been used in both numerical methods.The has been refined in both methods they provided the close solutions in terms of local fields, energy and moment of induced currents.For FDM, the symmetry over the polar angle ( 22)-( 24) has been used first, and then four coupled 2-D equations have been solved.2-D equations have been approximated using the integral-interpolation method [19].The finite-differential schemes for the equations are similar to those reported in [20].

Hollow cylinder
Similarly to the cylinder in the axial field the calculation formulas for the transverse field can be applied in the case of the hollow cylinder of the inner radius r i , outer radius r 0 and axial length 2z 0 .The equations for the components of vector potential (32),(39),and (40) differ only by the integration range.The lower zero limits in integrals (32),( 39) and (40) should be replaced by r i .Thus the grid along the radial coordinate is built in the interval [r i ,r 0 ].The expression for the scalar electric potential contains an additional Bessel function of second kind K 1 ( k r) with the constant coefficients This is because the scalar potential is considered for r>0.Thus the solution (53) that includes infinite functions K 1 ( k r) at r=0 is finite everywhere in the hollow cylinder.The boundary condition on the outer surface of the hollow cylinder is the same as (44).Similarly, on the inner surface we have The coefficients G k and F k are determined from the boundary condition system ( 44) and (54) as For Combining ( 56)-( 58) we obtain the flux density on the axis of the cylinder Other components are zero on the axis of the cylinder, thus B x gives us the magnitude of the field on z-axis.
Similarly to the solid cylinder the formulas for the hollow cylinder have been verified by the comparison with the induced currents computed by numerical methods (Fig. [11][12][13][14] in the range of frequencies.The steps in the radial grid have been selected based on the frequency of the applied field.For example, at least 30 steps (N=31) over the radial coordinate is used for the frequency of 2 kHz in the case study with the results presented in Fig. 14.
Figure 11: Radial component of current density in YZsection of a hollow cylinder of r 0 =1 cm, r i =0.5 cm z 0 =1 at f=300 Hz, B 0 =10 mT: solid line = real j z , analytics, dash line =imaginary j z , analytics;  = real j z , FEA,  = imaginary j z , FEA.
Figure 12: Peripheral component of current density in XZsection of a hollow cylinder of r 0 =1 cm, r i =0.5 cm, z 0 =1 cm at f=300 Hz, B 0 =10 mT, data labels are the same as in Fig. 11.
Figure 13: Axial component of current density in YZ-section of a hollow cylinder of r 0 =1 cm, r i =0.5 cm, z 0 =1 at f=300 Hz, B 0 =10 mT at f=300 Hz, B 0 =10 mT, data labels are the same as in Fig. 11.
Figure 14: Field on axis of the hollow cylinder of r 0 =1 cm, r i =0.5 cm z 0 =1 cm at different frequencies: solid line = analytics,  =FEA.

Cylindrical shell
A cylindrical shell is a special case of the hollow cylinder when we can neglect a radial dependence of current density [2].First, consider how the thin shell approximation can simplify the solution in case of hollow cylinder.Since the thickness of the shell typically is significantly less than its radius it is practical to consider a one-layer (N=1) approximation for the derived solution.
Indeed if we neglect the dependence of electric scalar on the radial coordinate, from (41) we can derive This simple approach has few limits.A one-layer approximation totally ignores the difference in the current density across the shell thickness although it accounts for the radial component of vector potential, and hence for the radial current density associated with its time derivative.Also the one-layer model may lead to the ill-defined matrix for the harmonics of vector potential starting at high frequencies.It is illustrated in Fig. 15 for a 0.35 mm thick copper shell of 1 cm radius.When the frequencies are low the field on the axis is in a good agreement with the accurate approximation of the shell by five layers (N=6).When the frequency increases the matrix condition number becomes too large, and the approximation (60) at N=1 starts to give us an inaccurate field on the axis.It should be noted that the calculation of q k in (60) using a straight expression (41) gives a more accurate result since in case of just one layer the derivative dA r /dr in the expression for A div  is hard to determine accurately.
The alternative solution for the thin shell can be derived using a concept of stream function [2] and a scalar magnetic potential [21][22] .The magnetic scalar potential satisfies the Laplace equation everywhere but the shell.The whole space along the radial coordinate is divided by the shell into two regions: r > r 0 (marked by upper index +) and r < r 0 (marked by upper index -).The scalar potential inside (-) and outside (+) the shell can be written in terms of Bessel functions as where C m -and C m + are the constants.Because the normal component of the field is continuous everywhere including the shell surface the constants in (61) relate as: The field tangential component has a jump by the amount of surface current density when crossing the shell thickness from inside to outside The scalar magnetic potential has also a jump [18] on the shell surface u=U + -U -.After expressing the field through the scalar magnetic potential in (63) one can write where index s means that the nabla operator  s works on the surface , n  is the unit vector normal to the surface.
In terms of components (64) can be rewritten as In accordance with [2] the stream function(,z) relates to the current density as or in terms of components The stream function satisfies the equation where H r is the radial , i.e. component of the field normal to the cylindrical surface.Comparing (65) with (67) we can conclude that the potential jump u satisfies the equation similar to (68), namely The second equation for the coefficients C m -and u m can be determined from the potential junction at r=r After integration (72) over [0,z 0 ] and substitution C m + from (62) we obtain Equations ( 71),(73) form the system of linear equations for the coefficients C m -and u m that can be solved directly.However (71),(73) can be solved by iterations with some under relaxation factor.We assume C m -=0 for the first iteration, find u m from (71) and update C m -using (73) with under relaxation.We iterate (71) and (73) until the convergence tolerance is achieved.The iterative method implements the perturbation of the normal field H r by the induced field starting with H 0 as a first approximation for the field on the shell surface.The iterative algorithm gives a more steady solution at high frequencies when the matrix of system (71),(73) becomes ill defined .
The magnetic field on z-axis is calculated as (74) The thin sheet approximation using the scalar potential works well even at high frequencies (Fig. 15).The discrepancy from the accurate solution (N=6) increases at frequencies above 15 kHz when the shielding capacity of the shell is close to saturation.The solution (71),(73) becomes unresponsive to the frequency because of difficulties in the accurate estimation of the normal component of the resultant field H r on the shell surface (right part of (69)).

Magnetic moment and power loss
The magnetic moment M x of currents induced in the cylinder is The results of calculation of power loss and magnetic moment of the copper solid cylinder using (75a)-(76) in test models agree well with the data from FEA and FDM (Fig. 16).

Conclusions
The analytical solutions for the eddy current problem in the conductive solid and hollow cylinders in cases of axial and transverse fields have been derived.The magnetic field outside the cylinder or inside the bore is expressed in terms of Bessel functions.Formulas have been verified by the comparison of the calculation results in test models with the FDM method and modern FEA codes.Formulas for the power loss and magnetic moment of currents induced in the cylinder have been also found.These have been validated over a wide range of frequencies within the limits of quasistatic approximation.

Figure 1 :
Figure 1: Conductive cylinder of radius r 0 and length 2z 0 in uniform axial field B 0 directed in z-axis.where ) 2 /( ) 1 2 (    z m m   , z  is the axial coordinate where the field satisfies the far-field boundary condition (2a).
outside the cylinder the harmonics of potential are calculated in terms of harmonics on the grid within the cylinder as:

Figure 9 :
Figure 9: Peripheral component of current density in XZplane across the cylinder of r 0 =1 cm, z 0 =1 cm at f=300 Hz, B 0 =10 mT, data labels are the same as in Fig. 8.

Figure 10 :
Figure 10: Axial component of current density in YZ-plane across the cylinder of r 0 =1 cm, z 0 =1 cm at f=300 Hz, B 0 =10 mT, data labels are the same as in Fig. 8.
. The magnetic field is split into the applied field, 0 H  and the field from the eddy currents, e H  where the normal component of x-oriented field on the shell surface is H 0r =H 0 cos.Outside the cylindrical shell we can use the scalar magnetic potential U to calculate the field from the eddy currents U H e    potential jump u exists only at r=r 0 and 0  z  z 0 and it can be presented it in the form of series 70) in (69) and integrating over [0,z 0 ] the coefficients u m are determined as follows: