The Skin effect. Impedance of a straight wire at radio frequencies


  1. Basic theory
  2. Octave script
  3. Bessel's in c-code
  4. Ad hoc formula
  5. Coaxial tube
  6. Hankel: A solution for the outer tube
  7. 9. Bessel-bignum

1. Basic theory, SI units

From Maxwell's equations one can get the wave equation for the electrical E-field in a conductive medium:

Form_1 Maxwell

Deduction

The E-field is a vector, but in the discussion below only the component parallel to the axis is of interest.

The usual technique to solve such an equation is to extract a fixed frequency time function and write the equation as


Form_2

For a metal wire the second term in this equation can be dropped, unless we are interested in the optical range.

For a circular wire it is convenient to use cylindrical coordinates:


Form_3.png cylindrical coords

In cylindrical coordinates, if we assume a rotational symmetry around the z-axis, and z-independence, the equation can be written


Form_4.png J0

This is Bessel's differential equation with n=0, and solutions are the Bessel functions J0 and Y0. The one of interest here is


Form_5.png E(R)

The other solution Y0, has a singularity in the center and is not of interest for a solid wire, but useful when we need the solution for e.g. the outer conductor of a coaxial tube.

More about Bessel functions.

With this solution as a base, the field can be written as proportional to it:


Form_6.png A_definition

The expression gives the electric field inside the wire, and the current density is proportional to the field according to Ohms law.

For simplicity, let A denote the root-expression:


./Form_7.png ac_currrent

Note that A is complex and E, J0, J1, g and z are complex entities.


./Form_7.png ac_currrent


To get he total current in the wire integrate the current density from radius 0 to radius R. R denotes the boundary radius, i.e. half the wire diameter.


Form_9.png E

The interesting solution is the one that fulfills the boundary condition at the wire surface r = R. When one volt is connected over a wire of one meter length, the field at the boundary is one V/m.


Form_9.png zeta

Assuming this boundary field we get the current i_ac in the wire as the first formula below, and the wires impedance z is simply 1/g, as in the second formula below. (There is a recurrence formula for the Bessel functions implying that the integral J0(Ar)*r*dr is simply r/A*J1(r). [Standard Mathematical Tables]


Form_11.png zeta


Form_12.png zeta

Also check a nice animation (for a plane wave) at <http://fermi.la.asu.edu/w9cf/skin/index.html >

Markus Meisinger has pointed out to me that the wire conductance 1/z, can also be represented as a sum over roots to the Bessel functions.

Please read more about roots here.


2. Octave script

The following octave script calculates the resistance and the inner inductance per unit length of a straight metal wire.
Here you can find an Octave script Bessel that does the job.
For a wire of thickness 1mm it would produce following output, in case the script doesn't work on your computer.



freq,  delta,   L=imag(z)/omega,   real(z)/R_DC,  tube model,  Col4/Col5
Hz      [m]          [H]              [ratio]      [ratio]      [ratio] 
1e+00   0.066006   -5.000000e-08     1.000000      1.000000    1.000000
1e+01   0.020873   -5.000000e-08     1.000000      1.000000    1.000000
1e+02   0.006601   -4.999998e-08     1.000001      1.000000    1.000001
1e+03   0.002087   -4.999829e-08     1.000069      1.000000    1.000069
1e+04   0.000660   -4.982952e-08     1.006822      1.000000    1.006822
1e+05   0.000209   -3.918752e-08     1.451263      1.513669    0.958772
1e+06   0.000066   -1.315192e-08     4.049728      4.055194    0.998652
1e+07   0.000021   -4.173172e-09    12.231119     12.232539    0.999884
1e+08   0.000007   -1.320079e-09    38.126500     38.126924    0.999989

The following plots are produced from this script are showing the wire impedance, R/R_dc and L, for various frequencies.
These entities are plotted versus beta*r, where

beta = sqrt(omega*sigma*u0/2) = 1/skin_depth ;



Bessel,skin depth

3. Bessel's in c-code

The Bessel functions are part of libm (#include <math.h>), and are there named j0, j1, jn.
These functions only take real arguments, but we need to use complex arguments to calculate wire impedance.

Fortunately they can be expressed as a series that can be expanded with complex arguments.
./jn_series.png
Here is a c-routine: jn_series.c for this.

The series expansion method is problematic for arguments > 20 because the terms grow.

For arguments > 20 it is better to use the formula
Form_13.png jn_expression

and do the integration numerically.
Here is a c-routine: jn_exp_int.c for this.
Unfortunately this routine is not very fast, but it will only take a few milliseconds to run on a computer.
Find also a program including functions that computes a wires impedance per meter.

4. Ad hoc formula

In order to get a useful and fast c-function to calculate wire impedance, I have written a reasonably fast ad hoc formula that works over the whole range. For beta*R below 16 the series expansion method is used in this formula. For higher values of beta*R the results from the Bessel's integral ( N=9999), are used to adapt a polynomial of degree 8 to two or three ranges of beta*R. Well, actually it is a polynomial of 1/R as function of 1/(beta*R), that way it works a lot better.

The formula is accurate to eleven digits over the whole range, although the the accuracy remains to be proved for beta*R-values over 500. There is no reason to believe that something strange should happen there, however.

It would be possible to make the formula a little faster by adapting polynomials to lower ranges of beta*R, but the number of coefficients grow. The formula above is a reasonable compromise.

Those who do not believe in these calculation can find a raw integration of the wave equation. You can compare the result, and for high values of integration steps (N) the results should approach each other. Check it out.

5. Coaxial tube

The inductance per unit length of a coaxial tube, if the outer tube is superconducting, can be written
L_coax = L_internal + L_external.
The term L_internal refers to the inductance associated with the magnetic field lines located inside the wire, penetrating the metal.
L_internal is what is shown in the plot above. It is what the Octave script says.
L_external is easy to compute: (u0/(2pi)) ln(r2/r1). It can be found in any book about electronics. Beware of r2. If the outer tube is not superconducting, there may me a skin effect also in the outer tube.
L_internal is sometimes given as u0/8pi. That is an exact value in the low-frequency limit, where the current distribution in the wire is uniform.

It is amazing to observe that the internal inductance is not only exactly u0/8pi, it is also exactly 50 nH/m. The only deviation from exactness is the fact that the metal might have some relative permeability, if so we use L_internal = u0.ur/8pi, and the fact we do not live long enough to reach zero frequency.

6. Hankel: A solution for the outer tube

If the outer tube is not superconducting, we must search for a solution for Bessel's differential equation also there.

Assume we have a circular symmetric coax tube with radii R1, R2 and R3.
R1 is the radius of the inner conductor, R2 is the inner radius of the tube and R3 is the outer.

Since there is no restriction regarding the center, the Bessel function J0 is not the only solution to the equation. The general solution is a linear combination of J0 and Y0.

As the Y0 reach high values for big abs(arg), it would seem more natural to choose the Hankel function as a base. The Hankel H0 of first kind is a certain combination of J0 and Y0 that reach low values for big abs(arg). It represents an outgoing wave.

So our ansatz will be:

  1. Ez= c0*J0 + c1*H0
  2. Find the values of c0 and c1 that meet the proper conditions.
As just said, the natural choice would be a Hankel function of first kind, since its values behave as we would expect when the radius grows.

The boundary conditions are not directly available, but let's assume that the coax tube is short compared to the wavelength and that the tube is not grounded at both ends. Under this assumption the current flowing into the central conductor will return through the outer tube. It is not strictly true, because part of the current can flow as a displacement current in the outer space, but if the tube i short, this fraction can be neglected. Let us arbitrarily set the current through the tube to one ampere.

So the following two conditions are to be met:

  1. The current in the tube is one ampere
  2. The magnetic B-field at R3 is zero
The second conditions may need some comments. The reason for this is that the current in the inner conductor will exactly cancel the current in the tube. The implication is that the B-field is zero also at the surface R3, and even in the metal infinitely close to R3. So the Maxwell's equation curl(E) = - dB/dt implies that curl(E) is zero at the surface R3.


./Form_14.png maxwell3

Now, the theta component of curl(E) is dEr/dz - dEz/dr. The first term is zero at the boundary R3, so the condition gives that dEz/dr is also zero at the surface R3. Accepting this, we can solve the matrix and obtain c0 and c1.

Fortunately, as seen in section 1, it is easy to take the derivative of the Bessel functions, and also to integrate them.

Form_15.png hankel

It is not so easy to compute the Bessel and Hankel functions, but they are available in the Octave program. You can find a calculation of the coax impedance in the Octave-script Coax. The program is also available in c-code in the file hankel.c.gz, which also includes c-code algorithms to compute Bessel and Hankel functions.
Compile it with the command 'gcc -lm hankel.c'.

Here is an example of what that program can produce:


 The impedance of a coax. Loss and delay
 Coax dimensions in mm: R1=5.00e-04 R2=1.00e-03 R3=1.50e-03
 DC-values:
 R=2.627966e-02   C=8.025903e-11
 AC-values:
 f[Hz]    L[H/m]      R[ohm/m]    Z0[ohm]     loss[dB/m]  delay[s/m]
 1.0e+00  2.2137e-07  2.6280e-02  5.2519e+01  4.3463e-03  4.22e-09
 1.0e+01  2.2137e-07  2.6280e-02  5.2519e+01  4.3463e-03  4.22e-09
 1.0e+02  2.2137e-07  2.6280e-02  5.2519e+01  4.3463e-03  4.22e-09
 1.0e+03  2.2137e-07  2.6283e-02  5.2518e+01  4.3469e-03  4.22e-09
 1.0e+04  2.2094e-07  2.6582e-02  5.2467e+01  4.4005e-03  4.21e-09
 1.0e+05  1.9891e-07  4.3432e-02  4.9782e+01  7.5779e-03  4.00e-09
 1.0e+06  1.5838e-07  1.2883e-01  4.4422e+01  2.5189e-02  3.57e-09
 1.0e+07  1.4489e-07  3.9765e-01  4.2489e+01  8.1291e-02  3.41e-09
 1.0e+08  1.4061e-07  1.2483e+00  4.1856e+01  2.5905e-01  3.36e-09
 1.0e+09  1.3926e-07  3.9386e+00  4.1654e+01  8.2129e-01  3.34e-09
 1.0e+10  1.3883e-07  1.2446e+01  4.1590e+01  2.5993e+00  3.34e-09
 1.0e+11  1.3869e-07  3.9349e+01  4.1570e+01  8.2218e+00  3.34e-09



In the table above, the characteristic impedance Z0 is set to Z0 = sqrtl (L /C).
This may be wrong at low frequencies. If we set Z0 = sqrtl ((R + omega*L )/(G + omega*C)),
and assume G is infinity, we get the following table:


Hankel:
The impedance of a coax. Loss and delay
Coax dimensions in mm: R1=5.00e-04 R2=1.00e-03 R3=1.50e-03
DC-values:
R=2.627966e-02   C=8.025903e-11
AC-values:
f[Hz]    L[H/m]      R[ohm/m]    Z0[ohm]     loss[dB/m]  delay[s/m]
1.0e+00  2.2137e-07  2.6280e-02  7.2191e+03  3.1619e-05  4.22e-09
1.0e+01  2.2137e-07  2.6280e-02  2.2834e+03  9.9965e-05  4.22e-09
1.0e+02  2.2137e-07  2.6280e-02  7.2380e+02  3.1537e-04  4.22e-09
1.0e+03  2.2137e-07  2.6283e-02  2.3426e+02  9.7452e-04  4.22e-09
1.0e+04  2.2094e-07  2.6582e-02  8.9577e+01  2.5775e-03  4.21e-09
1.0e+05  1.9891e-07  4.3432e-02  5.7789e+01  6.5280e-03  4.00e-09
1.0e+06  1.5838e-07  1.2883e-01  4.7210e+01  2.3702e-02  3.57e-09
1.0e+07  1.4489e-07  3.9765e-01  4.3407e+01  7.9572e-02  3.41e-09
1.0e+08  1.4061e-07  1.2483e+00  4.2151e+01  2.5724e-01  3.36e-09
1.0e+09  1.3926e-07  3.9386e+00  4.1748e+01  8.1944e-01  3.34e-09
1.0e+10  1.3883e-07  1.2446e+01  4.1620e+01  2.5974e+00  3.34e-09
1.0e+11  1.3869e-07  3.9349e+01  4.1579e+01  8.2199e+00  3.34e-09

When L and R are split into components for the inner conductor, the air gap and outer tube:


AC-values:
f        L_inner     L_gap       L_outer       R_inner     R_outer
[Hz]     [H/m]       [H/m]       [H/m]         [ohm/m]     [ohm/m]
1.0e+00  5.0000e-08  1.3863e-07  3.2741e-08    2.1900e-02  4.3799e-03
1.0e+01  5.0000e-08  1.3863e-07  3.2741e-08    2.1900e-02  4.3799e-03
1.0e+02  5.0000e-08  1.3863e-07  3.2741e-08    2.1900e-02  4.3800e-03
1.0e+03  4.9998e-08  1.3863e-07  3.2739e-08    2.1901e-02  4.3815e-03
1.0e+04  4.9830e-08  1.3863e-07  3.2480e-08    2.2049e-02  4.5325e-03
1.0e+05  3.9188e-08  1.3863e-07  2.1089e-08    3.1782e-02  1.1650e-02
1.0e+06  1.3152e-08  1.3863e-07  6.5956e-09    8.8688e-02  4.0138e-02
1.0e+07  4.1732e-09  1.3863e-07  2.0871e-09    2.6786e-01  1.2979e-01
1.0e+08  1.3201e-09  1.3863e-07  6.6006e-10    8.3496e-01  4.1336e-01
1.0e+09  4.1746e-10  1.3863e-07  2.0873e-10    2.6285e+00  1.3101e+00
1.0e+10  1.3201e-10  1.3863e-07  6.6006e-11    8.3001e+00  4.1459e+00
1.0e+11  4.1746e-11  1.3863e-07  2.0873e-11    2.6235e+01  1.3114e+01



7. Bessel-bignum

Sometimes we want to compute Bessel funtions with better precision. Particularly when the Hankel functions are computed the numerical loss of precision is annoying for some arguments. I have a program called bessel_gmp.tgz (14 k), that could be used to compute Bessel functions to arbitrary precision. It is based on the surprisingly fast GNU MP Library. Note that bessel-bignum also include a reasonably fast routine to compute pi (3.14159....) to arbitrary precision. I have used it to get pi to one million places. It takes about half a minute. To compute Bessel functions to one million places would take a very long time, but to get 1000 takes less than a second.

Note: If you would like to do programming based on bessel-bignum, please first take a look at e.g. MPFR , a library containing a lot of useful functions. Most of them faster and more refined than the slower counterparts used in bessel-bignum. However, to my knowledge, it does not contain Bessel functions, but faster log, exp and euler-gamma.


E-mail: <sven.wilhelmssonNOSPAMswipnet.se>
Back to home-page, Updated Apr 30 2010