**4. Some applications of the hysteresis modeling in electrical engineering**

In this section we present some applications of hysteresis modeling in electrical engineering. The first two are motivated by the hysteresis observed when modeling the BH curve of a ferromagnetic material. The third one has to do with the hysteresis present in the State of Charge (SoC)-Open Circuit Voltage (OCV) mapping of Li-ion batteries.

#### **4.1 Hysteresis power losses computation**

For the performance analysis of electrical machines, an important factor to be considered is the power losses that occur in the ferromagnetic materials that make up the core of the machine. These losses, usually known as iron losses, can generally be divided into eddy current, hysteresis and excess losses. Eddy current losses are resistive losses due to the currents induced in the magnetic material by the time varying magnetic induction; its magnitude strongly depend on the size of the continuous conductive regions. That is why the core of the machine is usually assembled from thin steel sheets, insulated from each other by a non-conductive coating on the surface.

The hysteresis and excess losses are essentially related to the microscopic properties of ferromagnetic materials. These consist of small magnetic domains, which tend to align along the external (time – varying) applied magnetic field. As a consequence of the molecular friction produced in this continuous movement, heat losses, generally referred to as magnetic hysteresis losses, are produced. Given a material point *x*, the density of hysteresis losses in the time interval ½ � *t*1, *t*<sup>2</sup> can be computed from the integral:

$$\int\_{t\_1}^{t\_2} \frac{\partial \mathbf{B}}{\partial t}(\mathbf{x}, t) \cdot \mathbf{H}(\mathbf{x}, t) \, dt. \tag{18}$$

Nevertheless, these losses are especially difficult to compute, since for ferromagnetic materials, the magnetic induction in each point depends on the intensity of the present magnetic field to which it is exposed, and also on previous exposures to magnetic field intensity of each volume element. This causes differences in the magnetization curve under increasing and decreasing fields; therefore, hysteresis loops arise as the one illustrated in **Figure 9** (left).

In the literature there are several simplified analytical expressions that are used to approximate the different components of the losses, such as those proposed by Bertotti [14], which are among the most widely used. However, the assumptions under which it is possible to apply these formulas are not met in most practical situations. In this context, numerical simulation is a viable option to overcome these limitations.

As an example, we consider an application consisting of the computation of hysteresis and eddy current losses in a laminated medium with toroidal geometry, as sketched in **Figure 9**. It consists of *N* circular sheets of rectangular section of thickness *d* surrounded by a coil. We denote by *R*<sup>1</sup> and *R*<sup>2</sup> the internal and external radius of the core, respectively. We will assume that the coil is infinitely thin so that it can be modeled as a surface current density (A/m).

As shown in [24], in this situation, it is possible to reduce the computational domain to the meridian section of one single sheet for which it is necessary to derive appropriate boundary conditions. More precisely, the axisymmetric transient eddy current problem to be solved reads:

*Find H r*ð Þ , *z*, *t and B r*ð Þ , *z*, *t such that*

$$\frac{\partial B}{\partial t} - \frac{\partial}{\partial r} \left( \frac{1}{\sigma r} \frac{\partial (rH)}{\partial r} \right) - \frac{\partial}{\partial \mathbf{z}} \left( \frac{1}{\sigma} \frac{\partial H}{\partial \mathbf{z}} \right) = \mathbf{0} \quad \text{in } \Omega \times (0, T), \tag{19}$$

$$B(r,z,t) = \mu\_0(H(r,z,t) + [\mathcal{F}\_\mathcal{S}(H,\xi)](r,z,t)) \quad \text{in } \Omega \times (0,T), \tag{20}$$

$$H(r, z, t) = h(r, z, t) \quad \text{on} \quad \Gamma \times (0, T), \tag{21}$$

$$B(r, t, \mathbf{0}) = B\_0(r, z) \quad \text{in} \quad \Omega,\tag{22}$$

where Ω ¼ ½ �� *R*1, *R*<sup>2</sup> ½ � 0, *d* , Γ its boundary, and *σ*ð Þ *r*, *z*, *t* , *h r*ð Þ , *z*, *t* , *ξ*ð Þ *r*, *z* and *B*0ð Þ *r*, *z* are given functions. We have employed the usual notation: *H* is the magnetic field, *B* the magnetic induction and *σ* the electrical conductivity. Moreover, it can be shown that (see [24] for details):

$$h(r, z, t) \coloneqq \frac{n\_\epsilon I(t)}{2\pi r},\tag{23}$$

**Figure 9.**

*Toroidal laminated media (left) and its meridian section (center); hysteresis loop measured and approached (right).*

#### *Preisach Hysteresis Model – Some Applications in Electrical Engineering DOI: http://dx.doi.org/10.5772/intechopen.99590*

where *I t*ð Þ is the current intensity flowing through the coil at time *t* and *ne* denotes the number of winding turns. In order to model the B-H relation, the classical Preisach hysteresis operator was employed. The Everett function was approximated from experimental data. **Figure 9** (right) shows a comparison between the B-H hysteresis curves measured experimentally and those computed from the Everett function thus assessing the validity of the approximation.

The energy losses can be evaluated by computing the magnetic field and magnetic induction solutions to (19)–(22) with appropriate numerical techniques, such as Finite-Difference Time-Domain (FDTD) [33, 34] or Finite-Element method (FEM) (see, for instance [35–42] where scalar and vector hysteresis models are considered, respectively). In particular, we have applied the latter and a fixed-point iteration (see [24], Section 3) and also [43, 44] for different fixed point iteratios applied to hysteresis models).

In order to compute the losses, we consider the energy balance in the axisymmetric setting, which reads [24]:

$$\underbrace{\int\_{\Omega} \left( \int\_{t\_1}^{t\_2} \frac{\partial B}{\partial t} H dt \right) r \, dr dz}\_{\mathcal{L}\_H} + \underbrace{\int\_{t\_1}^{t\_2} \left( \int\_{\Omega} \frac{|\mathbf{curl}(H\mathbf{e}\_\theta)|^2}{\sigma} r \, dr dz \right) dt}\_{\mathcal{L}\_E} = \underbrace{\int\_{t\_1}^{t\_2} \left( (rH)\_{|\Gamma} \frac{\partial}{\partial t} \left( \int\_{\Omega} B \, dr dz \right) \right) dt}\_{\mathcal{L}\_T} \tag{24}$$

In this expression, the terms L*E*, L*<sup>H</sup>* and L*<sup>T</sup>* represent the eddy current, the hysteresis and the total losses, respectively. From the numerical point of view, these losses have been computed (in J*=*m<sup>3</sup> ð Þ) as

$$\mathcal{L}\_{\rm E}^h \coloneqq \frac{1}{\pi d \left(R\_2^2 - R\_1^2\right)} \left\{ \int\_0^T \left(2\pi \int\_{\Omega} \frac{|\operatorname{curl}(H\_h \mathbf{e}\_\theta)|^2}{\sigma} r dr dz\right) dt \right\},\tag{25}$$

$$\mathcal{L}\_H^h \coloneqq \frac{1}{\pi d \left(R\_2^2 - R\_1^2\right)} \left\{ 2\pi \int\_{\Omega} \left(\int\_0^T \frac{\partial B\_h}{\partial t} H\_h dt\right) r dr dz \right\},\tag{26}$$


*a Total losses computed experimentally.*

*b Relative error between experimental losses and* L*<sup>h</sup> <sup>E</sup>* <sup>þ</sup> <sup>L</sup>*<sup>h</sup> H. <sup>c</sup>*

*Relative error between experimental losses and* L*<sup>h</sup> T.*

#### **Table 1.**

*Total losses. Numerical versus experimental results (J/*m<sup>3</sup>*).*

$$\mathcal{L}\_T^h \coloneqq \frac{1}{\pi d \left(R\_2^2 - R\_1^2\right)} \left\{ 2\pi \int\_0^T R\_2 H\_h \left| \begin{matrix} \left(\frac{d}{dt}\right)\_\Omega B\_h dr dz \\ \end{matrix} \right|\_\Omega B\_h dr dz \right\}.\tag{27}$$

The measurements reported in **Table 1** were performed on an Epstein frame considering a material sheet of width 30 mm and thickness 0.5 mm. The sheet is subjected to sinusoidal flux excitation with peak induction levels *Bm* equal to 0*:*5, 0*:*9 and 1*:*4 T and frequencies *f* equal to 25 and 150 Hz. For each of these peak levels and frequencies, the physical measurements were the total electromagnetic losses per cycle and unit volume and the magnetic field on the boundary of the sheet. The following geometrical data were considered to simulate the experimental setting with our axisymmetric model: *R*<sup>1</sup> ¼ 100 m, *R*<sup>2</sup> ¼ 100*:*3 m, *d* ¼ 0*:*0005 m, *σ* ¼ 4064777 (Ohm/m)�<sup>1</sup> . **Table 1** summarizes the results obtained for different frequencies and magnetic induction peak levels. Let us recall that equality L*<sup>E</sup>* þ L*<sup>H</sup>* ¼ L*<sup>T</sup>* holds at the continuous level and thus, the difference among the losses L*<sup>h</sup> <sup>E</sup>* <sup>þ</sup> <sup>L</sup>*<sup>h</sup> <sup>H</sup>* and L*h <sup>T</sup>* (columns 6 and 7 in **Table 1**) is due to numerical approximation of the axisymmetric model.

This table shows that the computed losses are close to the losses measured in the experiments. The largest differences are obtained for *Bm* ¼ 0*:*5 T and are probably due to the fact that our Everett function is unable to generate accurate B-H cycles for "small" values of *B* (see **Figure 9**, right), making the numerical approximation of the hysteresis losses more inaccurate. These discrepancies can be explained by the *congruency properties* of the classical Preisach model which states that B-H cycles between two fixed extreme magnetic fields H are independent of the induction level B. To avoid these discrepancies, a Moving Preisach hysteresis model can be applied. In this model the congruency property of minor loops is relaxed and it is expected to reproduce lower symmetric hysteresis loops (see [12]).

#### **4.2 Hysteresis modeling in magnetic inspection particle processes**

The technique known as Magnetic Particle Inspection (usually MPI for its acronym in English) is a non-destructive method for detecting flaws located on or near the surface of ferromagnetic parts. It exploits the fact that when a ferromagnetic sample is exposed to the influence of a magnetic field, the induced magnetic flux density accumulates inside the material. Then, if there is a crack, the magnetic field will be distorted, causing local magnetic leakage around the defect. Therefore, if fluorescent magnetic particles are sprayed on the magnetized sheet, they will concentrate on the cracks and produce deposits that are easy to identify under ultraviolet light. The purpose of MPI is to find faults in any direction. Because the breakings are easier to detect when orientated perpendicular to the specimen's magnetic field, it is common to apply the

*Preisach Hysteresis Model – Some Applications in Electrical Engineering DOI: http://dx.doi.org/10.5772/intechopen.99590*

magnetization in two orthogonal orientations, such as circular and longitudinal orientation (see **Figure 10**, left). Most times, after inspection, the workpieces must be demagnetized since residual magnetism could interfere with later processing.

In [25], the authors introduced two models for the numerical simulation of the entire magnetization and demagnetization procedures entailed in an MPI test. In particular, the remanent flux density in specimens with cylindrical symmetry is computed. The method is based on scalar models that include magnetic hysteresis and is carried out in three steps: long-term magnetization, circular magnetization and final demagnetization. To achieve circular magnetization (and demagnetization), the workpiece is clamped between two electrical contacts between which an electric current is circulated (see **Figure 10**). For modeling purposes, all fields are assumed to be *θ*independent. It is also assumed that the current density has the form *J x*ð Þ� , *t Jz*ð Þ *r*, *z*, *t ez*, that is, that it traverses the piece along its axial direction. Consequently, the magnetic field stands *H x*ð Þ� , *t Hθ*ð Þ *r*, *z*, *t eθ*. As proven in [25], the problem can be spatially reduced to the meridional section Ω of the workpiece, and reads

$$\frac{\partial \mathcal{F}\_{\mathcal{S}}(H\_{\theta}, \xi)}{\partial t} \mathbf{e}\_{\theta} + \mathbf{curl} \Big( \frac{1}{\sigma} \mathbf{curl} (H\_{\theta} \mathbf{e}\_{\theta}) \Big) = \mathbf{0} \quad \text{in} \ \Omega \times [t\_{0}, T], \tag{28}$$

$$H\_{\theta}(\mathbf{0}, z, t) = \mathbf{0} \quad \text{on} \quad (\mathbf{0}, L) \times [t\_0, T], \tag{29}$$

$$H\_{\theta}(R\_{S}(\mathbf{z}), \mathbf{z}, t) = \frac{I(t)}{2\pi R\_{S}(\mathbf{z})} \quad \text{on} \quad (\mathbf{0}, L) \times [t\_{0}, T], \tag{30}$$

$$\frac{\partial H\_{\theta}}{\partial \mathbf{z}}(r, \mathbf{z}, t) = \mathbf{0} \quad \text{on} \quad (\Gamma\_1 \cup \Gamma\_2) \times [t\_0, T], \tag{31}$$

where Γ1, Γ<sup>2</sup> denote the boundaries clamped to the electrical contacts, *L* is the workpiece length in the axial direction, S*<sup>z</sup>* any cross-section transversal to this direction and *R*Sð Þ*z* its radius. As in previous sections, F*<sup>S</sup>* is the classical Preisach hysteresis operator and *ξ* the initial magnetization state.

For the longitudinal magnetization and demagnetization processes, the piece is placed inside a conducting coil carrying an alternating current in the azimuthal direction (see **Figure 10**). To avoid the use of a vector hysteresis law, an infinitely thin conducting surface Γ<sup>S</sup> of radius *R*<sup>S</sup> carrying a surface current density *J*Sð Þ� *x*, *t J*Sð Þ *r*, *t e<sup>θ</sup>* is employed to approximate the coil. Let Ω*<sup>c</sup>* be the part to be inspected, *Rc* representing its radius; moreover, let us denote Ω<sup>0</sup> the air surrounding it which is assumed to be artificially bounded in the *r*-direction, *R*<sup>∞</sup> being its outer radius. Then problem reduces to (see [25])

$$\sigma \frac{\partial A\_{\theta}}{\partial t} - \frac{\partial}{\partial r} \left( \mathcal{F}\_{\mathbb{S}}^{-1} \left( \frac{1}{r} \frac{\partial}{\partial r} (r A\_{\theta}), \xi \right) \right) = 0 \quad \text{in} \quad (\mathbf{0}, \mathbb{R}\_{\epsilon}) \times [\mathfrak{t}\_{0}, T], \tag{32}$$

$$\frac{\partial}{\partial r}\left(\frac{1}{\mu\_0}\frac{1}{r}\frac{\partial}{\partial r}(rA\_\theta)\right) = 0 \quad \text{in} \quad [(R\_\circ, R\_\circ) \cup (R\_\circ, R\_\circ)] \times [t\_0, T],\tag{33}$$

$$\left[\frac{1}{\mu\_0} \frac{1}{r} \frac{\partial}{\partial r} (r A\_\theta)\right]\_{r=R\_\\$} = f\_\\$ \quad \text{on} \quad \{R\_\\$\} \times [t\_0, T], \tag{34}$$

$$A\_{\theta}(\mathbf{0}, t) = \mathbf{0} \qquad \text{and} \qquad \frac{1}{r} \frac{\partial}{\partial r} (r A\_{\theta})(R\_{\circ \circ}, t) = \mathbf{0}, \forall t \in [t\_0, T], \tag{35}$$

**Figure 11.**

*Left – Major hysteresis loop and (blue) and initial magnetization curve (red). Center – Demagnetizing source. Right – Circular demagnetization. Remanent flux density at z* ¼ 0*:*2 *m vs. radius. Number of cycles comparison.*

where *A<sup>θ</sup>* denotes the azimuthal of component magnetic vector potential. In (34) ½ �� *<sup>r</sup>*¼*R*<sup>S</sup> represents the jump across the surface *<sup>r</sup>* <sup>¼</sup> *<sup>R</sup>*S. Notice that the source is the surface current density *J*<sup>S</sup> which is equal to the jump of *H* � *n* at *r* ¼ *R*S.

It is worth noting that when solving (32)–(35), the effective computation of the inverse of the hysteresis operator is avoided by using an iterative algorithm of a fixed point type, based on the properties of the Yosida regularization for maximal monotone operators (see [45]). Moreover, the fields *H* and *B* have only one non-null component in the **e***<sup>θ</sup>* direction in the circular case, and in the **e***<sup>z</sup>* direction in the longitudinal one. This allows using a scalar constitutive law for the nonlinear material, in both the circular and the longitudinal formulations.

For the sake of brevity, only the numerical results corresponding to the circular case are included. The ferromagnetic piece is characterized by the initial magnetization curve and the major hysteresis loop, as shown in **Figure 11**, which was obtained by using an artificial weight function *p*. The magnetization stage is performed by using a harmonic surface current density on the cylindrical surface Γs. When the source is turned off, the workpiece retains a residual field or remanence. In order to remove this remanent field, a surface current density on Γ<sup>s</sup> is considered with amplitude equal to the one used for magnetizing the piece and continuously reversing while gradually decreasing to zero (see **Figure 11**). The numerical results give important information about the source needed to successfully demagnetize the piece. In particular, the efficiency of the demagnetization process strongly depends on the source and the number of cycles. In particular, the greater the number of source cycles, the better the demagnetization.

#### **4.3 Hysteresis model applied to battery modeling**

It is well-known that hysteresis has a significant impact on the ability to monitor battery performance since it affects the voltage during charging and discharging cycles for different battery chemistries. Therefore, it must be considered for diagnosis algorithms that use the so-called open circuit voltage (OCV) as a measure of the state of charge (SoC) of the battery. SoC indicates the residual charge of the battery with respect to its maximum nominal and it is usually expressed as a percentage of the battery capacity. Since the SoC provides a measure of the available energy of the battery, as well as of its instantaneous power capacity, an accurate estimate is important to monitor battery performance. Unfortunately, it cannot be directly measured but must be estimated from other battery quantities (see, for instance [46–50], to mention some recent works). However, the non-linear electrochemical reactions involved in the

*Preisach Hysteresis Model – Some Applications in Electrical Engineering DOI: http://dx.doi.org/10.5772/intechopen.99590*

**Figure 12.**

*Measurements along the major loop and the approximation given by the Preisach model (left). Experimental inner loops and approximation of the Preisach model for different SoC values (right).*

battery system, the effects of temperature, aging and, specially, hysteresis effects make this task especially difficult.

For a battery cell, hysteresis means that the battery cell reaches different OCV values at the same SoC value and temperature between charge and discharge, depending on the previous charge–discharge history. In particular, from experiments it is shown that, after discharging, the OCV always relaxes below the OCV after charging for the same SoC (see **Figure 12**, left). As a consequence, the SoC-OCV relation is a bundle of curves enclosed in a major loop. The loop consists of two SoC-OCV curves obtained by fully charging and discharging the battery, first to minimum voltage and then to maximum voltage while recording the battery voltage and accumulated ampere hours, for a complete battery cycle. Minor loops lying inside the major SoC-OCV loop can be obtained by conducting similar experiments but with partial cycles (see **Figure 12**, right). A simple approach is to compute the average of the major loop and to relate the battery's OCV with its SoC through this average single-valued curve (see **Figure 12**, left), but usually this approximation suffers from significant errors in the real SoC estimation.

Based on the work [51], in this section we consider a Preisach-type hysteresis model to approximate the SoC-OCV hysteresis loop (see [52–54] for different approaches to handle hysteresis in batteries). The battery dynamics is modeled by using an equivalent circuit model (ECM) as extensively used in the literature [55]. In particular, we consider a general *j*-th order circuit (see **Figure 13** (left) for the case *J* = 2). It consists of the so-called equivalent series resistance and a number of *J* resistorcapacitor *Rj*,*Cj* elements in series. The latter allow us the modeling of the diffusion voltages and are a good approximation of the so-called Warburg impedance. The notations employed in the sequel are the usual ones: *Q* is the maximum charge (or capacity) of the battery, *i t*ð Þ is the current intensity, *η*ð Þ*t* the Coulombic efficiency,

**Figure 13.** *Sketch of the 2RC model and parameters values.*

**Figure 14.**

*Percentage error in the voltage approximation (left). Measured and computed voltage with ECM and Preisach model (right).*

*V t*ð Þ the voltage, *z t*ð Þ the SoC, *Rj* the *j*�th polarization resistance, *Cj* the *j*�th polarization capacitance, and *vRj* ð Þ*t* is the voltage between the ends of the the *j*-th resistor *Rj*. Thus, the problem to solve reads:

*Given functions OCV, i t* d ð Þ *and η*ð Þ*t , constant parameters Q*, *R*0, *Rj and Cj, and initial conditions z*<sup>0</sup> *and vRj*,0 *, j* ¼ 1, … , *J, find functions z t*ð Þ*,V t*ð Þ *and vRj* ð Þ*t , satisfying,*

$$\frac{dz}{dt}(t) = -\frac{1}{Q}\eta(t)i(t) \tag{36}$$

$$\frac{d\boldsymbol{v}\_{R\_{\boldsymbol{\gamma}}}}{dt}(t) + \frac{1}{R\_{\boldsymbol{\gamma}}\mathbf{C}\_{\boldsymbol{\gamma}}}\boldsymbol{v}\_{R\_{\boldsymbol{\gamma}}}(t) = \frac{1}{\mathbf{C}\_{\boldsymbol{\gamma}}}\boldsymbol{i}(t), \quad \boldsymbol{j} = \mathbf{1}, \ldots, \boldsymbol{J} \tag{37}$$

$$\boldsymbol{\nu}\_{\mathbb{R}\_{\rangle}}(\mathbf{0}) = \boldsymbol{\nu}\_{\mathbb{R}\_{\geq 0}} \quad j = \mathbf{1}, \ldots, J \tag{38}$$

$$\mathbf{z(0)} = \mathbf{z\_0} \tag{39}$$

$$V(t) = \widehat{OCV}(z(t)) - \sum\_{j=1}^{J} v\_{R\_j}(t) - R\_0 i(t) \tag{40}$$

*where function OCV gives the OCV from the state of charge z* d .

From a macroscopic standpoint, the SoC-OCV hysteresis relation can be considered rate-independent; this implies that the OCV depends on the SoC history and not on the velocity (battery current rate) of SoC. Our aim is now to apply a hysteresis model to a lithium-iron-phosphate (LFP) cell. It is known that the classical Preisach operator cannot be directly applied to model the SoC-OCV hysteresis loop. Thus, the hysteretic behavior is modeled with a modified Preisach operator in which the SoC is assumed to be the independent variable. In particular we consider

$$\widehat{OCV}(z(t)) = \int\_{S\_x^+(t)} p(\rho)d\rho + \text{OCV}\_{\text{min}} \tag{41}$$

where *S*<sup>þ</sup> *<sup>z</sup>* is defined in (9) and *OCVmin* is the minimum OCV value. As usual, the computation of the hysteretic term in (41) is done by means of the Everett function instead of the Preisach hysteresis function *p*. This function is identified from experimental data considering only the charging SoC-OCV curves of the battery, which is less time consuming than other approaches proposed in the literature. The SoC-OCV curves computed with the Preisach model (41) are compared with experimental data. The approximations are reported in **Figure 12**, where the experimental data and values simulated with the Preisach model are plotted. From this figure, it can be seen that the model approximates the measurements of both the major loop and the inner loops.

*Preisach Hysteresis Model – Some Applications in Electrical Engineering DOI: http://dx.doi.org/10.5772/intechopen.99590*

Next, we consider a realistic vehicle driving profile to quantify the accuracy of the ECM with hysteresis. We have solved (36)–(39) and compared the computed voltage (40) with the experimental values. The values of OCV are obtained with the average open circuit voltage curve (see **Figure 12**, left) and with the Preisach hysteresis model (cf. (41)). **Figure 14** shows the results obtained when a second-order circuit model (2RC) is considered. More precisely, this figure shows, on the left, the measured voltage and on the right the error between the experimental voltage and the voltage computed with the average curve (in red) and the Preisach model (in blue), respectively. The result has been excellent when the Preisach model is considered, having a relative error of 0*:*4%.
