image thumbnail

Skineffect Calculation

by

 

10 Jul 2012 (Updated )

Calculation of integral parameters for round wires by theoretical electrical engineering.

skineffect

Contents

function tab = skineffect
% This is a demo function for calculations with the Bessel-Function in the
% elegtromagnetic science of cylindrical wires. This demo consist of 3 parts.
%
% In PART I the specific resistance and inductance parameters are calculated
% for a couple of frequencies for a copper wire. As well the skindepth and
% the ratios of different resistancemodels are also calculated. The values
% are written in a cellarray, which is the result [tab] of this function.
%
% In PART II the currentdensity is calculated as a function of radius of the
% wire for a frequency of 50 Hz and a current of 10 Amps. The current
% distribution is plotted as 3D-surface plot as well as over the radius.
%
% In PART III the electric field strength for the wire of a length of 1000m
% and the voltage drop are calculated. The parameter for resistance and
% inductance are also computed.
%
%
% Bugs and suggestions:
%    Please send to Sven Koerner: koerner(underline)sven(add)gmx.de
%
% Programmed by Sven Koerner: koerner(underline)sven(add)gmx.de
% Date: 2012/07/10

PART I - Calculate integral parameter of wire for different frequencies

% Inputdata for a copper wire of diameter 1mm
u0          = 4*pi*1e-7;   % Permeability constant in [Vs/Am]
ur          = 1;           % relative Permeability of Material   e.g. Copper
rho         = 1.72e-8;     % Resistivity in [Ohm*m],             e.g. Copper
D_wire_m    = 1e-3;        % Wire-Diameter 0.001m
f           = [1e-3, 1e0, 16+2/3, 50 , 60,  1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8]';       % frequency vector

% Calculate additonial parameter
R_m         = D_wire_m/2;  % Wire-Radius in [m²]
area        = R_m^2*pi;    % Circle_area in [m²]
sigma       = 1/rho;       % Conductivity, [S/m]

R_DC = 1/(sigma*area);     % DC-Resistance load per unit length in [Ohm/m]
omega  =  2*pi*f;                                  % angular frequency in [Hz]
delta  = 1./sqrt(omega.*sigma.*u0.*ur./2);         % skindepth in [m]

% calculate impedance of this wire;
z_math = Z_wire (omega, R_m,  sigma, u0, ur);

% R_DC_Ratio by an tubemodel-calculation (DC-Resistance in depence of skindepth and wire-radius)
R_dc_tube_ratio = tube_model_func(delta ,R_m ,sigma);


% *Collect the data for the table*
tab_header = {'frequency', 'skindepth', 'L lpul', 'R_DC lpul', 'R_DC-ratio ', 'R_DC-tubemodel_ratio'; ...
              'Hz',        'm',          'H/m',   'Ohm/m' ,    '[1]',          '[1]' };
tab2       = num2cell([ f, ...                                              % 1. column frequency
                        delta, ...                                          % 2. column skindepth
                        imag(z_math)./omega,  ...                           % 3. column impedance  (H/m)
                        real(z_math), ...                                   % 4. column resistance (Ohm/m)
                        real(z_math)/R_DC,...                               % 5. column resistance-Ratio [1]
                        R_dc_tube_ratio   ]);                               % 6. column resistance tube_model ratio [1]
tab       = [tab_header; cellfun(@num2str,(tab2), 'UniformOutput',0)];
% lpul: load per unit length

PART II - Plot currentdensity on wire for 50 Hz

f_circle_Hz = 50;     % frequency in Hz
current_A   = 1;      % current in wire in Amps
dr          = 200;    % number of steps (dicretization

% create a grid for that circle
[X,Y] = meshgrid(-R_m: R_m/(dr-1) :R_m , -R_m: R_m/(dr-1) :R_m  );

% Calculate radius for each Point:
r = sqrt(X.^2+Y.^2);

% calculation of currentdensity:
A      =  sqrt(-1j*2*pi*f_circle_Hz*sigma*u0*ur);       % Nomalized argument for Bessel function
[J0]   =  besselj ( 0 , A.* r );                        % first bessel-function of zeroth order
[J1]   =  besselj ( 1 , A.*R_m );                       % first bessel-function of first order
J_vec  =  A.*current_A.*1./(2.*pi.*R_m).*J0./J1;

% Recalc only the wire with radius
J_vec((X.^2+Y.^2)>=R_m^2)   = NaN;
X(isnan(J_vec))             = NaN;
Y(isnan(J_vec))             = NaN;

% Plot results
figure;
% Plot the Current density
surf(X,Y,real(J_vec),'LineStyle','none');
colorbar;
axis square
view([0 0 90]);
title (['Currentdensity in wire (R=', num2str(R_m), ' m, I=',num2str(current_A),' A, f=', num2str(f_circle_Hz), ' Hz)'  ]);

figure;
% Plot the Current density over radius
plot(Y(X==0), J_vec(X==0));
grid on;
xlabel('radius in m');
ylabel('J in A/m²');
title (['Currentdensity distribution over radius in wire (R=', num2str(R_m), ' m, I=',num2str(current_A),' A, f=', num2str(f_circle_Hz), ' Hz)'  ]);

PART III - Calculate the electrical field strength and the voltage drop for the parameters of PART II

length_m    = 1000;  % lenghth of wire;

E_Vpm       = current_A .* A./(2 .* pi.* R_m.* sigma)  .* besselj ( 0 , A.* R_m )./ besselj ( 1 , A.*R_m )   % electrical field_strength in Volt per m
deltaU_V    = E_Vpm.*length_m                              % complex voltage drop in Volt
absU_V      = abs(deltaU_V)                                % absolute voltage drop in Volt
R_Ohm       = real(deltaU_V/current_A)                     % Resistance in Ohm
Li_H        = imag(deltaU_V/(2*pi*f_circle_Hz*current_A))  % inductance in Henry
end

R_DC tubemodell (function)

function  R_dc = tube_model_func( delta ,r_outer, sigma )
  R_rod   = pi/2/sigma*r_outer^2;                       % rod resistance
  r_inner = max(r_outer - delta ,0);                    % inner radius of tube belong to skindepth
  R_tube  = pi./2./sigma.*(r_outer.^2 - r_inner.^2) ;   % tube resistance
  R_dc    = R_rod./R_tube;                              % Calculated DC-Resistance
end

Impedancemodell with Besselfun. (function)

function [z_math] = Z_wire (omega, R, sigma, u0, ur)
  A = sqrt(-1j*omega*sigma*u0*ur);              % Nomalized argument for Bessel function
  J0  =  besselj ( 0 , A.* R );                 % first bessel-function of zeroth order
  J1  =  besselj ( 1 , A.* R );                 % first bessel-function of first order
  z_math     = A./(sigma*2*pi.*R).*J0./J1;      % specific impedance
end

Contact us