Appendix

Mathematica package sphere.nb

Download this Mathematica notebook.

Eddy currents in spherical metal samples

Plot Bessel function of first kind and first order

     Plot[{Re[BesselJ[1, x * (-1 + I)]], Im[BesselJ[1, x * (-1 + I)]]}, {x, 0, 10}]

     -Graphics-

Define phi-component of magnetic vector potential A
(T. Morisue, M. Fukumi, "3-D Eddy Current Calculation Using the Magnetic Vector Potential." IEEE Trans. Magn., vol. 24, no. 1, 1988)

     mu0 = N[4 * Pi * 10^(-7)]

               -6
     1.25664 10

     aphi[r_, rs_, theta_, omega_, mu_, sigma_] := a[rs, kx = k[omega, mu, sigma], mu] * BesselJ[1, kx * r] * Sin[theta]

     k[omega_, mu_, sigma_] := (-1 + I) * Sqrt[omega * mu * sigma/2]

     a[rs_, k_, mu_] := 3 * rs/(2 * mu0 * d[rs, k, mu])

     d[rs_, k_, mu_] := BesselJ[1, k * rs]/mu0 + (k * rs * BesselJ[0, k * rs] - BesselJ[1, k * rs])/mu

     b[r_, k_, mu_] := r^3 BesselJ[1, k * r]/(d[r, k, mu, d] * mu0) - (k * r * BesselJ[0, k * r] - BesselJ[1, k * r])/(2 * mu)

Define material parameters, frequency, external field

Copper sphere

                                                             6
     StyleBox[{omega = N[109.89 * 2 * Pi], sigma = 56.82 * 10 , mu = 1 * mu0, rs = 0.00365, B = 5.17}, ShowStringCharacters -> True]

     690.459

             7
     5.682 10

               -6
     1.25664 10

     0.00365

     5.17

Iron sphere

                                                     7
     StyleBox[{omega = N[50 * 2 * Pi], sigma = 1 * 10 , mu = 20 * mu0, rs = 0.05, B = 1}, ShowStringCharacters -> True]

Calculate eddy current density

     j[r_, rs_, theta_, omega_, mu_, sigma_, B_] := -I * B * omega * sigma * aphi[r, rs, theta, omega, mu, sigma]

     js[r_] := j[r, rs, N[Pi/2], omega, mu, sigma, B]

Plot eddy current density

     Plot[Abs[j[r, rs, N[Pi/2], omega, mu, sigma, B]], {r, 0, rs}]

     -Graphics-

Calculate magnetic moment and magnetization

     moment[rs_, omega_, mu_, sigma_, B_] := NIntegrate[r * r^2 * Pi * Abs[j[r, rs, theta, omega, mu, sigma, B]], {r, 0, rs}, {theta, 0, Pi}]

     M[rs_, omega_, mu_, sigma_, B_] := moment[rs, omega, mu, sigma, B]/(4 * Pi * rs^3/3)

     M[rs, omega, mu, sigma, B]

     604310.

     Plot[M[rs, omega, mu, sigma, B], {omega, 0, 1000}, PlotPoints -> 5, PlotDivision -> 1]

     -Graphics-

     Plot[M[rs, omega, mu, sigma, B], {omega, 0, 5000000}, PlotPoints -> 10, PlotDivision -> 1]

     -Graphics-

Simple integration formula for linear dependence of eddy current density on radius

     Abs[j[rs, rs, N[Pi/2], omega, mu, sigma, B]] * 6 * rs/20

     604611.

Print results in table form

                                                                                                                                                                                              1                1                1
     TableForm[{atab = Table[a, {a, 0, rs, N[rs/108]}], Abs[jres = Map[js, atab]]/10^6, Re[jres]/10^6, Im[jres]/10^6}, TableDirections -> {Row, Column}, TableHeadings -> {{"r", "Abs[j] (MA/m )", "Re[j] (MA/m )", "Im[j] (MA/m )"}}]

David Meeker, Finite Element Method Magnetics. User's Manual

This is the manual for David Meeker's FE package [3].

Go to David Meeker's homepage and get his free magnetics finite element package for Win95/98/NT.

Up: Calculation of eddy currents Previous: List of models Contents