Skin effect for the single long copper conductor

Theoretical background

Based on Maxwell equations for quasi-stationary fields (sinusoidal currents and voltages) one can derive the following formulas:

$\mathrm{rot} \mathbf{\overline{\underline{H}}} = \mathbf{\overline{\underline{J}}}$ (1)

$\mathbf{\overline{\underline{J}}} = \gamma \cdot \mathbf{\overline{\underline{E}}}$ (2)

$\mathrm{rot} \mathbf{\overline{\underline{H}}} = \gamma \cdot \mathbf{\overline{\underline{E}}}$ (3)

$\mathrm{rot} \mathbf{\overline{\underline{E}}} = - \dfrac{\partial \mathbf{\overline{\underline{B}}}} {\partial t} = -j \omega \mathbf{\overline{\underline{B}}} $ (quasi-stationary field) (4)

Performing rotation on both sides of equation (3):

$\mathrm{rot} (\mathrm{rot}\mathbf{\overline{\underline{H}}}) = \gamma \cdot \mathrm{rot}\mathbf{\overline{\underline{E}}}$ (5)

$\mathrm{grad} (\mathrm{div} \mathbf{\overline{\underline{H}}}) - \nabla^2 \mathbf{\overline{\underline{H}}} = -j \gamma \omega \mathbf{\overline{\underline{B}}}$ (6)

remembering that:

$\mathbf{\overline{\underline{B}}} = \mu \mathbf{\overline{\underline{H}}} $ (7)

$\mathrm{div} \mathbf{\overline{\underline{B}}} = \mathrm{div} \mathbf{\overline{\underline{H}}} = 0 $ (8)

We get:

$\nabla^2 \mathbf{\overline{\underline{H}}} = j \gamma \omega \mu \mathbf{\overline{\underline{H}}} $ (9)

Defining k as:

$k = \sqrt{j \gamma \omega \mu } $ (10)

The final equation takes form known as Helmholtz equation:

$\nabla^2 \mathbf{\overline{\underline{H}}} = k^2 \mathbf{\overline{\underline{H}}}$ (11)

Assumptions

The conductor in the analyzed case can be represented as:


Assumptions based on this very case are the following:

  • The conductor is very long, which means that d/dz = 0.
  • Conductor width a is much smaller than its height b, which means d/dy = 0.
  • Current flows in the direction of the axis z.
  • The specific solution of the Helmholtz equation in such case can be written with the following equation:

    $\mathbf{\underline{H(x)}} = \mathrm{A}e^{kx} + \mathrm{B}e^{-kx} $ (12)

    Coefficients A and B can be obtained from the boundary conditions.

    Boundary conditions

    The red lines on the picture below represent the lines the magnetic field. The red arrows show the direction of the field.

    One can see that the first boundary condition can be written as:

    $\mathbf{\underline{H}} \left( \frac{-a}{2} \right) = - \mathbf{\underline{H}} \left( \frac{a}{2}\right)$

    Using equation 12, one gets:

    $\mathrm{A}e^{-k\frac{a}{2}} + \mathrm{B}e^{k\frac{a}{2}}= \mathrm{-A}e^{k\frac{a}{2}} - \mathrm{B}e^{-k\frac{a}{2}}$

    After cleaning up the equation can be rewritten as:

    $\mathrm{A}\left(e^{-k\frac{a}{2}} + e^{k\frac{a}{2}}\right)= \mathrm{-B}\left(e^{-k\frac{a}{2}} + e^{k\frac{a}{2}}\right)$

    which in the end gives:

    $\mathrm{A} = \mathrm{-B}$ (13)

    Adding 13 to 12 gives:

    $\mathbf{\underline{H(x)}} = \mathrm{A}\left(e^{kx} - e^{-kx} \right) $

    Rembering that :

    $\sinh{(x)} = \frac{e^{x} - e^{-x}} {2}$

    We finally get:

    $\mathbf{\underline{H(x)}} = 2\mathrm{A}\sinh{(kx)}$ (14)

    The second boundary condition can be written that at x = a/2 the magnetic field has a specific value:

    $ \mathbf{\underline{H}} \left( \frac{a}{2}\right) = \mathbf{H_0}$ (15)

    Using Ampere's law we can write:

    $ \oint_l \mathbf{\overline{\underline{H}}} \circ \overline{dl} = \underline{I} = H_0 \cdot 2b$ (16)

    where:
    $\underline{I}$ is a current flowing through the conductor
    $H_0$ is magnetic field intensity value at the surface of the conductor at x=a/2

    Knowing that one can rewrite equation 15 to:

    $ \mathbf{\underline{H}} \left( \frac{a}{2}\right) = \mathbf{H_0} = \frac{\underline{I}}{2b} = 2\mathrm{A}\sinh{\left(k \frac{a}{2} \right)}$

    which finally gives:

    $\mathrm{A} = \frac{\underline{\displaystyle I}}{\displaystyle 4b \sinh{\left(k\frac{a}{2}\right)}}$ (17)

    and which leads to:

    $\mathbf{\underline{H(x)}} = \frac{\underline{\displaystyle I}}{\displaystyle 2b} \frac{\displaystyle \sinh{(kx)}}{ \displaystyle \sinh{\left(k \frac{a}{2}\right)}}$ (18)

    Remembering (1):

    $\mathrm{rot} \mathbf{\overline{\underline{H}}} = \mathbf{\overline{\underline{J}}} = \overline{1_z} \frac{ \partial{\mathbf{\underline{\displaystyle H(x)}}}}{\partial{\displaystyle t}}$

    We can obtain the formula for the current density:

    $\mathbf{\underline{J(x)}} = \frac{\displaystyle k}{\displaystyle 2b} \underline{I} \frac{\displaystyle \cosh{(kx)}}{\displaystyle \sinh{\left(k \frac{a}{2} \right)}}$ (19)

    
    
    In [1]:
    #importing all required modules
    #important otherwise pop-up window may not work
    %matplotlib inline 
    import numpy as np
    import scipy as sp
    import matplotlib as mpl
    import matplotlib.pyplot as plt
    import math
    import cmath
    import seaborn
    
    
    
    
    /Users/rafal/anaconda/envs/py35/lib/python3.5/site-packages/matplotlib/__init__.py:872: UserWarning: axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter.
      warnings.warn(self.msg_depr % (key, alt_key))
    
    
    
    In [2]:
    #Defining conductor dimensions [mm]
    a = 10
    b = 100
    
    a = a/1000
    b = b/1000
    
    #Conductor material properties
    gamma = 58*10**6
    u = 4*math.pi*10**-7
    
    #frequency 
    f = 50
    omega = 2*math.pi*f
    
    #skin depth
    delta = math.sqrt(2/(omega*gamma*u))
    
    #parameter k
    k = cmath.sqrt(1j*omega*u*gamma)
    
    #Current flowing through conductor [A]
    I = 1000
    
    #defining x
    x = np.linspace(-a/2, a/2, 20)
    
    #Current density [A/mm2]
    J = [k/(2*b)*I*cmath.cosh(k*w)/cmath.sinh(k*a/2)/10**6 for w in x]
    
    
    
    In [5]:
    #definition of module of current density J (otherwise only real component will be plotted)
    modJ = [abs(J[i]) for i in range(0, len(J))]
    rJ = [J[i].real for i in range(0, len(J))]
    iJ = [J[i].imag for i in range(0, len(J))]
    
    #plot preparation
    fig, ax = plt.subplots(2, figsize=(10,10))
    
    ax[0].plot(x*10**3,modJ, label="|J| [A/mm2]", color = "green")
    ax[0].plot(x*10**3,rJ, label="Re{J} [A/mm2]", color = "red")
    ax[0].set_ylabel("Current density [A/mm2]")
    ax[0].set_xlabel("Conductor width [mm]")
    ax[0].set_title("Current density distribution along conductor width")
    #ax.axis([-5, 5, 0.9, 1.3])
    ax[0].legend()
    
    ax[1].plot(x*10**3,iJ, label="Im{J} [A/mm2]", color = "blue")
    ax[1].set_ylabel("Current density [A/mm2]")
    ax[1].set_xlabel("Conductor width [mm]")
    ax[1].legend()
    plt.tight_layout()
    
    
    
    

    Power loss calculation

    By definition power loss is given by the following formula:

    $P = \int_{V} \dfrac{J^2}{\gamma} dV$

    By performing the integral one gets:

    $P = \dfrac{lk^2}{4\gamma b} I^2 \dfrac{1}{sinh^2\left(k\dfrac{a}{2}\right)}\left(\dfrac{1}{2k}\sinh{(ka)}+\dfrac{a}{2}\right)$

    
    
    In [4]:
    #conductor length [m]
    l = 1
    
    #Power loss for ac current
    P = l*k**2/(4*gamma*b)*I**2/(cmath.sinh(k*a/2))**2*(1/(2*k)*cmath.sinh(k*a)+a/2)
    
    #Power loss for dc current
    R = l/(gamma*a*b)
    P_dc = I**2*R
    
    abs(P), P_dc
    
    
    
    
    Out[4]:
    (17.117004128257975, 17.24137931034483)
    
    
    In [ ]: