Exploiting the T(x) function in fast hysteresis models for transient circuit simulations

To consider non-linear inductors in circuit simulation, several approaches have shown good results in the past. The simplest implementation only covers saturation by using a non-linear BH-curves without any hysteresis effects. In (dos Santos et al., 2017) several mathematical models are compared and their accuracy is estimated. It is shown that exponential and hyperbolic functions are good candidates to represent the behavior of magnetic components.

By considering the magnetization such a simple model can be extended to include the magnetic hysteresis. The Jiles–Atherton model presented in (Jiles and Atherton, 1984) consists in its original form only of five parameters and gives reasonable results for major and minor loops. Here, the anhysteretic magnetization is calculated using a hyperbolic function – the Langevin function, shown in equation (2). This allows for much faster calculation than the classical Preisach model in (Preisach, 1935), where a large number of independent elements need to be evaluated.

Especially for numerical simulation, energy-based models can provide benefits. Both the Preisach and the Jiles-Atherton models are not naturally extending in three dimensions of space, and e.g. dissipated and stored energies are not available at every instance in time. A discussion of further differences is given in (Jacques et al., 2017).

The T(x) model on the other hand, extensively discussed in (Jenő’s (2003) study, is not fully based on existing theories. It can be understood as a very elegant mathematical description of the magnetic hysteresis. Therefore, its parameters are considered dimension-less and have no physical meaning.

While being not very well known, it is certainly very interesting for transient simulations. The Brillouin-function, B(x), as well as its simplified form, the Langevin-function, L(x), both suffer from poles at x = 0 as can be easily seen in equations (1) and (2). In contrast to these two functions, the formulation of the T(x)-function in equation (3) uses the hyperbolic tangents. Hence, here, x = 0 does not require any special treatment:

A comparison of all three functions is given in Figure 1. Here, for B(x) the constants C₁, C₂ and C₃ were set to 1. The constants of the other two functions were fitted to match the B(x) function as accurately as possible.

Another benefit of the T(x)-function is the existence of an exact inverse function.

Furthermore, it is shown (Jenő, 2003) that with proper choices of constants, the results of models based on the T(x)-function can be expected to closely match those of the widely accepted Preisach model.

The basic T(x) function is suitable for soft magnetic materials only. However (Dospial et al., 2014), shows an extension to hard magnetic materials by using T(x) function as a base function of a series.

The basic principles of this model are discussed using a normalized description of the magnetic hysteresis in order to make it easier for the reader to follow the equations. By shifting the T(x) function, the basic magnetic hysteresis loops are formed:

with

The functions f₊ and f − are the ascending and descending branch branches, respectively. x_m defines the minor loop and a₀ corresponds to the magnetic remanence. To obtain the virgin curve, simply, the mean value of these two functions is calculated:

An example of the resulting curves is given in Figure 2.

Reversal curves, needed to describe the hysteresis behavior in case of an arbitrarily changing exciting magnetic field, are more complex. Some examples for the resulting movements are shown in Figure 3.

Starting from x = x₀, the descending branch of the minor loop defined by x_m = 2 is followed. At x = x₁ the direction is reversed. The starting point of this new ascending curve needs to coincide with the reversal point at x = x₁. Furthermore, this new curve needs to reach the starting point x₀ at x = x_m. This can be accomplished by a scaled constant. Equation (4) will therefor get the form:

where c is the difference of the ascending and descending branch values at x = x₁ and the function H() is a scaling function ranging from 1 (at x = x_r) to 0 (at x = x_m) with the parameters x, x_r and x_m.

To define H() (Jenő, 2003), suggests a hyperbolic function.

This finally leads to:

with

This reversal curve can now be followed to x = x₂ = x₀ – the start of the minor loop. Increasing of x beyond x₂ to x = x₃, will cause a change to the virgin curve. Any change of direction on the virgin curve gives a new minor loop to be followed. However, if the reversal curve leading from x₁ to x₂ is interrupted at $x=x2′$⁠, a new value for x_m needs to be calculated to start moving on a new minor loop. To accomplish that, equation (10) has to be solved with regard to x_m:

To simplify the equations we define:

This now gives:

The minor loop defined by this new value of x_m leads to the points $x3′$ and $x4′$⁠. A much more detailed discussion of the construction of first- and higher-order reversal curves is given in chapter 5 of (Jenő, 2003). Below the T(x) model is applied to magnetic hysteresis. The magnetic flux density B corresponds to the value of the T(x) function whereas the magnetic field H is its parameter x (i.e. for the ascending branch simply B(H) = f₊ with x = H).

The circuit simulation code adopting the presented hysteresis model implements a simple, yet effective strategy.

Instead of a classic sparse matrix approach, as implemented in SPICE, each model is presented with vectors containing the node voltages (V), the residual currents (R), the results of an integration operation as well as vectors holding generic values representing e.g. model states. The model updates elements in R, representing the nodes it is connected to, and sets the derivative of its associated integration variables. After all models have been executed, each element of R represents the sum of all branch currents of the particular node. By the means of a quasi-newton algorithm on the node voltage vector, the residual current vector R is minimized using a least-squares approach to fulfill Kirchhoff’s current law. While losing the versatility of the classic SPICE approach (i.e. not all branch relations can be implemented), the source code is simplified a lot. However, considering real components instead of ideal ones, this shortcoming nearly vanishes since series resistors will be present in all models.

Integration methods implemented include forward-Euler, backward-Euler and bilinear algorithms. These can be selected at run-time. For the presented simulations, a bilinear integration was used for accuracy and stability reasons.

As discussed above, some branch relations cannot get implemented with the available simulation code. By adding a resistor, representing the losses in the copper wire, in series to the inductor, as shown in Figure 4, this problem is overcome. For the purpose of this paper, further parasitic properties, such as the capacitance of the winding and effects showing up at higher frequencies, are neglected.

The T(x) function provides a relation between the magnetic flux density B and the magnetic field strength H. For this reason, an application of Faraday’s Law of Induction presents a straightforward approach to model the inductor. Additionally, this gives useful results for the process of parameter identification from measurements as shown later on.

In datasheets of magnetic cores, usually the effective cross-section area A_eff and the effective length of the magnetic path l_eff are given. This allows for an application of the common approximations:

and

where Φ is the magnetic flux and N is the number of turns of the winding.

Applying these approximations to Faraday’s Law gives:

This result is, besides the use of $dBdH$ as the magnetic permeability μ, well known and represents the approximation applicable to long, thin inductors. However, the crucial parameters A_eff and l_eff are given in datasheets and the applications notes to different magnetic materials point toward using the given approximations. This suggests that model errors were already considered in these parameters. Equation (15) is further rewritten to give the current i in its differential form:

By utilizing the integration facilities of the simulator, the current I through the inductor can be obtained for every time-step. However, in order to do so, the term $dBdH$ needs further discussion, since, for obvious reasons, each curve will need to be treated differently. In contrast to the discussion above, the presented model implements the not normalized form of the T(x) function:

with

Here, B(H)_asc, B(H)_des, B(H)_vir represent the ascending and descending branch as well as the virgin curve, respectively. The constants A₀, a_o, B₀, C₀ have no direct physical meaning. Their values are results of a parameter fitting process using measurement values. In contrast to the previous constants, H_m is changing during the simulation and defines the maximum of the magnetic field strength of the minor loop.

From these definitions, the needed derivatives can be obtained. This gives for B(H)_vir:

As discussed above, the functions defining the ascending and descending branches need to be modified to represent reversal curves. Applying the introduced constant c and the scaling function H() to the not normalized functions gives:

and

This adds the constants c, the difference of the magnetic flux density between the ascending and descending branch at the reversal, and H_r, the magnetic field strength at the reversal, to the list of changing constants during the simulation.

Calculating the term $dBdH$ for these branch functions gives:

and

To finally calculate the derivative of the current through the inductor, the voltage across it needs to be defined. Since the current I, the result of the last integration operation, is an input parameter to the model, we can simply apply Kirchhoff’s voltage law to define the last unknown v_l:

The different constants used in this model can be separated into two categories:

1. constants changing during execution; and

2. real constants.

Real constants are the parameters of the T(x) function (A₀, a_o, B₀, C₀) as well as physical parameters of the core and winding (R, N, l_e, A_eff). In contrast to these, on certain events, the magnetic field strength defining the minor loop (H_m), the magnetic field strength at reversal (H_r) as well as the constant used to shift the reversal curve (c) are changing. These conditions need to be discussed next since they also dictate the choice of the function for the term $dBdH$⁠.

The behavior during the simulation is governed by a finite state machine. Its basic structure is shown in Figure 5.

The three states represent the virgin curve (𝕍) as well as the ascending (𝔸) and descending (𝔻) branches of the hysteresis loop. When entering one of these states, the function associated with the particular state is selected for the term $dBdH$⁠. The edges between these states represent changes of direction and changes of exciting magnetic field strength H_m defining the minor loop.

On the virgin curve, the current must monotonically increase or decrease. Otherwise, the state changes to either the descending (edge 1) or the ascending (edge 2) branch. This change is indicated by a change in the sign of $didt$ with respect to the current I. For positive currents, a negative $didt$ would call for a change of the state while for negative currents a positive $didt$ would do so. Since $didt$ is one of the results calculated by the model, this check can be efficiently implemented by simply evaluating $I·didt>0$⁠. When changing to a different state, the current magnetic field strength gives a new value for the variable H_m.

The edges 3 and 5 represent changes back to the virgin curve. As explained above this only occurs if |H| > |H_m|. Both values are available to the model.

Changes from the ascending to the descending branch and vice versa, are similar to the virgin curve. Again, evaluating $I·didt>0$ indicates a change along the edges 4 (ascending to descending) and 6 (descending to ascending). If a change happens before reaching H_m, the constants c and H_r need to be updated. This has already been explained above for equation (9).

Figure 5 suggests that the model can start not only at the virgin curve but also has the ability to start on the ascending and descending branches as well. This is useful to take, e.g. the residual magnetization into account. In this case, the initial minor loop parameters are calculated from the inverse functions of ascending and descending branches during the initialization phase of the model. For positive residual flux density, the simulation starts on the descending branch and for negative ones on the ascending branch respectively.

Therefore H_m for the state D in case of positive residual magnetization is:

A negative magnetization results in a start in state A and H_m gets:

To verify the presented model, a small inductor was prepared consisting of two ungapped E/25/13/7 cores out of 3C94 material and 10 turns of 1 mm magnet wire. Additionally, a winding consisting of 3 turns of RG213 coaxial cable was used to indirectly measure the magnetic flux density. The outer conductor of the coaxial cable was connected to ground on one end to form an electrostatic shield. Connecting the winding formed by the inner conductor to an R-C network forming a lowpass filter, the magnetic flux density can be measured after some considerations.

The construction of the inductor is shown in Plate 1. It is to be noted that in contrast to image, the inductor used for the simulations discussed below has no air gap. A 3d-printed frame compresses both halves of the core to ensure repeatable measurement results.

Faraday’s law of induction suggests that no constant magnetic field strength can be measured with this technique. In addition to that, the integrator, formed by the low pass filter introduces a phase error. The number of turns, the cross-section area of the core and the attenuation of the filter need to be considered to calculate the magnetic field B from this signal. To minimize the errors, the load at the measurement winding needs to be minimized; hence, the resistor value in the low pass filter needs to be as high as possible.

For practical reasons, the lowpass was built from a 1MΩ resistor and a 10nF capacitor. Bigger resistors need special care since contaminations of the circuit due to e.g. touching it may result in huge measurement errors. The capacitor was chosen so minimize the phase error, thus this was further neglected. Another criterion to be considered for the capacitor it the voltage dependency of the dielectric material. Here a C0G type capacitor was used. An amplification by a factor of 100 was used to a get signal to noise ratio of the measured signal. The circuit is shown in Figure 6. A 4-quadrant linear amplifier, capable of providing up to 40 A and 10 V, was used to excite the inductor. The test signals were generated by a programmable arbitrary-waveform-generator (AWG).

With the effective cross section area A_eff given in the datasheet of the core we get:

where Φ̂ is the signal after the integration and C represents the combined attenuation of the lowpass filter, the amplifier and the number of turns of the measurement winding. In addition to the voltage at the measurement winding, the current through the inductor was measured using a current clamp. The magnetic field strength is calculated using:

where l_e is the effective length of the magnetic path and N is the number of turns of the main winding.

Both signals were captured by a digital oscilloscope while applying a sinusoidal current, high enough to fully saturate the inductor. The model parameters for the branches of the major loop were fitted to these results. Some experiments showed that a genetic algorithm yields to good results while being very robust against the introduced measurement errors and noise. The choice of the optimization algorithm is uncritical since the parameters for the simulation are only calculated once. Figure 7 shows the measured data as well as the fitted major loop. The last missing model parameter, the series resistance, was measured using an LCR meter.

In addition to the current through the main winding of the inductor, the voltage across its terminal was also measured when collecting data for the parameter identification process. The measured current was post-processed to remove as much measurement noise as possible. To do so, first, a median filter over 11 samples was used to reduce outliers. After that, a first-order lowpass filter with a corner frequency of 10 times the frequency of the exciting current was applied. In the simulation, the measured current was injected into the inductor. Between the samples, an interpolation using a Lanczos filter kernel was used. This was done to further reduce the effects of measurement noise and to prevent oscillations in the simulation code (Figure 8).

As shown in Figure 10, simulation and measurement show a very good agreement. The noise on the simulated waveform, especially visible in the non-saturated region, was tracked down to the remaining noise in the measured current waveform.

In a further test, the more complex waveform in Figure 9 was used to excite the inductor. The previous simulation showed that measurement noise greatly degrades the simulation results. Therefore, in this test, the waveform data provided to the AWG was used in the simulation instead of the measured data. As a result, the simulation in Figure 10 has less noise and runs faster.

To show the capability of simulating minor loops, the major and a minor loop was extracted from the transient simulation results. Figure 11 shows the BH-curves and the corresponding measurements. Here, verification was not possible due to the measurement noise.

This paper demonstrates the feasibility of using the T(x) function in transient simulations.

It was shown that a rather simple model can be used to achieve simulation results in a good agreement to measurements.

An implementation of a more generalized form to support hard magnetic materials, as presented in (Dospial et al., 2014), an investigation whether the presented model also leads to accurate estimations of losses in the magnetic core and benchmarking against other hysteresis models present topics for future work.