The Use of Hypergeometric Functions in Hysteresis Modeling

Abstract

Accurate hysteresis models are necessary for modeling of magnetic components of devices such as transformers and motors. This study presents a hysteresis model with a convenient analytical form, based on hypergeometric functions with one free parameter, built upon a class of parameterized curves. The aim of this work is to explore suitability of the presented model for describing major and minor loops, as well as to demonstrate improved agreement between experimental and modeled hysteresis loops. The procedure for generating first order reversal curves is also discussed. The added parameter, introduced into the model, controls the shape of the model curve, especially near saturation. It can be adjusted to provide better agreement between measured and model curves. The model parameters are nonlinearly dependent; therefore, they are determined in a nonlinear curve fitting procedure. The choice of the initial approximation and a suitable set of constraints for the optimization procedure are discussed. The inverse of the model function, required to generate first order reversal curves, cannot be obtained in analytical form. The procedure to calculate the inverse numerically is presented. Performance of the model is demonstrated and verified on experimental data obtained from measurements on construction steel sheets and grain-oriented electrical steel samples.

1. Introduction

Whenever a cyclic variation of magnetic field is applied to a ferromagnetic material, the so-called hysteresis loop is obtained. Usually, symmetric closed curves are considered, cf. Figure 1, yet it should be mentioned that in some applications it is crucial to predict the shapes of magnetization curves for arbitrary input signals. The practical importance of hysteresis loop stems from the fact that it allows one to determine the values of some significant properties of magnetic materials like remanence induction and coercive field strength (shown in the figure for the limiting (saturating) loop). Moreover, integration of loop area allows one to determine energy dissipated per cycle and per unit volume [1]. Therefore, the study of hysteresis phenomenon still attracts the attention of the engineering community.

The choice of material for a specific application may be facilitated using databases containing measurement results for different materials. However, the abundance of possible excitation scenarios makes the use of such libraries prohibitive from the practical point of view; therefore, much attention is paid to the development of mathematical descriptions of hysteresis loops. Advanced analyses and computations concerning e.g., distribution of magnetic fields and losses in cores of electric devices are often built upon these data.

A convenient model of hysteresis should be simple to use in calculations, described with a small set of parameters, and provide a good fit to the experimental data. Various mathematical models have been developed to describe hysteresis loops, such as the wide-spread Jiles–Atherton (JA) [2,3,4] and Preisach formalisms [5,6,7,8]. However, parameter determination for the JA model can be complicated, whereas the Preisach model can be computationally demanding. Significant effort is being invested in improving the accuracy of hysteresis models. Therefore of special interest are models based on convenient choice of analytical functions, whose parameters are easily determined. For example, in [9,10,11,12], a hysteresis model based on the hyperbolic tangent function is discussed, whereas, in [13,14], similar descriptions, involving the arctan function, are presented.

In this paper, we present an improved hysteresis model, based on hypergeometric functions, which defines a class of parameterized curves suitable for modeling of major and minor hysteresis loops, as well as first order reversal curves (FORCs). The method for generating FORCs, presented in this paper, is based on methods described in [15,16,17], with a modification based on numerical calculation of the inverse model function. Many authors have explored translation and scaling to generate FORCs and minor loops, to mention as representative examples Refs. [18,19,20,21]. Numerical procedures for parameter determination and FORC construction, specific to our model, are presented. The good accuracy of the model is verified on experimental data obtained from measurements on commercially available steel sheets and grain-oriented electrical steel.

2. Model Description

An essential factor in choosing a good analytical model is how accurately it describes the physical behavior under the applied magnetic field, when the domain structure changes due to domain wall motion, pinning, and wall rotation. Starting from a simple model for the magnetic moment of the atom and the quantum interpretation of spin, the energy of the magnetic dipole is expressed via probabilities; see, for example [9,22,23]. This approach leads to the model of an anhysteretic magnetization curve expressed with the hyperbolic tangent function

$$M=M_S\tanh \left(\frac{mH}{k_{B}T}\right)=M_S\frac{e^{2x}-1}{e^{2x}+1}=M_S\left(1-\frac{2}{e^{2x}+1}\right),$$

(1)

where $M_S$ is saturation magnetization, $m$ is magnetic moment, $k_B$ is the Boltzmann constant, and $T$ is temperature. $x=mH/k_{B}T$ is the reduced field strength expressed in dimensionless units.

Hysteresis loop branches can be constructed by shifting the anhysteretic curve (1) horizontally by the value of $h_c$, dependent on the coercive field. In this way, only irreversible magnetization occurring due to domain wall motion is included in the model. Reversible magnetization due to domain rotation is represented by the term $qH$, which exhibits an approximately linear dependence on the applied magnetic field.

It is worth mentioning that the well-known Landau–Lifshitz–Gilbert (LLG) equation can be used to describe the hysteresis phenomenon with a simplified micromagnetic model that can be combined with a magnetic circuit simulation [24], yet it requires a relatively long computation time. Reference [25] points out that some solutions to the LLG equation can be explicitly expressed with confluent hypergeometric functions, which are also included in the present model.

Analytical hysteresis models based on hyperbolic functions [10] usually depend on physical parameters such as magnetization at saturation, coercive force, and temperature. Hysteresis loops are thus mathematically expressed as $B=B\left(H\right)$ or $M=M\left(H\right)$. We refer the readers to the book [9] for further details.

Consider the confluent hypergeometric function [26,27,28]

$${}_1{F}_1\left(a;b;x\right)={\displaystyle \sum _{n=0}^{\infty }\frac{{\left(a\right)}_n}{{\left(b\right)}_n}\frac{x^n}{n!}}$$

(2)

where ${\left(a\right)}_n=a\left(a+1\right)\cdots \left(a+n-1\right)$ for $n>0$ and ${\left(a\right)}_0=1$, with ${\left(b\right)}_n$ defined in the same way.

For certain combinations of $a$, $b$, and $x$, the hypergeometric function reduces to some of the special functions. One such case is $e^x={}_1{F}_1\left(a;a;x\right)$ for any arbitrary $a$. Taking tanh from (1) and replacing $e^{2x}$ with

$${}_1{F}_1\left(a;1;2x\right)={\displaystyle \sum _{n=0}^{\infty }\frac{{\left(a\right)}_n}{{\left(n!\right)}^2}{\left(2x\right)}^n}$$

we get

$$f(a,x)=1-\frac{2}{1{+}_1{F}_1\left(a;1;2x\right)}$$

(3)

Let us now define

$$\begin{array}{l}{B}_{\mathrm{upper}}\left(H\right)=B_{s}f\left(a,\left(H+h_c\right){\tau }^{-1}\right)+qH,\\ B_{\mathrm{lower}}\left(H\right)=-B_{s}f\left(a,\left(-H+h_c\right){\tau }^{-1}\right)-qH,\end{array}$$

(4)

to describe the descending and ascending branches of hysteresis, respectively, where $h_c$ is coercive force, $B_s$ is saturation magnetic flux density, and $\tau$ is the scaling factor. For simplicity, we denote the model (4) with $F\left(a,H\right)$ unless a specific branch is addressed, and, for $a=1$, we use the shorthand notation $T\left(H\right)=F\left(1,H\right)$. The term $B_{s}f\left(a,\left(H+h_c\right)/\tau \right)$ in (4) is related to irreversible magnetization, and $qH$ to reversible magnetization. This form is well suited for symmetric major and minor loops. The constant $d$ can be added to (4) to account for vertical displacement occurring in asymmetric curves and FORCs, and also to compensate for adjustments needed to make the model curves coincide at loop tips.

The hysteresis model based on the hyperbolic tangent has a physical explanation, but it does not include all physical effects occurring in materials. By adjusting the parameter $a$ in (3), the model curve is adjusted to better represent measured data when the material approaches saturation under magnetization. As the distance between $a$ and 1 increases, the difference between $f\left(a,x\right)$ and the hyperbolic tangent becomes more prominent, cf. Figure 2. This property of (3) has been shown to lead to improved agreement between the model and measured data. However, the difference shown in Figure 2 must be accounted for in the parameter determination procedure.

3. Parameter Determination

Estimates of the model parameters $\tau ,B_s,h_c,q,d$ can be determined from the measured data; this method is usually supplemented with a numerical fitting procedure. The quality of the fit can be assessed by examining numerical indicators such as the relative error of the residual vector expressed via the Euclidean norm

$$R_V=\sqrt{{\displaystyle \sum _{i=1}^n{\left(F\left(a,h_i\right)-b_i\right)}^2}}/\sqrt{{\displaystyle \sum _{i=1}^n{b}_i{}^2}}$$

(5)

or graphically, by plotting the error expressed as a normalized difference between the model and experimental data. Here, $\left(h_i,b_i\right)$, $i=1,2,\dots ,n$ denotes an array of measured data. Curve fitting involves iterative solving of a non-linear least squares problem using methods such as Newton or Levenberg–Marquardt. As the process involves complex calculations, i.e., evaluating the Jacobian at each iteration, a good choice of initial values is important. Wolfram Mathematica software was used for numerical calculations.

3.1. Choice of the Initial Values

Initial values for the optimization procedure are determined from the physical properties of the curve in a manner similar to the one described in [10]. Model parameters can be obtained directly for $F\left(1,H\right)$ and then used as initial values in the fitting procedure with the additional parameter $a$ in $F\left(a,H\right)$. Depending on the shape of the curve and the choice of initial parameters, the optimization procedure may converge to a less than satisfactory stationary point, or not converge at all. If necessary, the initial values can be adjusted using interactive calculation possibilities of the Mathematica GUI by the Manipulate function. In this way, it is possible to visually track the effect of parameters on the model curve and determine constraints for the fitting procedure. From our experience, a reasonable constraint on the parameter $a$ prevents excess distortion of the model curve, for example:

$$0.5\le a\le 1.4$$

(6)

3.2. Coincidence Point Adjustment

The symmetrical $\tanh$ function alone cannot represent asymmetry in the magnetization and demagnetization parts of hysteresis. To resolve this, a displacement parameter $d$ is usually added to the model $T\left(H\right)$ to shift the upper and lower branches. An example is presented in Figure 3, with model parameters $h_c=382.4$, $\tau =272.9$, $B_s=0.554$, $q=2.17\times {10}^{-4}$, and $d=-0.117$. This simple approach cannot be applied to $F\left(a,H\right)$ due to the nonlinear dependence of parameters.

Let us observe the example in Figure 4, where the measured data is approximated with $T\left(H\right)$ for $h_c=518.1$, $\tau =1017.1$, $B_s=1.72$, $q=3.5\times {10}^{-5}$, and $d=-8.6\times {10}^{-5}$. Let us note that the model $T\left(H\right)+d$ provides a good fit for the loop in Figure 3. However, for the loop in Figure 4 the best fit value of $d$ is practically zero, i.e., $d$ has no favorable influence on the model curve, which visibly deviates from the measured data. This is a known limitation of the model $T\left(H\right)$.

On the other hand, the presented model (4) is based on an inherently asymmetrical function $f\left(a,H\right)$ which better describes asymmetric branches. Let us remark that parameter values for $F\left(a,H\right)$, $a\ne 1$ cannot simply be carried over from $T\left(H\right)$. Instead, they should be used as initial values in a new optimization procedure, which should be performed on all parameters simultaneously, including $a$.

Plotting $F\left(a,H\right)$ with the same set of parameters as in Figure 4 and some $a\ne 1$, the situation similar to the one in Figure 5 is obtained, where the upper and lower branches are apart by $\mathsf{\Delta }B$ and do not coincide at $\left(H_{tip},B_{tip}\right)$, Figure 5. To enforce coincidence, $B_s$ is adjusted by introducing the additional constraint

$$B_{upper}\left(H_{tip}\right)=B_{lower}\left(H_{tip}\right)=B_{tip}$$

(7)

which causes $B_s$ to deviate somewhat from its original physical interpretation. Furthermore, additional constraint is required to limit the deviation, for example

$$0.8b_{\max }\le B_s\le 1.2b_{\max }$$

(8)

Although at the first sight it may seem that this situation could have been resolved by simply scaling the saturation value or by shifting the branches by $\mathsf{\Delta }B$, the effects of the parameters $a$, $B_s$ and $d$ are nonlinearly dependent and as such must be treated simultaneously. Thus, the constraints (6), (7) and (8) must be introduced into the fitting problem.

With the introduced constraints, a new fitting problem with one additional parameter, $a$, is formed and solved. The resulting model curve for $a=0.65$, $h_c=540.1$, $\tau =489.8$, $B_s=1.53$, $q=8.7\times {10}^{-5}$, and $d=0.63$ shows good conformance to experimental data (Figure 6).

The comparison of relative errors and mean relative errors between models is shown in Figure 7 for the upper branch of the loop; the results are symmetrical in $H$ for the lower branch. The model $F\left(a,H\right)$ yields a lower maximum relative error of 4.5% versus 13% for $T\left(H\right)$ and a lower mean error of 2% versus 4.5%.

4. Measurements

Experimental hysteresis data was obtained in a laboratory, at room temperature, using ring specimens equipped with primary and secondary windings. A principal scheme of the measurement setup is shown in Figure 8. The value of magnetic field strength H is proportional to current $i_1$ and the voltage drop $u_1$ on resistor $R_1$. The value of magnetic flux density B is proportional to integrated voltage induced in the secondary winding of the sample. The values of voltage drop $u_1$ and output voltage of the integrator are simultaneously sampled on the data acquisition equipment (DAQ) and stored. The setup excites the sample with magnetic field strength H, generated by the programmable generator of the control signal. Measurements may be carried out using sinusoidal excitation for frequencies from 10 Hz to 1 kHz.

Measurements were conducted using the measurement setup shown in Figure 9. Excitation for the primary coil was supplied from the programmable GF-1 function generator. A Tektronix TDS5032 digital oscilloscope was used for simultaneous data acquisition on two channels. Voltage on a precision 1 Ω shunt of negligible inductance, which was connected in series with the primary coil of the toroidal sample, was sampled on the first channel. Output from the electronic integrator was sampled on the second channel. Measurements were performed for excitation frequencies from 10Hz, 50Hz up to 450Hz and for current intensities from tens to hundreds of milliamperes. For the purpose of modeling, we have used the lowest possible frequency in order to avoid the distorting effect of eddy currents induced in the conductive core material on the shape of hysteresis loop.

Sample 1 was a specimen of commercially available 0.4mm steel sheets used for manufacture of device cases, whereas Sample 2 was a specimen of grain oriented electrical steel (M103-27p), cf. Figure 10.

5. Numerical Results

The accuracy and practical suitability of the presented model was tested on experimental data obtained from laboratory measurements presented in the previous section. Hysteresis data was obtained by measuring primary current and integrated output voltage from the secondary coil of the test core, with different excitation currents at a low excitation frequency. Calculations were performed in Mathematica, using the NonlinearModelFit function.

5.1. Major Loops

The major loop for Sample 1 is shown in Figure 4 and Figure 6, modeled with $T\left(H\right)$ and $F\left(a,H\right)$ respectively. Model errors are compared in Figure 7. This example was used to discuss the parameter determination procedure in Section 3, and the advantages of the presented model are also discussed there.

The major loop for Sample 2 is shown in Figure 11 and Figure 12, modeled with $T\left(H\right)$ and $F\left(a,H\right)$ respectively. Parameter values are $h_c=68.75$, $\tau =5.66$, $B_s=1.44$, $q=2.5\times {10}^{-4}$, $d\approx 0$ for $T\left(H\right)$ and $a=0.52$, $h_c=70.78$, $\tau =6.29$, $B_s=1.53$, $q=1.86\times {10}^{-4}$, $d=0.055$ for $F\left(a,H\right)$. It can be seen that the presented model follows the measured data exceptionally well during the magnetization stage of the hysteretic cycle. The comparison of relative errors and mean relative errors between the two models (Figure 11 and Figure 12) is shown in Figure 13. With the mean error of 10%, $F\left(a,H\right)$ is clearly better than $T\left(H\right)$ yielding 21%.

5.2. First Order Reversal Curves

First order reversal curves can be calculated in various ways from the model for the major loop. We present a technique based on the method established in [17], which in turn is based on earlier papers [15] and [16]. In [17], hysteresis loops are expressed in terms of magnetizing current versus flux linkage, $i_m\left(\lambda \right)$ in the $i\lambda$ plane. A first order reversal curve is defined as $i_m=f\left(\lambda \right)=F\left(\lambda +x\right)$, where $i_m=F\left(\lambda \right)$ is the equation of ascending branch of major hysteresis loop and $x=X\left(\lambda \right)$ is a function describing the displacement of the major branch along the $\lambda$ axis. To apply the considered method to the model (4), formulas must be applicable in the HB plane.

A first order reversal curve is generated between an arbitrary starting point on the descending major branch and the positive loop tip (Figure 14). The curve is obtained by displacing the major branch by an amount expressed with the function

$$X\left(b\right)=\frac{1}{2}{x}_1\left(1-{\left(\frac{b-b_1}{B_{tip}-b_1}\right)}^p+{\left(\frac{B_{tip}-b}{B_{tip}-b_1}\right)}^w\right).$$

(9)

where $p=1-\alpha k$, $w=1+\beta k$, $k=\left(B_{tip}+b\right)\left(B_{tip}+b_1\right)B_{tip}^{-2}$. Parameters $\alpha$ and $\beta$ depend on ferromagnetic material and can be adjusted to influence the slope of the magnetic path. Here $\alpha =0.1$ and $\beta =0.3$.

For a chosen starting point 1 with coordinates $\left(h_1,b_1\right)$, the value of the parameter $x_1$ is calculated as the difference $x_1=\left|B_{lower}\left(h_1\right)-B_{upper}\left(h_1\right)\right|$. It represents the initial displacement between points 1′ and 1.

An arbitrary point P on the reversal curve is calculated as

$$P=\left(B_{lower}^{-1}\left(b-X\left(b\right)\right),b\right)$$

for $b\in \left[b_1,B_{tip}\right]$. Since the ascending branch in the model (4) is expressed with $B_{lower}\left(H\right)$, the inverse function $H=B_{lower}^{-1}\left(B\right)$ is needed for computation. Analytical form of the inverse hypergeometric function is not feasible, therefore one must resort to numerical computation.

The value of the inverse function is computed by numerically solving the equation $B_{lower}\left(h\right)=b_0$ for any given value $b_0=b-X\left(b\right)$. The initial value $h_0$ for the root-finding method is chosen as follows. The hysteresis branch is approximated with several linear segments, constructed by a linear fit. The value $h_0$ is obtained from the inverse of the closest linear segment (Figure 15). The root-finding method then calculates the numerical solution $h_n$ in $n$ iterations, which is taken as the value of the inverse $B_{lower}^{-1}$. Thanks to the shape of hysteresis branches, convergence conditions for the iterative method are satisfied for all values of $b_0$. In our experiments, using five linear segments (Figure 15), the iterative method took an average of 3.6 iterations to calculate and plot $B_{lower}^{-1}$.

5.3. Minor Loops

One of the important properties of hysteresis models is the ability of the model to represent minor loops. There are various approaches to minor loop construction, e.g., by shifting and scaling major loop branches and adjusting the loop tips. Several minor loops for Sample 1, approximated by $T\left(H\right)$ and $F\left(a,H\right)$ are shown in Figure 16 and Figure 17, respectively. Comparing relative errors ($R_V$) as the minor loops approach the major loop, advantage of the model (4) becomes more pronounced (Figure 17).

6. Conclusions

By employing model curves based on a family of parameterized hypergeometric functions, a more accurate hysteresis model with one free parameter is constructed. The presented model has a convenient analytical form that makes it suitable for practical applications. The conditions for parameter determination using a nonlinear fitting procedure are discussed, as well as the applied modification of the method for generating first order reversal curves. Using experimental data, it is shown that the model can be used with good accuracy to represent major and minor loops, as well as simulate first order reversal curves.

The ease of use and the flexibility of simple analytical T(x) model is retained in the case of the extended description, based on hypergeometric functions with a single free parameter. The generalized model is able to describe both major and minor loops of representative soft magnetic materials accurately. This fact is important from the practical point of view, since the examined description is much simpler in implementation than the Jiles-Atherton or Preisach models. The application of hypergeometric functions allows one to consider different classes of magnetic materials. Since the hypergeometric function may be reduced to hyperbolic tangent in the limiting case and the description based on hyperbolic tangent i.e., the T(x) model is particularly suited for materials with strong uni-axial anisotropy, we believe that our generalized model may represent soft magnetic materials with differing anisotropy level. The examination of physical meaning of hypergeometric series shall be the subject of forthcoming research.