MATLAB Examples

SRACS example

This file demonstrates the usage of SRACS (self resonant air coil
solver). You should use the 'Publish' functionality of matlab, or 'run
and advance', to get through this file step-by-step.
Copyright 2015 Bernd Breitkreutz

Contents

Designing an Archimedean spiral coil

A simple coil geometry is the Archimedean spiral. Its windings are in one plane, and they all have the same distance.

We define the geometry of the coil by storing its properties in a struct. We call it arch, for Archimedean spiral.

A coil has a number of windings, W, and a diameter of the wire, diameter:

arch=struct; %overwrites any existing variable called arch
arch.W=6.33;
arch.diameter=1e-3;

The exact shape is defined by two functions, r(phi) and z(phi). phi is the winding angle. It starts at 0 and ends at W*2*pi. Therefore, the positions phi0 and phi0+n*2*pi are 'neighbors'.

In SRACS, r and z are defined by a function [r,z]=w_description(phi,parameters). To describe an Archimedean spiral, we take the function w_linear. We chose a winding distance of 1cm in r-direction, 0cm in z-direction (because it is not a cone or a screw), and an inner diameter of 2cm.

arch.ana.w_function=@w_linear;
arch.ana.w_params={1.e-2, 0.e-2, 2.e-2};

Now that the geometry is defined, we have to take care of a sufficient discretization for the numeric computation. We divide the spiral in N segments.

arch.N=10;

SRACS treats every segment as a set of concentric rings to define the electric properties. The wire is divided in W*N pieces, so make sure that W*N is sufficiently large to display the current distribution of the highest mode number you are interested in. We have 63 wire segments, so the 6th mode should be still okay.

You might want to use a high value N for smoother plots. But even the low N solutions, wich have a very edgy plot, are still pretty accurate. This is because SRACS calculates with pieces of rings, and not straight lines.

Computation of the resonances

Now we can start the simulation:

tic
archOut=sracs_solver(arch);
toc
display(sprintf('The spiral has %g windings.', archOut.W))
display('First seven resonant frequencies:')
display(sprintf('%f MHz\n',  archOut.sol.f(1:3) / 1e6))
Elapsed time is 0.162679 seconds.
The spiral has 6.3 windings.
First seven resonant frequencies:
57.861352 MHz
136.709766 MHz
212.952260 MHz

The elapsed time should be less than a second. Try some different values for N and see how the calculation time and the resonant frequencies change. Important is the fact that the value of W may change. It is rounded to a multiple of 1/N. Therefore, if you have a non-integer amount of windings, make sure you have a finer discretization.

arch=archOut;
arch
arch = 

           W: 6.3000
    diameter: 1.0000e-03
         ana: [1x1 struct]
           N: 10
         set: [1x1 struct]
         geo: [1x1 struct]
         num: [1x1 struct]
         sol: [1x1 struct]

The struct extended by a couple of new sub-structs. In set, there are some default values for calculation options.

arch.set
ans = 

          solNumber: 1
          MagSolver: 'noprox'
          EleSolver: 'noprox'
    MagSolverParams: 4
    EleSolverParams: 4

You may want to change your arch.set.solNumber to the mode you are interested in. This will be important for the quality factor.

arch.geo
ans = 

     phiI: [1x64 double]
       rI: [1x64 double]
       zI: [1x64 double]
     phiV: [1x64 double]
       rV: [1x64 double]
       zV: [1x64 double]
       Wh: [7 7 7 6 6 6 6 6 6 6]
    rImat: [7x10 double]
    zImat: [7x10 double]
    rVmat: [7x10 double]
    zVmat: [7x10 double]
     wire: [1x64 double]

In arch.geo is the discretized geometry. Current and voltage are not defined at the same positions, this is why we have I and V values. The variable wire consists of the length of the wire from the start to each point where the current is defined.

arch.num is only relevant for the elecromagnetic properties of the single wire pieces.

The results of the calculation are stored in arch.sol

arch.sol
ans = 

    f: [63x1 double]
    I: [64x63 double]
    V: [64x63 double]

arch.sol.f are the frequencies of the single modes in Hz, and arch.sol.I and arch.sol.V the corresponding currents and voltages.

But first, lets have a look at the geometry and the current distribution

sracs_plot_geo(arch)

We see the current of the arch.set.solNumber's mode, which should still be 1. The ends of the wire are black (no current), and the middle is red (maximum current).

sracs_plot_geo(arch, 3)

The third mode has one maximum and two minimums, which are displayed in blue.

sracs_plot_sol(arch, 2)

This function displays the current and voltage along phi. If you are interested in the distribution against the wire position, you have to plot it manually:

figure
plot(arch.geo.wire, arch.sol.I(:,2))
xlabel('position on wire in meters')
ylabel('current amplitude of second mode')

You may notice the sharp edge in the chart. This is because SRACS neglects stray capacitances at the ends of the wire. However, the effect on the frequencies and quality factors is rather small.

more complicated geometries

A strength of SRACS is the flexibility of the winding function. You can define pretty complicated resonators, as long as they are 'spiralish'. Don't try to define a screw with only 3 windings and a diameter of 10cm, that is a couple of meters long. This is a dipole antenna and not a coil!

Lets have a look at a ballish resonator:

ball=struct;
ball.N=20; %slower, but smooth plots
ball.W=10;
ball.diameter=1e-3;
ball.ana.w_function=@w_christmasTreeBall;
ball.ana.w_params={1e-2};

ball=sracs_solver(ball);
sracs_plot_geo(ball,2);

This is a 3D-geometry, but still pretty simple, because its windings are sorted (i.e. the second winding is between the first and the third one). Let's now put the windings of a spiral around the surface a torus, and let us surround it twice:

torus=struct;
torus.N=30;
torus.W=14.75;
torus.diameter=1e-3;
torus.ana.w_function=@w_torus;
torus.ana.w_params={10e-2, 2e-2, 2};
torus=sracs_solver(torus);
sracs_plot_geo(torus,2)

Now the windings are geometrically mixed, inner and outer ones close together.

Tuning

Usually one is interested given resonant frequencies. SRACS comes with a function to tune a spiral via the number of windings. One should define the winding function in a way that tuning leads to the desired change. In this example, we tune to 50MHz and the outer radius of the spiral will become bigger.

archTuned=arch;
archTuned.set.solNumber=1;
archTuned.set.freqGoal=50e6;
archTuned.set.freqTol=0.01;
archTuned=sracs_tuner(archTuned,1);

display(sprintf('\nOriginal: N=%i, W=%2.2g, f=%2.3gMHz', arch.N, arch.W, arch.sol.f(1)/1e6))
display(sprintf('Tuned: N=%i, W=%2.2g, f=%2.3gMHz', archTuned.N, archTuned.W, archTuned.sol.f(1)/1e6))
small N -> look for interval. W=3, 10, 20, 30, ... 
W=10
N=4----------------------------------------------------
Start bisection
W1=6_1/4  f1=58.5888MHz  |  W2=10_0/4  f2=26.5955MHz 
W1=6_3/4  f1=51.6499MHz  |  W2=10_0/4  f2=26.5955MHz 
Result
W1=6_3/4  f1=51.6499MHz  |  W2=7_0/4  f2=48.6774MHz 
0.5*delta_f/fGoal=0.029725 Tol=0.01
Error=0.032998   -0.026452
N=17----------------------------------------------------
Verify interval
W1=6_13/17  f1=51.4968MHz  |  W2=7_0/17  f2=48.7040MHz 
Start bisection
W1=6_15/17  f1=50.0680MHz  |  W2=7_0/17  f2=48.7040MHz 
Result
W1=6_15/17  f1=50.0680MHz  |  W2=6_16/17  f2=49.3781MHz 
0.5*delta_f/fGoal=0.0068982 Tol=0.01
Error=0.0013592   -0.012437

Original: N=10, W=6.3, f=57.9MHz
Tuned: N=17, W=6.9, f=50.1MHz

Quality factor

Now that we have the resonant frequencies and current distributions (and the electric properties in '.num') We can calculate the quality factor. Because it can be a little time consuming, it is only computed for one single mode. We will choose the first one:

arch.set.solNumber=1;

The qualitiy factor depends on the conductivity of the metal. We use copper:

arch.conductivity=58e6;

tic
arch=sracs_quality(arch)
toc
Warning: No proximity method DUT.set.ProxSolver was set. Default is
WireModel_extrapol 

arch = 

               W: 6.3000
        diameter: 1.0000e-03
             ana: [1x1 struct]
               N: 10
             set: [1x1 struct]
             geo: [1x1 struct]
             num: [1x1 struct]
             sol: [1x1 struct]
    conductivity: 58000000
              pp: [1x1 struct]

Elapsed time is 0.459077 seconds.

The struct was extended by a new sub-struct, which contains the postprocessing results for arch.set.solNumber, i.e. the stored energy, radiated and heat power and the quality factors with respect to ohmic losses, radiation losses and both of them.

arch.pp
ans = 

         W: 9.9375e-07
      Prad: 0.0327
    Pohmic: 0.3297
        Q0: 1.0957e+03
      Qrad: 1.1032e+04
      Qtot: 996.7087

The settings have been extended as well:

arch.set
ans = 

           solNumber: 1
           MagSolver: 'noprox'
           EleSolver: 'noprox'
     MagSolverParams: 4
     EleSolverParams: 4
          ProxSolver: 'WireModel_extrapol'
    ProxSolverParams: [20 21]

To calculate the quality factor, the skin-effect is always considered. But if the single windings are close to each other, the proximity-effect is significant as well. The defaults consider the proximity-effect, that is why the computation of Q is a little slow. If the gap between your windings is sufficiently huge, you can change the options to 'noprox' and 4 to speed up the computation.

Coupled resonators

SRACS is able to calculate systems of inductively(!) coupled resonators. Capacitive coupling is not implemented, because the fields are mostly concentrated between the windings.

Be aware that the magnetic coupling is implemented with a simple Neumann-formula. Therefore, you will propably need bigger N than for the uncoupled frequencies. Please try some different values and have a look at the convergence.

We define a new struct that defines a system of coupled coils. We want to know the coupling coefficient between two different Archimedean spirals, but with same resonant frequency.

arch.set.solNumber=1;
arch.set.freqGoal=50e6;
arch.set.freqTol=0.01;
arch=sracs_tuner(arch);

arch2=arch;
arch2.ana.w_params={1e-2, 0.3e-2, 8e-2};
arch2=sracs_tuner(arch2);

twoSpirals=struct;
twoSpirals.res=[arch arch2];

The first spiral stays at the origin of the coordinate system (x,y,z). The second one is shifted by 5cm, and tilted by 20 degrees.

twoSpirals.res(2).pos.z=10e-2;
twoSpirals.res(2).pos.beta=-20;

tic
twoSpirals=sracs_system_solver(twoSpirals)
toc
twoSpirals = 

      res: [1x2 struct]
      nor: 2
     segs: [117 81]
    first: [1 118]
     last: [117 198]
     dphi: [0.3696 0.2618]
      geo: [1x1 struct]
      num: [1x1 struct]
      sol: [1x1 struct]

Elapsed time is 0.085568 seconds.

The results are similar to those of a single spiral, but contain the solutions for all resonators. Therefore, the first and last index of every spiral is given.

A look at the results:

sracs_system_plot_geo(twoSpirals)
sracs_system_plot_geo(twoSpirals,2)
sracs_system_plot_sol(twoSpirals, [1,2])

Now that we have a system of two resonators, the first two eigen solutions represent the odd and even mode of the first self resonant frequency. We can easily calculate the coupling coefficient from the frequency values:

f1=twoSpirals.sol.f(1);
f2=twoSpirals.sol.f(2);

k=(f2^2-f1^2)/(f2^2+f1^2);
display(sprintf('The coupling is k=%1.3g', k))
The coupling is k=0.0873

Of course, we are not limited to two spirals.

N=8;
R=15e-2;

arch.pos.y=0;
arch.pos.z=0;
arch.pos.beta=0;

sys=struct;

for ii=1:N
    phi=(ii-1)*2*pi/N;
    sys.res(ii)=arch;
    sys.res(ii).pos.y=R*sin(phi);
    sys.res(ii).pos.z=R*cos(phi);
    sys.res(ii).pos.beta=phi/pi*180+90;
end

sys=sracs_system_solver(sys);
sracs_system_plot_geo(sys,3)