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):

PowerCornerCentral slopeSlope
3 3.01001.0E-029.2500
55.10011.0E-0448.812
77.35211.0E-06221.83
99.85261.0E-08947.11
1112.6971.0E-103,923
1315.9901.0E-1215,995
1519.8571.0E-1464,660
1724.4381.0E-16260,170
1929.9041.0E-181,044,100

An examination of some practical values for linear permeabilities follows:

MaterialUINITUMAXBSat(gauss)
Ferrite-2.7k150
Mumetal30k130k8.5k
Orthonol-75k2k
Permalloy C16k75k8k
Permalloy-180k10k
Radiometal2.2k22k16k
Permalloy B2k15k16k
Permalloy A12k90k11k
Cr-Permalloy12k60k8k
Mo-Permalloy20k75k8.5k
104040k100k6k
Megaperm3.3k68k9.3k
Hipernick3k70k15.5k
45 Permalloy2.7k23k16k

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.

MaterialUINITUSATBSATHSATRHOEPSILON
A1/3b 500340340010100/10150k/50k
A2/3c2 800350350010100/10150k/50k
A3/3c 800380380010100/10150k/50k
A4/3a 1000360360010100/10150k/50k
B1/4a 6502902900105e6/2e7400/15
B2/4b 2503203200205e6/2e7400/15
B3/4c 110902700305e6/2e7400/15
B4/4d 46422500605e6/2e7400/15
B5/4e 18241900805e6/2e7400/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):
PowerCornerCentral slopeSlope
34.01001.0100E+0010.250
56.10011.0001E+0049.812
78.35211.0000E+00222.83
910.8531.0000E+00948.11
1113.6971.0000E+003,924
1316.9901.0000E+0015,996
1520.8571.0000E+0064,661
1725.4381.0000E+00260,170
1930.9041.0000E+001,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:

MaterialFreqZ1Z3
Ferrite50k6.35e-65.67e-12
(G/F/J/W/P)20k1.61e-51.28e-10
10k2.10e-55.23e-10
5k2.77e-52.50e-09
Orthonol10k1.51e-5-1.10e-11
(tape-wound 2mil)6k2.41e-5-5.00e-11
3k3.49e-5-2.40e-10
1k7.13e-5-4.19e-09
0.4k1.46e-4-4.96e-08
Permalloy50k1.19e-42.88e-12
(tape-wound 2mil)20k2.49e-42.51e-11
10k4.77e-41.28e-11
5k8.03e-42.53e-09
2k2.47e-39.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.