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}]

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))]

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

omega = N[109.89*2*Pi]sigma = 56.82*(10^6)mu = 1*mu0rs = 0.00365B = 5.17

Iron sphere

omega = N[50*2*Pi]sigma = 1*(10^7)mu = 20*mu0rs = 0.05B = 1

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}]

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]

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

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

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

Print results in table form

\(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\^1\))"", 
            \*""Re[j] (MA\!\(/m\^1\))"", \ 
            \*""Im[j] (MA\!\(/m\^1\))""}}]\)