MATLAB Skin effect derivation and plotting in Matlab

[RGann]

This is driving me crazy. The derivation of the current distribution in a long cylindrical wire is extremely straightforward, giving
J(r) = J(a) \frac{J_0(k r)}{J_0(k a)}
where J is the current density, a is the radius of the wire, and k is the complex wave vector, which in a metal (with nearly no permittivity) is given by
k^2 \approx -i \omega \mu \sigma
J_0 is the Bessel function of order zero. These expressions match several books I've checked. But when I try to plot the current distribution for, say, copper in Matlab, it doesn't look like the plot in my book. The code is

mu = 1.2566290e-6;
sigma = 5.96e7;
omega = 2*pi*1e4; %10 kHz
a = .05;
k = sqrt(-1i*omega*mu*sigma);
r=0:.0001:.01;
J = besselj(0,k*r)/besselj(0,k*a);
rej = real(J);
plot(r,rej)

The plot is attached. Is this correct? It doesn't seem to match the Wikipedia plot either. Does the current actually dip negative? It is otherwise qualitatively right, in that all of the current is concentrated near the edge, but I thought the max was at the very edge. What am I missing here?

Thanks

 

[jasonRF]

Two things:
1) you are only plotting out to r=0.01, instead of all the way to the outer edge of the wire (r=a)
2) you are plotting only the real part of the current. If you plot the magnitude of the current I think you will find that it is monotonically increasing out to the edge of the wire.

jason

 

[RGann]

Good catch, you have identified my conceptual problem. When I start out, I have the Maxwell equations, namely

\nabla \times \vec{H} = \vec{J} + \frac{\partial \vec{D}}{\partial t} \quad \quad \nabla \times \vec{E} = -\frac{\partial \vec{B}}{\partial t}

Then I write that \partial\vec{B}/\partial t = i \omega \vec{B}. That would imply that I need to take the real part of J at the end, because I'm using e^{i \omega t} instead of \cos \omega t. But the correct answer is to take the magnitude, that is |J|. Could you explain why this is the correct thing to do? What significance does that have?

 

[jasonRF]

The magnitude tells you how much current is flowing; the phase tells you phase shift as a function of radius. I would expect the phase to go through 2 pi phase shift for every wavelength in the wire (check out the real part of k to determine the effective wavelength in the wire). Note that the average power dissipated due to Ohmic losses is
<br /> \frac{1}{2} \Re \int \mathbf{J \cdot E^\ast}\, dV = \frac{1}{2 \sigma} \int |\mathbf{J}|^2 \, dv<br />
since E and J are in phase. So it really is the magnitude of J that you probably care most about.

jason

 

[jasonRF]

Just thought about this post again - there are cases where you DO care about the phase of J. In particular, if you want to compute the impedance (per unit length) then you need the complex current density. In the case of the wire you do find that the impedance is complex - since it has a resistance and inductance per unit length.

jason

 

[RGann]

ok, so the wrinkle is that J(r) gives information not just about current but about the phase of the current at that location, and that at a given r it is not the same at adjacent points. This is a result of dropping the time dependence out explicitly. Brilliant, thanks.