Harmonic and DC Bias Hysteresis Characteristics Simulation Based on an Improved Preisach Model

Transformers, reactors and other electrical equipment often work under harmonics and DC-bias working conditions. It is necessary to quickly and accurately simulate the hysteresis characteristics of soft magnetic materials under various excitation conditions in order to achieve accurate calculations of core loss and the optimal design of electrical equipment. Based on Preisach hysteresis model, a parameter identification method for asymmetric hysteresis loop simulation is designed and applied to the simulation of hysteresis characteristics under bias conditions of oriented silicon steel sheets. In this paper, the limiting hysteresis loops of oriented silicon steel sheets are obtained through experiments under different working conditions. The first-order reversal curves(FORCs) with asymmetric characteristics is generated numerically, and then the Everett function is established under different DC bias conditions. The hysteresis characteristics of the oriented silicon steel sheets under harmonics and DC bias are simulated by improving FORCs identification method of the Preisach model. By comparing the results of simulation and experiment, the effectiveness of the proposed method is verified, so as to provide an important reference for material production and application.


Introduction
The concept of the Preisach hysteresis model was proposed by the German physicist F. Preisach in 1935. After improved by several researchers, the classical Preisach model with mathematical expressions was later widely used [1][2][3]. The Preisach model can be classified into continuous models and discrete models. The former needs to be realized by identifying all parameters in the distribution function which identification process is complex and computationally expensive [4][5][6][7][8], the latter requires the identification of the Everett function matrix by first order reversal curves (FORCs) [9]. However, in the actual measurement process, due to the asymmetric excitation, the DC bias component is brought, and the measurement result has a large deviation. In order to avoid the large error caused by the measurement, Naidu converted the Preisach distribution function into a form that only requires the limiting hysteresis loop to obtain the FORCs and the Everett function [10]. The classical Preisach model takes the magnetic field strength as the independent variable and the magnetic flux density as the dependent variable, it is not easy to apply in finite element analysis. Dlala deduces the specific method of determining the FORCs according to the limiting hysteresis loop and the reversal points for the inverse Preisach model with the magnetic flux density as the input and the magnetic field strength as the output. The Everett function of magnetic flux density is obtained by interpolating the FORCs [11].
Hans Hauser proposed the energetic hysteresis model based on energy balance and magnetic domain statistics theory in 1994 [12][13][14]. The model employs an energy equation to describe the hysteresis phenomenon in ferromagnetic materials, in which the magnetic energy of the crystal is decomposed into an irreversible and a reversible part, and uses a probability function to represent the interaction of domain wall motion and stray fields at the pinning center. The parameters of magnetic characteristics are related to spontaneous magnetization, magnetocrystalline anisotropy, magnetostriction, and microstructure, and can account for the effects of stress, temperature, and magnetization direction on hysteresis. In addition, compared with the J-A model, the energetic model also has the advantages of faster solution rate and fewer parameters.
In addition to the above-mentioned commonly used scalar hysteresis models, there are also some vector hysteresis models, such as the Enokizono and Soda (E&S) model [15]. The vector hysteresis models can be constructed by superimposing the scalar models in different directions, such as the vector Preisach model [16], the vector J-A model [17].
Silicon steels have been widely made as the core of transformers and motors. However, with the development of power systems, the working environment of power equipmentalso becomes more complex, especially at high-frequency, high temperature, stress, nonsinusoidal excitation. The commonly-used hysteresis models and loss calculation methods are difficult to simulate magnetic properties under complex conditions. Therefore, it is necessary to establish accurate and efficient models suitable for various complex excitation conditions to quickly simulate the magnetic properties of materials.

Dynamic Preisach Hysteresis Model
The basic idea of the classical Preisach model is that ferromagnetic materials are composed of a large number of magnetic dipoles with hysteresis properties, and in the case of an output value u(t) at a given moment, the hysteresis characteristics of ferromagnetic materials produce the sum of the hysteresis characteristics of all magnetic dipoles. The output value of the magnetic dipole γ αβ u(t) is assumed to be either +1 or −1 values, µ(α,β) is a distribution function of α, β, also known as the Preisach function ( Figure 1). The classical Prcisach model can be represented by the following equation: (1) magnetic domain statistics theory in 1994 [12][13][14]. The model employs an energy equation to describe the hysteresis phenomenon in ferromagnetic materials, in which the magnetic energy of the crystal is decomposed into an irreversible and a reversible part, and uses a probability function to represent the interaction of domain wall motion and stray fields at the pinning center. The parameters of magnetic characteristics are related to spontaneous magnetization, magnetocrystalline anisotropy, magnetostriction, and microstructure, and can account for the effects of stress, temperature, and magnetization direction on hysteresis. In addition, compared with the J-A model, the energetic model also has the advantages of faster solution rate and fewer parameters. In addition to the above-mentioned commonly used scalar hysteresis models, there are also some vector hysteresis models, such as the Enokizono and Soda (E&S) model [15]. The vector hysteresis models can be constructed by superimposing the scalar models in different directions, such as the vector Preisach model [16], the vector J-A model [17].
Silicon steels have been widely made as the core of transformers and motors. However, with the development of power systems, the working environment of power equipmentalso becomes more complex, especially at high-frequency, high temperature, stress, non-sinusoidal excitation. The commonly-used hysteresis models and loss calculation methods are difficult to simulate magnetic properties under complex conditions. Therefore, it is necessary to establish accurate and efficient models suitable for various complex excitation conditions to quickly simulate the magnetic properties of materials.

Dynamic Preisach Hysteresis Model
The basic idea of the classical Preisach model is that ferromagnetic materials are composed of a large number of magnetic dipoles with hysteresis properties, and in the case of an output value u(t) at a given moment, the hysteresis characteristics of ferromagnetic materials produce the sum of the hysteresis characteristics of all magnetic dipoles. The output value of the magnetic dipole γαβu(t) is assumed to be either +1 or −1 values, µ(α,β) is a distribution function of α, β, also known as the Preisach function ( Figure 1). The classical Prcisach model can be represented by the following equation:  Based on the domain theory and statistical principles, the loss separation theory proposed by Bertotti assumes that the total loss of soft magnetic materials is divided into hysteresis loss, classical eddy current loss and excess loss, as follows: where: W is total loss, W hy is hysteresis loss, W cl is eddy current loss, W ex is excess loss. According to the theory of ferromagnetism, hysteresis loss is caused by the irreversibility of domain wall movement and domain rotation. The hysteresis loss of the material is not affected by the frequency of magnetic field. However, eddy current loss and excess loss are greatly affected by frequency. Thus, hysteresis loss is usually calculated at an extremely low frequency.
Eddy current loss is caused by eddy current during the alternating magnetization process. The skin effect of the magnetic field can be ignored at low frequencies, and the magnetic field inside the magnetic material is uniformly distributed. Eddy current loss can be solved according to the specific shape of magnetic materials. The single-cycle eddy current loss per unit volume W cl and the corresponding loss power per unit volume P cl can be expressed as: where: σ is the electrical conductivity of the material, in S/m, d is the thickness of a single, in m, f is the magnetization frequency, in Hz. Excess loss is caused by the dynamic rotation and translation of the domain inside the material under the applied magnetic field, which is caused by the microscopic eddy current. The expressions of the single-cycle excess loss per unit volume W ex and the loss power P ex per unit volume are obtained as: where: S is the cross-sectional area of the magnetic material, in m 2 ; G is a dimensionless coefficient (G = 0.1375); V 0 is a statistical parameter characterizing the local magnetic field distribution. The total loss W can be expressed as: (7) where: H hy , H cl and H ex are magnetic field strength components corresponding to hysteresis loss, eddy current loss and excess loss, respectively; l is the B-H closed path in a single period.
A dynamic Preisach hysteresis model with magnetic flux density B as the independent variable can be obtained: When the magnetic field excitation is sinusoidal, the peak value of the magnetic flux density is B p , eddy current loss W cl (B p ) and excess loss W ex (B p ) can be obtained by integrating (3) and (5), respectively: Hysteresis loss W hy is not affected by frequency f. The total loss minus eddy current loss (W − W cl ) is represented by: It can be seen that (W − W cl ) is proportional to the square root of the frequency. The intersection of the corresponding function and the ordinate axis is the hysteresis loss, and the slope k is 8.76(σSGV 0 ) 0.5 B 1.5 p . Therefore, the slope k can be obtained from the total loss measured by the experiment and the calculated eddy current loss, and then the statistical parameter V 0 can be obtained.
In this paper, the magnetic loss of silicon steel sheets under sinusoidal excitation at different frequencies from 5 to 500 Hz is measured. Then, the relationship between (W − W cl ) and the square root function of frequency is obtained by calculating eddy current loss as shown in Figure 2. The solution of the statistical parameter V0 under any peak value of magnetic flu density is interpolated at discrete points, as shown in Figure 3. The solution of the statistical parameter V 0 under any peak value of magnetic flux density is interpolated at discrete points, as shown in Figure 3. The solution of the statistical parameter V0 under any peak value of ma density is interpolated at discrete points, as shown in Figure 3.

Methodology
The accuracy of electromagnetic calculations depends not only on the a electromagnetic simulation, but also on the accurate characterization of th properties of materials. The magnetic properties of materials measured unde conditions are not suitable for non-standard conditions. Therefore, it is necessa ure the magnetic properties under non-standard conditions. In this paper, a s tester (SST, 500 × 500 mm) is used to measure the magnetic properties of silicon B27R090 under different excitation conditions. Figure 4 shows a magnetic me system, including SST, power supply, feedback control system and PC.

Methodology
The accuracy of electromagnetic calculations depends not only on the algorithm of electromagnetic simulation, but also on the accurate characterization of the magnetic properties of materials. The magnetic properties of materials measured under standard conditions are not suitable for non-standard conditions. Therefore, it is necessary to measure the magnetic properties under non-standard conditions. In this paper, a single sheet tester (SST, 500 × 500 mm) is used to measure the magnetic properties of silicon steel sheets B27R090 under different excitation conditions. Figure 4 shows a magnetic measurement system, including SST, power supply, feedback control system and PC.  A specific waveform signal is input to the power supply using the cor software of the measurement system, and the corresponding excitation of the p coil is applied after amplification by the power amplifier. According to the cuital law, the expression for the magnetic field strength H is A specific waveform signal is input to the power supply using the corresponding software of the measurement system, and the corresponding excitation of the primary side coil is applied after amplification by the power amplifier. According to the Ampere circuital law, the expression for the magnetic field strength H is (16) where: N 1 is the number of turns of the primary side coil, l m is the equivalent magnetic circuit length, i dc is the DC excitation current; i ac (t) is the AC excitation current and it is necessary to apply an AC voltage u(t) containing harmonics to the primary side coil and the expression u(t) is as follows: where: U 1m is the amplitude of voltage at the fundamental frequency; U km is the kth harmonic amplitude of voltage; f is the fundamental frequency set as 50 Hz in the paper; θ k is the phase difference between the kth harmonic and the fundamental wave. In this paper, the harmonic content η is defined as: where: B km is the amplitude of the magnetic flux density of the kth harmonic, B 1m is the amplitude of voltage of the fundamental frequency. The magnetic flux density can be obtained by integrating the voltage of secondary side, and the expression of the magnetic flux density is as follows: The main parameters of the magnetic measurement system used in the measurement are shown in Table 1. In this paper, the loss characteristics of silicon steel sheets under different bias currents are measured. The harmonic content of each harmonic is 30% of U 1 , and all the phase angles are set to be zero. Figure 5 shows the loss characteristics under different bias currents and third harmonics. As the magnetic field strength H dc increases from 0 A/m to 100 A/m, the loss also increases. However, as the magnetic field strength increases, the rate of the loss change gradually decreases.
In this paper, the loss characteristics of silicon steel sheets under different bias currents are measured. The harmonic content of each harmonic is 30% of U1, and all the phase angles are set to be zero. Figure 5 shows the loss characteristics under different bias currents and third harmonics. As the magnetic field strength Hdc increases from 0 A/m to 100 A/m, the loss also increases. However, as the magnetic field strength increases, the rate of the loss change gradually decreases.

Fitting Results under Single Harmonic Excitation
The dynamic Preisach model established above is used to simulate the dynamic hysteresis loops of silicon steel sheets B27R090. Figure 6 shows the comparison of dynamic hysteresis loop simulation and measurement results with varying the harmonic contents and the amplitude of the fundamental magnetic flux density. As can be seen from Figure  6a, when the third harmonic content η is 20%, the hysteresis loop is only distorted inwardly without local hysteresis loops. However, when the third harmonic content is 40%, the hysteresis loops appear as two identical local hysteresis loops at symmetrical positions, as shown in Figure 6c. From the fitting accuracy of the hysteresis model to the dynamic hysteresis loop, when Bp is 1.7 T, the simulated hysteresis loop is approximately the same as the measured result. Comparing between Figure 6b,d, it can be seen that the hysteresis model can still maintain the accuracy of the fitting results even with a local hysteresis loop in the hysteresis loop as η changes.

Fitting Results under Single Harmonic Excitation
The dynamic Preisach model established above is used to simulate the dynamic hysteresis loops of silicon steel sheets B27R090. Figure 6 shows the comparison of dynamic hysteresis loop simulation and measurement results with varying the harmonic contents and the amplitude of the fundamental magnetic flux density. As can be seen from Figure 6a, when the third harmonic content η is 20%, the hysteresis loop is only distorted inwardly without local hysteresis loops. However, when the third harmonic content is 40%, the hysteresis loops appear as two identical local hysteresis loops at symmetrical positions, as shown in Figure 6c. From the fitting accuracy of the hysteresis model to the dynamic hysteresis loop, when B p is 1.7 T, the simulated hysteresis loop is approximately the same as the measured result. Comparing between Figure 6b,d, it can be seen that the hysteresis model can still maintain the accuracy of the fitting results even with a local hysteresis loop in the hysteresis loop as η changes.  Figure 7 shows the comparison of dynamic hysteresis loop simulation and measurement results under various excitation conditions when k is 5 and the harmonic phase difference θ5 is 0°. It can be seen from Figure 7 that when the fifth harmonic content is 20%, there are four local hysteresis loops that are independent of each other, and when the fifth harmonic content is 40%, the four local hysteresis loops become two-crossed, as shown in Figure 7c. Compared with Figure 6, it can be seen that the hysteresis loops under the ex-   Figure 7 shows the comparison of dynamic hysteresis loop simulation and measurement results under various excitation conditions when k is 5 and the harmonic phase difference θ 5 is 0 • . It can be seen from Figure 7 that when the fifth harmonic content is 20%, there are four local hysteresis loops that are independent of each other, and when the fifth harmonic content is 40%, the four local hysteresis loops become two-crossed, as shown in Figure 7c. Compared with Figure 6, it can be seen that the hysteresis loops under the excitation of the fifth harmonic with the same content of 20% appear a local hysteresis loop. It can be seen that with the increase in harmonic order, the influence of harmonics on the shape of hysteresis loop becomes more complex. From the perspective of the fitting accuracy of the hysteresis model on the dynamic hysteresis loop, though the local hysteresis loops become crossed, the hysteresis model can still simulate the hysteresis loop accurately.
(c) k = 3, η = 40%, Bp = 0.8 T ( d) k = 3, η = 40%, Bp = 1.7 T Figure 6. Comparison between measured and calculated hysteresis loop under k = 3 (θ3 = 0°). Figure 7 shows the comparison of dynamic hysteresis loop simulation and measurement results under various excitation conditions when k is 5 and the harmonic phase difference θ5 is 0°. It can be seen from Figure 7 that when the fifth harmonic content is 20%, there are four local hysteresis loops that are independent of each other, and when the fifth harmonic content is 40%, the four local hysteresis loops become two-crossed, as shown in Figure 7c. Compared with Figure 6, it can be seen that the hysteresis loops under the excitation of the fifth harmonic with the same content of 20% appear a local hysteresis loop. It can be seen that with the increase in harmonic order, the influence of harmonics on the shape of hysteresis loop becomes more complex. From the perspective of the fitting accuracy of the hysteresis model on the dynamic hysteresis loop, though the local hysteresis loops become crossed, the hysteresis model can still simulate the hysteresis loop accurately. The hysteresis loops and losses are simulated by the dynamic hysteresis model under various excitation conditions. Table 2 shows the comparison between the model calculation loss result Pcal and the measured loss result Pmea under each excitation condition when the harmonic order k = 3 and 5. It can be seen from Table 2 that the loss increase as increasing the harmonic order and the harmonic content. Under these different excitation conditions, the relative error of the model calculation loss is within 4%, which meets the accuracy requirements.  The hysteresis loops and losses are simulated by the dynamic hysteresis model under various excitation conditions. Table 2 shows the comparison between the model calculation loss result P cal and the measured loss result P mea under each excitation condition when the harmonic order k = 3 and 5. It can be seen from Table 2 that the loss increase as increasing the harmonic order and the harmonic content. Under these different excitation conditions, the relative error of the model calculation loss is within 4%, which meets the accuracy requirements.

Simulation of Hysteresis Characteristics under AC-DC Hybrid Excitation
The fitting accuracy of the Preisach hysteresis model mainly depends on the Everett function. When the excitation contains DC component, the limiting hysteresis loop of the soft magnetic material will become asymmetrical, the classical method is no longer applicable and need to be improved. For any point P on the GPQ, the distance to the ascending branch is ΔH, then the magnetic field strength HP at point P is:

Hysteresis Model under DC Bias
During the movement of point P from point G to point Q, ΔH decreases continuously from ΔHg to zero, which can be expressed as a function of ΔHg and the limiting hysteresis loop width ΔHP = Ha(Bp) − Hd(Bp), and finally obtain HP as: According to the above improved method, the FORCs of B27R090 under the excitation of AC magnetic flux density peak value of Bacp = 1.9 T, AC frequency f = 5 Hz, and DC bias magnetization of Hdc = 50 A/m is obtained, as shown in Figure 9. The Everett functions based on the single curve of the descending branch and the ascending branch are obtained by solving, respectively, as shown in Figure 10. For any turning point G on the descending branch, the ascending branch and the distance ∆B g to B = B Q are given, respectively: For any point P on the GPQ, the distance to the ascending branch is ∆H, then the magnetic field strength H P at point P is: During the movement of point P from point G to point Q, ∆H decreases continuously from ∆H g to zero, which can be expressed as a function of ∆H g and the limiting hysteresis loop width ∆H P = H a (B p ) − H d (B p ), and finally obtain H P as: According to the above improved method, the FORCs of B27R090 under the excitation of AC magnetic flux density peak value of B acp = 1.9 T, AC frequency f = 5 Hz, and DC bias magnetization of H dc = 50 A/m is obtained, as shown in Figure 9. The Everett functions based on the single curve of the descending branch and the ascending branch are obtained by solving, respectively, as shown in Figure 10.
For any point P on the GPQ, the distance to the ascending branch is ΔH, t magnetic field strength HP at point P is:

H = H B -Ψ ΔH ,ΔH
According to the above improved method, the FORCs of B27R090 under th tion of AC magnetic flux density peak value of Bacp = 1.9 T, AC frequency f = 5 Hz, bias magnetization of Hdc = 50 A/m is obtained, as shown in Figure 9. The Everett fu based on the single curve of the descending branch and the ascending branch are o by solving, respectively, as shown in Figure 10.  The static hysteresis loop at Hdc = 50 A/m is fitted by the improved method, and the peak values of the AC magnetic flux density are 1.0 T and 1.5 T, respectively, and the results are shown in Figure 11. It can be seen that the model fitting results are approximately the same as the measurement results. The static hysteresis loop at H dc = 50 A/m is fitted by the improved method, and the peak values of the AC magnetic flux density are 1.0 T and 1.5 T, respectively, and the results are shown in Figure 11. It can be seen that the model fitting results are approximately the same as the measurement results.
The static hysteresis loop at Hdc = 50 A/m is fitted by the improved method, a peak values of the AC magnetic flux density are 1.0 T and 1.5 T, respectively, a results are shown in Figure 11. It can be seen that the model fitting results are ap mately the same as the measurement results. When the DC bias Hdc ranges from 10 A/m to 70 A/m, the difference between th loss and eddy current loss (W − Wcl) is calculated, and V0 is obtained by solving th of the linear regression function, as shown in Figure 12.
When the peak value of the AC magnetic flux density Bacp is less than 1.0 T creases with the increase in Hdc at the same Bacp, and when Bacp is greater than 1.0 T at different bias currents almost keep constant. Therefore, V0 is determined by bo and Hdc under AC-DC mixed excitation. When the DC bias H dc ranges from 10 A/m to 70 A/m, the difference between the total loss and eddy current loss (W − W cl ) is calculated, and V 0 is obtained by solving the slope of the linear regression function, as shown in Figure 12.  When the peak value of the AC magnetic flux density B acp is less than 1.0 T, V 0 increases with the increase in H dc at the same B acp , and when B acp is greater than 1.0 T, the V 0 at different bias currents almost keep constant. Therefore, V 0 is determined by both B acp and H dc under AC-DC mixed excitation.
Compared with the improved method, loss separation theory and V 0 fitting results under DC bias, the dynamic hysteresis loop under AC/DC hybrid excitation is simulated, and the results are shown in Figure 14, the peak value of the AC magnetic flux density B acp = 1.7 T in each excitation condition. Compared with no bias current, the dynamic hysteresis loop under AC/DC hybrid excitation is mainly distorted in the upper and lower parts, and the middle part is still greatly affected by harmonics. The comparison between the loss calculation result P cal and the measurement result P mea is listed in Table 3. It can be seen that the improved model is more suitable for simulating the dynamic hysteresis loop under the AC-DC mixed excitation. Compared with the hysteresis loop, it can be seen that the DC bias leads to the magnetization saturation only in the bias direction, which increases the magnetic loss of soft magnetic materials.  Compared with the improved method, loss separation theory and V0 fitting under DC bias, the dynamic hysteresis loop under AC/DC hybrid excitation is sim and the results are shown in Figure 14, the peak value of the AC magnetic flux den = 1.7 T in each excitation condition. Compared with no bias current, the dynamic h sis loop under AC/DC hybrid excitation is mainly distorted in the upper and lowe and the middle part is still greatly affected by harmonics. The comparison betw loss calculation result Pcal and the measurement result Pmea is listed in Table 3. It seen that the improved model is more suitable for simulating the dynamic hystere under the AC-DC mixed excitation. Compared with the hysteresis loop, it can be s the DC bias leads to the magnetization saturation only in the bias direction, w creases the magnetic loss of soft magnetic materials.

1.
Based on the asymmetric limiting hysteresis loop, the first-order reversal curve under bias magnetic field is generated, and the parameter identification of the Preisach hysteresis model is realized, which can simulate the harmonic bias hysteresis characteristics more accurately.

2.
Considering the influence of DC bias magnetic field and AC flux density peak on excess loss, the function formula of relevant parameters in excess loss is constructed, and the parameters are extracted based on the experimental results, so as to realize the accurate simulation of hysteresis loop under DC bias conditions. 3.
When the DC magnetic field strength increases from 0 A/m to 100 A/m, the loss also increases. However, the speed of loss change gradually decreases as the magnetic field strength increases.

4.
The improved numerical simulation method in this paper is suitable for the simulation of asymmetric hysteresis characteristics. However, due to the different types of silicon steel sheets, the influence of DC bias on total loss and different types of loss will also be different. This paper considers the hysteresis characteristics of DC bias, but not other working conditions, such as the influence of temperature and stress on the hysteresis loop that need to be studied and improved in the future.