 # Design Simulation Systems Ltd

Modelling Cored Inductors And Transformers
User Manual & Tutorial

Mark Sitkowski

### Overview

A magnetomotive force (MMF) sends a flux (PHI) through a magnetic circuit. The MMF is produced by current in a coil and is measured in ampere turns:

```     MMF = 0.4PI * N * I
```

One effect of MMF is magnetising force, H, where:

```      H = MMF / LP
= (0.4PI * N * I) / LP
```

The other effect of MMF is flux, which is analogous to current, and from which the flux density is defined as:

```      B = PHI / A
```
whence
```      PHI = B*A
```

or, in saturation,

```     PHI ~ 2*B*A
```

For a given material, a magnetising force, H produces a magnetisation, B, and the slope of the curve is mu, where:

```      mu = B / H
```

In the linear region, B changes appreciably with H, and mu has its proper value, while in the saturation region, it has a much lower value, since B changes little with H.

The inductance of a cored coil is:

```     L = (mu*N^2*0.4PI*A*1e-8) / LP
```

The voltage across a coil is defined both in terms of flux and inductance:

```     E = L*di/dt = N*dPHI/dt
```

Our modelling approach is taken from the paper "A Computer Model of Magnetic Saturation and Hysteresis for Use on SPICE2", by Pei & Lauritzen. However, their method for extracting model parameters has been ignored, since it requires curve-fitting techniques, which are not easily implemented. A major factor in the problem, is the fact that they use a single coefficient/power to account for the entire shape of the hysteresis curve. This sometimes leads to large coefficients, (1000e+33) which cause convergence problems.

We have split the hysteresis function into three separate, easily characterized sub-functions:

• Core saturation current
• Linear permeability
• Saturation permeability

### Saturation Permeability

To model a cored coil, we model the saturation effect by using a feedback current proportional to an odd power of the coil current.

For a high-order polynomial, the output is near zero for values of input less than unity, and very large for values greater than unity. This means that there is little or no feedback for values up to +/-1, and progressively larger values outside this. As the feedback current increases, the inductor current saturates at a rate proportional to the slope of the polynomial.

The slope of the polynomial after unity is equal to the saturation permeability of the material, while the slope of the polynomial between the values +/-1 is the linear permeability.

Although the polynomial saturates at an input value of unity, this is not the case for any real components.

To relate the turnover points of the hysteresis loop to values of current other than unity, a scaling controlled source is used, and the actual saturation current (which determines H) determines the gain of the controlled source. Thus, if the core saturates at a current of 4 mA, the gain of the source is set at (1 / 4e-3) = 250.0

There is a small error relating to this, as the polynomial slope for unity input is quite shallow for low order polynomials, but if this is the case, (1.5 / 4e-3) may be used.

The following table may be used to select the order of the polynomial to match the saturation permeability of a given material.

NOTE: Since a hysteresis loop is an odd-order polynomial on its side, all references to "slope" in the following actually mean "inverse slope".

The values of corner slope, (the slope between x=0.9 and x=1.1) and the slope between x=2 and x=1.5 are given below for odd powers between the orders of 3 and 19 (SPICE will not accept a polynomial of order higher than 20):

 Power Corner Central slope Slope 3 3.0100 1.0E-02 9.2500 5 5.1001 1.0E-04 48.812 7 7.3521 1.0E-06 221.83 9 9.8526 1.0E-08 947.11 11 12.697 1.0E-10 3,923 13 15.990 1.0E-12 15,995 15 19.857 1.0E-14 64,660 17 24.438 1.0E-16 260,170 19 29.904 1.0E-18 1,044,100

An examination of some practical values for linear permeabilities follows:

 Material UINIT UMAX BSat(gauss) Ferrite - 2.7k 150 Mumetal 30k 130k 8.5k Orthonol - 75k 2k Permalloy C 16k 75k 8k Permalloy - 180k 10k Radiometal 2.2k 22k 16k Permalloy B 2k 15k 16k Permalloy A 12k 90k 11k Cr-Permalloy 12k 60k 8k Mo-Permalloy 20k 75k 8.5k 1040 40k 100k 6k Megaperm 3.3k 68k 9.3k Hipernick 3k 70k 15.5k 45 Permalloy 2.7k 23k 16k

For comparison, the following table lists the characteristics of some Ferroxcube materials. BSAT is measured in gauss, and HSAT in oersteds, while RHO (the resistivity) is in ohm-cm, and EPSILON (the dielectric constant) is measured at 1kHz and 10MHz.

 Material UINIT USAT BSAT HSAT RHO EPSILON A1/3b 500 340 3400 10 100/10 150k/50k A2/3c2 800 350 3500 10 100/10 150k/50k A3/3c 800 380 3800 10 100/10 150k/50k A4/3a 1000 360 3600 10 100/10 150k/50k B1/4a 650 290 2900 10 5e6/2e7 400/15 B2/4b 250 320 3200 20 5e6/2e7 400/15 B3/4c 110 90 2700 30 5e6/2e7 400/15 B4/4d 46 42 2500 60 5e6/2e7 400/15 B5/4e 18 24 1900 80 5e6/2e7 400/15

Thus, for permalloy, a 13th order polynomial gives approximately the correct slope, while for orthonol, an 11th order is near enough. Ferrite can be modelled with a 7th order term. Since these are all saturation values, the tolerance is acceptable.

### Linear Permeability

A single term odd-order polynomial gives a very steep characteristic for input values less than unity, since fractional values raised to an odd power change little, so it is necessary to add a term proportional to the input. This increases the slope of the saturation characteristic as the following table shows for a unity linear term of form:

```     (poly(1)  0  1.0  0....  0  1.0):
```
 Power Corner Central slope Slope 3 4.0100 1.0100E+00 10.250 5 6.1001 1.0001E+00 49.812 7 8.3521 1.0000E+00 222.83 9 10.853 1.0000E+00 948.11 11 13.697 1.0000E+00 3,924 13 16.990 1.0000E+00 15,996 15 20.857 1.0000E+00 64,661 17 25.438 1.0000E+00 260,170 19 30.904 1.0000E+00 1,044,100

As may be seen, a unity coefficient gives a central slope of unity, while not affecting the saturation slope by more than a few percent. The inverse of the value of the linear permeability then becomes the coefficient of the first term of the polynomial. For example, a material with a linear permeability of 2700 (Ferrite), would have its hysteresis loop defined by:

```     poly(1) 0 3.704e-4 0 0 0 0 0 1.
```

In terms of plotting characteristics, it is not necessary to plot B and H explicitly. From above, it may be seen that H is related to the current in the coil by the relationship

```     H = 0.4PI*N*I / LP
```
Thus, a plot of current will be equivalent to a plot of H multiplied by a constant scaling factor.

Likewise, the voltage across the coil is related to flux, by the relationship E = N*dPHI/dt, and B = PHI / A. Thus, again, a plot of the integral of the coil voltage gives B, multiplied by a scaling factor. Performing an X-Y plot of a sinusoidal stimulus against coil current gives the correct shape of hysteresis loop.

Now that we have a useable hysteresis loop, the magnetising inductance of the transformer may be calculated from the permeability figures, and this value assigned to the inductor in the transformer model.

### Losses

The permeability is a function of frequency. As the frequency goes up, the hysteresis loop gets fatter, which means that the central slope goes down, corresponding to a drop in permeability. The drop off is fairly linear, at approximately 3dB/decade. A feature of this which can be obtained from a data sheet is the frequency at which the permeability is 3dB down - U3DB.

The core losses are actually of the form:

```     Rm / Lm = Zh*Bm*f + Zc*f + Ze*f^2
```

where

• Bm is maximum flux density
• Rm is core resistance
• Lm is core inductance
• Zh is a hysteresis loss coefficient (typ: 0.16, 0.1, 0.03)
• Zc is a residual loss coefficient (typ: 3.0, 2.6, 1.0)
• Ze is a eddy current loss coefficient (typ: 2.5, 0.002, < 0.002)

If the frequency dependence of the losses is irrelevant, they may be quite simply modelled by a fixed resistor across the inductor.

Pei and Lauritzen model the losses as a parallel voltage-controlled current source, whose transconductance, Gm, is defined as:

```     Gm = (Z1 * V) + (Z3 * v^3) + (Z5 * v^5)...
```

where V is the voltage across the coil, and Z1..Z5 are coefficients.

For a period T, in steady state, the core loss can be calculated as:

```     Core loss = (1 / T) * (integral 0-T (V * Gm). dt)
= (0.5*Z1*V^2) + (0.375*Z3*V^4) + (0.3125*Z5*V^6)
```

The first coefficient, Z1, is already included in the definition of the hysteresis curve, as the maximum permeability, while the third coefficient can be omitted, as its effect is negligible.

The coefficients must be determined numerically from the magnetic core loss characteristic curves and are only good at one frequency, so they must be recalculated each time.

This places the approach on the same level as using a fixed resistor. However, the shape of the hysteresis loop is correct, and the following table lists the coefficients for ferrite, orthonol and permalloy:

 Material Freq Z1 Z3 Ferrite 50k 6.35e-6 5.67e-12 (G/F/J/W/P) 20k 1.61e-5 1.28e-10 10k 2.10e-5 5.23e-10 5k 2.77e-5 2.50e-09 Orthonol 10k 1.51e-5 -1.10e-11 (tape-wound 2mil) 6k 2.41e-5 -5.00e-11 3k 3.49e-5 -2.40e-10 1k 7.13e-5 -4.19e-09 0.4k 1.46e-4 -4.96e-08 Permalloy 50k 1.19e-4 2.88e-12 (tape-wound 2mil) 20k 2.49e-4 2.51e-11 10k 4.77e-4 1.28e-11 5k 8.03e-4 2.53e-09 2k 2.47e-3 9.16e-09

### Non-linear Transformers

Every paper published in this field shows that the transformer model used by SPICE is unsuitable for saturating-core type transformers. The preferred model is due to Middlebrook, which is the one we have adopted. The model itself is a DC coupled one, using a combination of voltage and current controlled sources, to model the coupling from primary to secondary, together with the relationships of the voltages and currents.

Onto this model may be added the equivalent circuit of a transformer, with everything referred to the primary. This gives excellent performance in both time and frequency domain, with the minor disadvantage of having DC coupling. An example of the model is shown in this library entitled "windings2", which is expanded by duplicating the secondary VCVS, and adding the extra secondary currents to the polynomial defining the primary current source. A three winding (one primary, two secondary) version is entitled "windings3".

### But first, Linear Transformers

Consider a two-winding transformer, as shown below:

```                   1_______   _______3
+     ) (     +
VALUE = 1mH               ) (           VALUE = 1mH
LOCATION = L1       V1    ) (     V2    LOCATION = L2
PARAMS = 0.8              ) (
COUPLED = L2              ) (
2_______) (_______4
1 : RATIO
```

This is created in GEX, using two inductors (see GSP manual for details of properties etc) and appears in the netlist as:

``` L1 1 2 1mH
L2 3 4 1mH
K1 L1 L2 0.8
```

What this means is that each winding is defined as an inductance, and the "K" line specifies the coupling between them, and the coefficient of coupling. The transformer's winding polarity is defined by the "+" on each inductor, which belongs to the first node number which, in the above example, is nodes 1 and 3.

This can be extended to handle transformers with more than two windings but it takes a little more work.

Since an example is worth ten pages of theory, consider the following three winding transformer, with unknown coupling coefficients:

```                   1_______   _______3
+     ) (     +
VALUE = 445u              ) (           VALUE = 653uH
LOCATION = L1       V1    ) (     V2    LOCATION = L2
PARAMS = A B C            ) (
COUPLED = L2 L3           ) (
2_______) (_______4
8 : 10
_______5
(     +
(           VALUE = 521uH
(     V3    LOCATION = L3
(
(
(_______6
8 : 300
```

It will appear in the netlist as:

``` L1 1 2 445uH
L2 3 4 653uH
L3 5 6 521uH
K1 L1 L2 A
K2 L1 L3 B
K3 L2 L3 C
```

Where A, B, C are the respective coupling coefficients. Note that you need three coupling definitions, to indicate that there is coupling between ALL the windings. If you omit one, you'll get misconvergence, or a one-way transformer.

To calculate the coupling coefficients from the transformer data, we need a minimum of the data given above.

For an ideal transformer:

```   V1 / V2 = N2 / N1    and    L1 / L2 = (N1 / N2)^2
```

The above is definitely not an ideal transformer, since

```    653u / 521u != (10 / 300)^2
```

Therefore, it is necessary to determine the actual coupling coefficients from the winding ratios and inductance values. To do this, we start with the basic differential equations for a transformer:

```   v1 = L1.di1/dt + M12.di2/dt
v2 = L2.di2/dt + M21.di1/dt
```

where the mutual inductance

```  M12 = M21 = k.sqrt(L1.L2)     (0 < k <= 1)
```

Therefore, replacing M12 M21 by M, and combining gives:

```   v1 = L1.di1/dt - M/L2.v2 + M^2/L2.di1/dt
```

Now, using v2 = v1(N2/N1) and tidying up:

```   v1 = [(L1 - M^2/L2) / (1 + M.N2/L2.N1)].di1/dt
```

Omitting the boring algebra, we finally arrive at:

```   k12 = k21 = N2/N1 * sqrt(L1/L2)
```

Recklessly substituting the above values:

```   k12 = 10 / 8 * sqrt(445u/653u) = 1.25 * 0.8255 = 1.032
```

WHAT?! That's not possible! We already agreed that k falls somewhere between zero and unity. The answer lies in the (boring) algebra we omitted. If this situation arises, we simply look for k21, instead of k12, which gives:

```  k21 =  8 / 10  * sqrt(653u/445u) = 0.8000 * 1.2114 = 0.969
```

Similarly:

```  k13 = 10 / 300 * sqrt(653u/531u) = 0.0333 * 1.1195 = 0.0373
k23 =  8 / 300 * sqrt(521u/445u) = 0.0267 * 1.0820 = 0.0288
```

Our transformer now appears like this:

``` L1 1 2 445uH
L2 3 4 653uH
L3 5 6 521uH
K1 L1 L2 0.969
K2 L1 L3 0.0373
K3 L2 L3 0.0288
```

Which just goes to prove that indiscriminately using a unity coupling factor can produce results which are wildly inaccurate. Most of the problems engineers have with transformer simulation lie in the fact that many of them don't go to enough trouble to get the correct values. This can lead to terminal convergence problems ("internal timestep too small..." remember?) or, worse, the simulation runs, but produces garbage. Each time, SPICE gets the blame. 