OverviewThe techniques described below have now been superceded by a better approach. The tool Vcores is available for free download from the downloads page. The PDF manual is here 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 / Awhence 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*1e8) / 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 curvefitting 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 subfunctions:
Saturation PermeabilityTo 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 highorder 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 / 4e3) = 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 / 4e3) 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 oddorder 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):
An examination of some practical values for linear permeabilities follows:
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 ohmcm, and EPSILON (the dielectric constant) is measured at 1kHz and 10MHz.
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 PermeabilityA single term oddorder 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):
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.704e4 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 / LPThus, 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 XY 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.
LossesThe 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
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 voltagecontrolled 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 0T (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:
Nonlinear TransformersEvery paper published in this field shows that the transformer model used by SPICE is unsuitable for saturatingcore 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 TransformersConsider a twowinding 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 oneway 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.
