A comparative study of J–A and Preisach improved models for core loss calculation under DC bias

Power transformer is widely used in the power system with the rapid development of the power converters connected to the grid. When a transformer operates under DC bias conditions, its iron core loss increases significantly, causing local overheating and threatening the proper operation of the transformer. However, there are persistent difficulties in accurately assessing the core loss when the induction waveform is influenced by a DC bias. This paper first proposes improvements to the J–A and Preisach models to evaluate the core loss of the iron core under DC bias. Additionally, we incorporate the hysteresis models into the finite element method (FEM) by modifying the fundamental constitutive equations in the FEM model in order to perform a precise core loss/distribution calculation. To verify the accuracy of prediction, a transformer prototype with a laminated core is developed. The improved J–A-FEM and Preisach-FEM models were directly compared in terms of calculation accuracy, numerical implementation, and computational burden.

Power converters become more prevalent in the power system due to the large demand for renewable energy sources connected to the power grid [1][2][3] .The Power transformer is a crucial component in a power converter as it can provide galvanic isolation and voltage transformation.DC bias is a common operating state of power transformers where direct current penetrates the windings, causing the iron core to approach half-cycle saturation, distorting the magnetization waveform and resulting in increased core loss and subsequent overheating 4 .Predicting the magnetic hysteresis behavior and corresponding core loss under DC bias conditions is crucial for optimizing the design of the power transformer to achieve high efficiency and prevent localized overheating 5,6 .
Conventionally, the core loss is predicted by using empirical formulas based on the Steinmetz equation (SE) with parameters identified from experimental measurements under sinusoidal excitations.Although the modified empirical equations proposed by the existing works [7][8][9][10][11] provide better accuracy in core loss estimation, they are based on experimental measurements.While DC biased magnetization significantly affects and exacerbates the core loss, it is hard to obtain all of hysteresis loops and core loss under different DC biased conditions to fit the Steinmetz equation.The Steinmetz method is thus uneconomic and time-consuming and cannot be a general method for accurate calculation of core loss.
The hysteresis models can provide a more economic and efficient way of describing the core loss 12 .Among numerous hysteresis models, the Preisach and Jiles-Atherton (J-A) models are widely accepted 13 .The J-A theory is based on the inhibition of domain wall motion by pinning sites (i.e., inclusions, voids, crystal boundaries, and lattice defects), and is very helpful in describing the behavior of domain wall motion at the primary hysteresis stage 14 .The magnetization of ferromagnetic materials arises from the movement of magnetic domains.However, the J-A model only accounts for reversible and irreversible magnetization during the domain wall motion process and neglects the mechanism of domain rotation.Consequently, it is not enough to accurately predict core loss when ferromagnets are subjected to alternating fields 15 or more distorted excitations such as DC bias and harmonics.
On the other hand, the Preisach model describes the hysteresis phenomenon through an infinite set of magnetic dipoles, which relate to the mechanisms of rotation of aligned moments within a domain 16 .The Preisach model has a remarkable prediction power as far as the description of hysteresis behavior is concerned, however since the magnetic dipoles (rectangular hysteron) can only describe irreversible magnetization, the deviation would appear in the hysteresis behavior prediction with a lower amplitude of alternating flux density 17 .However, its identification procedure is still complex, a number set of reversal curves or symmetric minor loops should be measured to achieve an acceptable result, which is impractical under magnetization with DC bias from the engineering perspective.
Even though the precise hysteresis model is developed, the accurate core loss estimation is still challenging because the magnetic field inside the core is not always homogenously distributed in many practical applications such as solenoid cores, and gapped transformer cores.Therefore, the finite element method (FEM) is necessary to calculate the magnetic field distribution with a complex core geometry.Normally, in the FEM model, the hysteresis characteristic of the core is not considered and the core loss is approximately estimated by the pure empirical Steinmetz equation, which limits the accurate prediction of the core loss for many applications 18 .Therefore, it is necessary to consider combining the hysteresis model with finite element analysis.However, incorporating the hysteresis model into finite element analysis remains a challenge due to several reasons.Firstly, complex hysteresis models are required to accurately describe hysteresis behavior, which often includes multiple parameters that can be challenging to identify.Secondly, modifications to Maxwell's equations are required, increasing the difficulty of solving the system 19 .Lastly, due to the nonlinear characteristics of hysteresis models, solvers typically require multiple iterations to obtain stable solutions 20 .
This paper proposes improvements to enhance the capability of the J-A model and Preisach model for DC biased core loss calculation, as well as simplify the expression of the models and the parameter identification process.Furthermore, a method is presented to calculate core losses using the finite element method (FEM) incorporating the improved J-A and Preisach models.The accuracy and performance of the proposed methods under magnetization with DC bias are studied and compared.The fidelity of the calculation method under varying DC biased excitations is evaluated and validated using experimental test results obtained from the Epstein Frame and transformer prototype of the laminated core.Finally, a clear selection scheme of models for different application scenarios is presented.

Proposed models and implementation The J-A dynamic hysteresis model
In a magnetic core, the total core loss can be separated into hysteresis, classical eddy current and excess losses.The classical eddy current loss is proportional to the square of the time derivative of flux density, and for a laminated core, it can be expressed as where e is the thickness of the material, σ the conductivity of the material, and β is the geometric coefficient.
Due to the existence of magnetic domains, the local eddy currents are generated near the domain walls when the domain configuration changes under a dynamic excitation, resulting in excess loss.
According to Bertotti's statistical theory of losses 21 , the excess loss can be calculated by where n 0 is the effective number of active magnetic objects (MOs) when the excitation is DC, S is the crosssectional area of the steel sheet, G is the coupling constant without unit, and V 0 is the statistical coupling field parameter, determining the ability of applied field to increase n 0 with increasing frequency.From the assumption that (2) can be reduced as The J-A dynamic hysteresis model is constructed by combining the traditional J-A hysteresis model with the models of instantaneous eddy current and excess losses.
By the J-A theory, the magnetization M can be decomposed into an irreversible magnetization component, M irr , due to the discontinuous domain wall motion caused by the pinning effect, and a reversible magnetization component, M rev , due to the elastic domain wall bending, and the change in energy of ferromagnet can be expressed as where H e is the effective field, H e = H + αM, k the pinning parameter, δ the direction coefficient, δ = 1 for dH/ dt > 0, and δ = − 1 for dH/dt < 0. The term on left side represents the actual magnetostatic energy, and the right ( 1) (5) µ 0 MdH e =µ 0 M an dH e − µ 0 kδ dM irr dH e dH e side represents the energy input and energy lost to the pinning effect.The anhysteretic magnetization M an equation can be expressed as where M s is the saturation magnetization, a loop shape parameter, and α is the local field parameter.
To consider the eddy current and excess losses, the energy conservation formula of the traditional J-A model can be modified as where the terms on the right-hand side are the input energy and the hysteresis, eddy current and excess losses, respectively.
Because H e = H + αM and B = μ 0 (H + M), the inverse expression of the dynamic J-A model can be obtained as where c is the reversible magnetization parameter, k e = e 2 σ/2β, and k ex = (GSV 0 σ) 1/2 .Both k e and k ex can be deduced from the models of instantaneous eddy current and excess losses.
The influence of DC bias on the excess loss is significant and should be properly considered 22,23 .In 7,10 , an expression of excess loss parameter k ex considering the influence of DC bias fields is provided as follows: where k 1 , k 2 and k 3 are constants fitting to the measurement.
The solution method of the J-A dynamic hysteresis model is summarized in the flowchart of Fig. 1.Firstly, initialize the model the model parameters including Calculate dM an /dH e, dM irr /dH e , dM/dH e and dM/dH e in order to obtain dM/dB based on (8).Update M i by using M i = M i−1 + dB[dM/ dB] i , and H i according to the constitutive relation between B and H.If the step number doesn't reach the maximum, then apply the updated H and M to the initialization process.

Generalized Preisach model
The Preisach model describes the hysteresis behavior via an infinite set of magnetic dipoles that have rectangular elementary hysteresis loops with different switching values of magnetic field intensity 24 , which is related to the mechanism of domain rotation, as shown in Fig. 2a.
Based on this, a hysteretic relation M(H) by ( 6)  where S is the triangular region shown in Fig. 2b, μ(α, β) is the distribution function of the dipoles which can be obtained by interpolation of the first-order reversal curve.
In the view of the complexity of interpolation and much additional experimental work, a simplified method is presented to find the integral of the distribution function based on the measured limiting hysteresis loop 16,25,26 .
From the Preisach diagram on the descending trajectory shown in Fig. 3a, the magnetization for a given magnetic field strength can be expressed as In this case ascending trajectory as shown in Fig. 3b, the magnetization can be expressed as where H n is the field intensity of the n-th reversal point, and M(H n ) is the corresponding magnetization, T (α, β), as shown in Fig. 4, is a function of the ascending and descending magnetization trajectories, and can be calculated from a data table of the limiting hysteresis loop by where M asd (α) and M dsc (β) are the ascending and descending magnetizations of the limiting loop, α and β are the ascending and descending field intensity.
In order to circumvent the difficulty of determining the distribution function μ(α, β), the dependence of the α and β is separated, and ( 13) can be modified as The function F thus can be calculated as follows (10)    The reversible component can be expressed as where M an is anhysteretic magnetization which is expressed as (6).Obviously, the expression of M an contains 3 parameters which would make the parameter identification process very complicated.Here, a simple M an identification method only based on the limiting loop is applied 27 , which is expressed as where H asd (M) and H dsc (M) are the ascending branch and the descending branch of the limiting loop as shown in Fig. 6, and M an (H) can be obtained by the inverse H an (M).
Finally, the total magnetization M t within the material can be obtained as ( 15) When the field intensity H increases (ascending trajectory), the total magnetization M t is obtained by (18), while the field intensity H decreases (descending trajectory), M t is obtained by (19).However, in the numerical implementation, the flux density vector is directly obtained with the magnetic vector potential formulation, it is desired to express the model in its inverse form 28 .
The filed intensity in the inverse form can be obtained from flux density as where ΔB is the change of input flux density, i the time step, and dB/dH = μ 0 + μ 0 dM/dH.The first order derivatives of descending and ascending magnetizations with respect to H can be calculated from ( 11) and ( 12) by where the H n is the return value with the n-th reversal point of input flux density.
The derivatives of total magnetization M t with respect to H can be obtained as where c is the reversible magnetization parameter which is the same as J-A theory, M an (H) the inverse expression of H an (M).
The solution method of the inverse generalized Preisach model is summarized in the flowchart of Fig. 7. Firstly, the model parameters, ascending branch, descending branch and M an are initialized.After calculating △B i , check the stack to obtain B n and H n .After finishing the calculation of F(H), F(− H), F(H n ), and F(− H n ), the stack status needs to be checked.If the stack is empty, calculate dM/dH on an ascending trajectory when B increases or calculate dM/dH on a descending trajectory when B decreases.If the stack is not empty, then initialize the magnetization stack.Once obtain the dM t /dH and dB/dH, calculate H i+1 and save the data of H, M and B. If the iteration step doesn't reach the maximum, then reiterate to check the stack with the updated H.

Implementation in FEM
By the Maxwell's equations, the magnetic field in terms of the vectorial magnetic potential A in domain Ω is governed by where M is the magnetization obtained by the hysteresis model, µ 0 the permeability of vacuum, and n is the unit normal vector, and J is the excitation current density, Γ 1 and Γ 2 are boundaries of the first and second boundary conditions, which suggest that on the boundary in domain Ω, the vectorial magnetic potential A is zero and the tangential component of the magnetic field intensity H is zero.
The solution space is discretized by using the Galerkin weighted residual method as: where N i (i = 1,…,m) is the nodal shape functions, and m is the total number of nodes.For 2D plane fields, A and J only contain the component in the z direction, ( 25) is expressed as The average flux density in the core is sinusoidal.Based on triangular mesh, the integral form of ( 26) is discretized as follows: where N ei donate the number of elements, S e the element area, S ij e = 1/(4S e )[c i c j + d i d j ], c i = y j − y k , d i = x k − x j , (i,j,k = 1,2,3), x, y the vertex coordinates of elements.
The matrix form of (28) in one triangular element is The total solving domain matrix equation can be expressed as the following where [M] is the magnetization matrix obtained by the hysteresis model at every node, [K] the stiffness matrix K ij e = v 0 S ij e , [F] the excitation matrix, and subscripts n and n-1 denote the n-th and (n-1)-th iterations.The total core loss can be obtained as   where i denotes the number of finite elements, S i e is the i-th element area, S is the area of FEM model, T the period of excitation, ρ the density of the core material, and N the total number of elements.The simulation results can be obtained by using the commercial FEM software COMSOL in conjunction with MATLAB to consider the hysteresis characteristics of ferromagnetic materials.

Experimental measurement
In this section, two experiment schemes are carried out, both specimen and transformer prototype tests are used for verification of the proposed models, and a part of specimen test results are also used for parameter identification.

Specimen test
The specimen test has been carried out by Epstein Frame as shown in Fig. 8a.The B-H loops and corresponding core losses of the HGO silicon steel sheets, B30P105, were measured.Figure 8b shows some examples of B-H loops measured under different sinusoidal excitations of 1.1-1.8T at 50 Hz.The measured BH and loss density data of B30P105 are presented in the Supplementary file.The measured data of sinusoidal excitation are used to identify the model parameters to predict magnetic characteristics under DC bias.

Transformer prototype of the laminated core
To verify the accuracy of the two models in engineering applications, a single-phase transformer prototype of a laminated core as shown in Fig. 9 is used in the experiments.the design parameters of the laminated core and the detailed structural dimensions are listed in Table 1.
Figure 10 depicts the experiment setup for core loss measurement under normal and DC biased excitation.The DC power supply provides the DC current to simulate the DC biased condition; the AC power supply provides AC voltage to produce a sinusoidal induction, and the excitation frequency is 50 Hz; the power analyzer can observe core losses in real time.The flux density waveforms at specific points are measured by magnetic probes.Then the excitation current on the primary side and the induced voltage of the secondary side are measured.The magnetic field intensity, H, is calculated by the excitation current while the magnetic flux density, B, is calculated by the induced voltage.

Model parameter identification
The J-A dynamic model The J-A dynamic model totally has 7 parameters (Ms, α, a, c, k, k e , k ex .),all these parameters depend on the type of materials, which can be determined from measured data by Particle swarm optimization (PSO) fitting 29 .The detailed process of parameter identification of the J-A dynamic model has been provided in our previous work 30 , Some of the parameters of the J-A dynamic model are shown in Table 2.

The generalized Preisach model
The function F(H) and M an of the generalized Preisach model are expressed only in terms of limiting hysteresis loop, and thus only the limiting loop is required for the model implementation as shown in Fig. 11.
It can be argued that the identification procedure of the J-A dynamic model parameters is more challenging compared to the generalized Preisach model.

Model accuracy comparison based on specimen test
To verify the correctness and accuracy of the generalized Preisach model and dynamic J-A model, the hysteresis loop is simulated and compared with the measured data of Epstein Frame in this part.
Comparing the calculated and measured loops at B peak = 1.6 T, both improved models are more accurate than traditional ones as shown in Fig. 12.When comparing the two models, only a slight difference between the two is apparent as shown in Fig. 13, and the core loss error of the generalized Preisach model is − 1.35%, and that of the J-A dynamic model is 1.35%.Both are acceptable and valid in engineering applications.
To evaluate the performance of the model under DC biased excitations, the predicted DC biased hysteresis loops are compared with the measured data in Figs. 14 and 15.Observed deviations of the J-A model appear at the turning part of the DC biased loop, as the red box shown in Fig. 14a.The enlarged view of the red box is depicted in Fig. 14b to make the deviation clearer.
With the increase of the DC biased fields, a bigger deviation of the J-A model appears as shown in Fig. 17, while simulated results of the Preisach model basically agree with the experiment results.Table 3 tabulates the error between measured and predicted results by the two models.
The comparison suggests that both models have good accuracy under sinusoidal excitation, but errors of J-A increase when the model is subjected to DC biased excitations.Even if the parameters of the J-A model are modified as (9), obvious deviations have still occurred as shown in Figs.14b and 15b, which means the prediction ability of the J-A dynamic model is worse than that of the generalized Preisach model as far as the DC bias is concerned.This is because the two models are based on different physical mechanisms: The J-A theory is based on the inhibition of domain wall motion by pinning sites, while the Preisach theory is related to the mechanism of domain rotation.The two physical processes dominate different magnetization stages, from the origin up to the  "knee" point of the normal magnetization curve, in which the domain wall motion is the main magnetization process, and from the "knee" point up to the full saturation, in which the magnetization process is dominated by the domain rotation.Thus, the J-A model has good accuracy in simulating the magnetization process with low peak flux densities (0 T-1.6 T), but the Preisach model is more accurate with higher flux densities (1.7 T-1.9 T).Moreover, the anhysteretic magnetization of J-A (M an ) which determines the hysteresis loop shape is only a single-valued function without a physical background, which further limits its application.The comparison of calculated and measured results of the specimen test demonstrates the Preisach model outperformed the J-A model in simulating the half-cycle saturation process due to its intrinsic advantage.

Local flux distribution comparison based on transformer prototype
Due to the hysteresis characteristics, when the applied field crosses zero, the flux density will not immediately drop to zero, and thus the residual flux density is produced.The residual flux density is an important criterion to verify whether the hysteresis models are successfully combined with the FEM. Figure 16 depicts the variations of flux density distribution at different instants and the residual flux density does exist.The local residual flux density distribution in the corner regions is non-uniformed, and it is closely related to the accuracy of the models in describing the hysteresis characteristics.Figures 17, 18  Further insight regarding the difference is obtained by analyzing the hysteresis loops at the specific points in the corner region of the core.Figure 20 compares the measured and calculated hysteresis loops by traditional   The calculated hysteresis loops of point A are compared with the measured results, and the operating point at 30 ms (the same as 10 ms) is marked in Fig. 22.The results show that the deviation at the turning part of hysteresis loops leads to the difference between the two models in the flux density distribution of the corner region.The error of flux density calculated by the J-A dynamic model at 10 ms is − 15% and that of the generalized Preisach model is 18%.Although the simulated results do not match exactly with the measured results, the accuracy of the J-A model is a little higher than the Preisach model.
Figure 23 compares the calculated and measured DC biased hysteresis loops of point B (marked in Fig. 19) where significant differences between the two models are obtained.As shown, the deviation between the generalized Preisach model and measurement is much less than that of the J-A dynamic model when the excitation is subjected to DC bias.In this verification scheme, the calculated flux distribution of the J-A dynamic model is a little closer to the measurement results under normal operating conditions, but as far as the DC bias is concerned, the generalized Preisach model has a much better performance.

Core loss distribution comparison based on transformer prototype
Besides the local flux density prediction, the total core loss and its distribution should also be considered to evaluate the model performance.As demonstrated in 31 , in different magnetization directions, the magnetic hysteresis loops of materials will exhibit different shapes.However, in practical engineering, the core of a power transformer is made up of laminated sheets of Grain-oriented steel.In order to ensure optimal magnetic properties, the magnetic path always follows the grain-oriented direction of the GO steel sheets.Therefore, in the following FEM simulation, it is assumed that the main magnetization direction of the transformer cores is along with the grain-oriented direction of the GO steel sheets.The average calculation time of Preisach is about 100 s per cycle and J-A is the 60 s in this FEM simulation.
Figure 24 depicts the core loss distribution in the prototype at B AC = 1.6 T. The core loss distribution is reasonably uniform except at the four corners, and the inner side of the corner appears to have more core losses.The calculated results of the two models are very close in this case, the predicted core loss by the generalized Preisach model is 1.58 W/kg and the J-A dynamic model is 1.50 W/kg at the inner side of the corner.In the limb region, the predicted core loss by the generalized Preisach model is 1.11 W/kg and the J-A dynamic model is 1.12 W/kg.
Figure 25 illustrates the distribution of core loss in the core of the prototype under a DC bias of 102.8 A/m and an alternating induction of 1.6 T. It is observed that the core loss values have increased compared to those without DC bias, but the loss distribution remains similar.

Figure 1 .
Figure1.Flowchart for the solution method of the J-A dynamic model.

Figure 7 .
Figure 7. Flowchart for the solution method of the inverse generalized Preisach model.

Table 1 .Figure 10 .
Figure 10.Experiment setup for measurements of normal and DC biased core loss.

Figure 13 .Figure 14 .Figure 15 .
Figure 13.Comparison of measured and calculated hysteresis loop by generalized Preisach model and dynamic J-A model.
and 19 compare the nonuniform fields predicted by the two models at 10 ms (T = 20 ms) with an alternating flux of amplitude 1.6 T and DC biased fields 0 A/m, 25 A/m, and 100 A/m, respectively.The DC biased fields 25 A/m, and 100 A/m are produced by injecting the direct current 0.5 A and 2 A, as depicted in Fig.10.As shown, the flux distributions of the two models are quite different at 10 ms, and the difference becomes larger as the bias field increases.

Figure 16 .Figure 17 .
Figure 16.The flux density distribution in the corner region at different instants.
).As shown, Both the improved models are more accurate.Figure21depicts the measured H-, B-waveform (H reduced by 33 times) and the hysteresis loop at point A, as shown, the value of flux density at 10 ms (or 30 ms) corresponds to the point b of the hysteresis loop.

Table 2 .
Parameters of J-A dynamic hysteresis model.

Table 3 .
Comparison of calculated results with Epstein frame tests.

Table 4 .
Comparison of predicted and measured results.