MATLAB Examples

## 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 ```