Study on Dynamic Hysteretic and Loss Properties of Silicon Steel Sheet Under Hybrid Harmonic and DC Bias Excitation

The dynamic and static hysteresis properties as well as loss properties of grain oriented (GO) silicon steel sheet are measured under the AC-DC hybrid excitation containing both DC bias and harmonic components. The influence of the DC-biasing magnetic field, harmonic orders and phase shift angles on the iron loss is analyzed in depth based on the experimental data under different excitations. A dynamic hysteresis model and the parameter identification method is presented to simulate the asymmetrical hysteresis curve with minor loops. According to the equivalence between loss separation and field separation, the static part of magnetic field intensity are calculated by an improved Preisach model, while the dynamic parts are obtained by fitting the classical loss as well as excess loss. The consistency between simulation and measurement verifies the effectiveness and accuracy of the proposed method.


I. INTRODUCTION
Soft magnetic materials are widely used in power equipment and electronic devices because of their high permeability and low loss [1]. The iron core of large-scale motors, transformers, reactors and other equipment is usually made of silicon steel sheet and works under the power frequency sinusoidal excitation. However, with the rapid development of highvoltage direct current (HVDC) transmission technology and the comprehensive application of power electronic elements, the DC bias or rich harmonics generated in power system and the distorted magnetic flux in magnetic core leads to the increase of iron loss, aggravates local overheating, and finally endangers the safe operation of the power system [2], [3]. Hence it is of great significance to investigate the hysteretic behavior and loss property under complicated excitations for optimal design and safe operation of electrical equipment.
The Steinmetz's equation and loss separation theory are two main approaches for modeling iron loss of soft magnetic The associate editor coordinating the review of this manuscript and approving it for publication was Zhe Zhang . materials. The former one was originally presented to estimate loss under sinusoidal excitations [4]. Afterwards the Steinmetz's equation was improved for loss modeling under harmonic excitations [5]- [7]. Different improved models were compared in detail by analyzing the limitations and advantages in application [8]. The Steinmetz's premagnetization graph was proposed and developed to estimate the core loss under DC bias conditions at high frequencies and flux densities [9], [10]. Besides, the loss separation theory, where the total loss is expressed as the sum of static and dynamic components, is also an important way to estimate the iron loss of soft magnetic material. The static or dynamic components of iron loss are assumed to be polynomial functions with respect to flux density to investigate the loss properties of iron core under sinusoidal or DC-biased excitations [11]- [13]. Another approach to estimate iron loss under harmonic or DC-biased magnetization is the field separation [14]- [16], which is equivalent to the loss separation. The static loss components is calculated by hysteresis model, and the modeling of dynamic loss requires the shape controlling factor g(B) of dynamic loop [17] or the statistical VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ parameter V 0 associated with the distribution of local coercive field [18]- [20]. More and more electrical devices operate under the hybrid excitations including DC bias and harmonics, such as HVDC converter transformers and high-frequency transformers [21]- [23], which makes the iron loss estimation and prediction a great challenge. When the DC-biased and harmonic magnetizations are simultaneously applied, the asymmetric hysteresis curves due to DC bias and minor loops as a result of harmonic magnetization requires excellent performance of hysteresis models in curve fitting, especially when the iron core is highly saturated. In addition, more appropriate methods for modeling dynamic loss should be used with the influence of frequency, amplitude of flux density as well as DC magnetic field taken into account. However, up to now much work pays attention to magnetic and properties under either DC bias excitation or harmonic excitation, which make it difficult to design the electrical equipment operating under the aforementioned complicated excitations.
In this article, the hysteretic and loss properties of GO silicon steel under AC-DC hybrid excitation are measured to investigate the influence of DC bias and harmonics on the iron loss. In order to achieve a more accurate simulation of static hysteresis, a new parameter identification method for asymmetric hysteresis loop simulation is proposed and then an improved Preisach hysteresis model is established. Furthermore, with the influence of DC magnetic field strength and harmonic order on the excess loss parameter V 0 considered, a dynamic hysteresis model for estimation and prediction of magnetic and loss properties under AC-DC hybrid excitations is presented. The fitted results of static and dynamic hysteresis loops are compared with the experimental ones, which verifies the proposed dynamic hysteresis model and numerical identification procedure.

II. HYSTERESIS MODEL UNDER AC-DC HYBRID EXCITATIONS
According to the statistical theory of losses (STL) developed by Bertotti, the total iron loss per cycle P s,t consists of three components, namely, hysteresis loss P s,hy , classical loss P s,cl and excess loss P s,ex [24].
Because of the complexity of hysteretic characteristic of soft magnetic materials, the hysteresis loss under the sinusoidal excitation is generally expressed as where k h is the coefficient of hysteresis loss, α is assumed to be a polynomial function with respect to the peaked flux density B m [12], and f m is the frequency of excitation.. Although the loss estimation can be achieved by using the loss separation method in (1), there are some inevitable defects that may impede its further application in a boarder range. Firstly the parameter identification requires abundant experimental data, but the performance of the method in loss prediction is often unsatisfactory. Secondly the loss separation cannot incorporate hysteresis nonlinearity into FEM to solve the nonlinear magnetic field more accurately.
Assuming that the magnetic field is uniformly distributed along the inside of the soft magnetic material, the classical loss, calculated from Maxwell's equation for a perfectly homogeneous material with on-domain structure, can be represented as follows, where d is the thickness of the GO silicon steel sheet, and σ is the electrical conductivity. W s,cl is the classical loss over a time period T . It should be noted that the skin effect is negligible when the frequency of the harmonic excitation is not high enough [18], [25], [26]. According to the STL, the excess loss arises from the microscopic local eddy current induced by the motion of the active domain walls. It is difficult to scale the physical action of domain walls that is related to dynamic magnetization, therefore a mathematical expression of excess loss is presented as follows [19], where S is the cross-sectional area of the sheet, G is a dimensionless coefficient, and V 0 is a statistical parameter with the dimensions of magnetic field, related to the microstructure characteristics of the magnetic material. W s,ex is the excess loss over a time period T . The field separation can be derived from the loss separation in (1), thus the dynamic magnetic field intensity H s,t is divided into three parts [19], [20].
where B is the magnetic flux density under AC-DC hybrid excitation. H s,hy , H s,cl and H s,ex are the corresponding magnetic field intensity of hysteresis loss, classical loss and excess loss, respectively. The relationship between field and energy can be represented as follows, The small increment of classical loss and excess loss during an extremely short time period dt can be respectively expressed as dW s,cl in (10) and dW s,ex in (11).
Combining (9) with (11) and (10) with (12), it can be derived as follows where the switch parameter δ is controlled by the derivative of magnetic flux density with respect to time i.e. dB/dt, which indicates the variation tendency of the flux density B.

III. IMPROVED PREISACH MODEL UNDER AC-DC HYBRID EXCITATIONS
A. MEASUREMENT OF MAGNETIC CHARACTERESTIC As shown in Fig. 1, a magnetic performance testing platform provided by BROCKHOUS is composed of a standard Epstein frame, a visual operating system, and a multi-functional power supply, which is employed to carry out the AC-DC hybrid test on the GO silicon steel sheet (see Table 1) in order to obtain the experimental data of hysteresis loops and iron loss under different superimposed magnetizations.

FIGURE 1.
Magnetic performance testing platform.
The operating system is responsible for data processing through the control instructions input by a PC terminal, and then the corresponding excitation is applied. The adjustable field excitation of primary winding can be realized by the signal provided by the power amplifier. The relationship between the exciting current i(t) and magnetic field intensity H (t) agrees with Ampere circuital theorem.
where N is the number of coil turns in Epstein frame. i(t) is the exciting current of the primary winding. l m is the equivalent length of the magnetic circuit. The DC magnetic field H dc is established in the sheet by applying direct current i dc into the primary winding. As shown in (14), H dc is linearly correlated with i dc . The harmonic voltage u(t) given in (15) is also applied to the primary winding, and then the induced voltage u 2 (t) can be obtained from the secondary winding.
where f , the fundamental frequency, is set to be 50 Hz to obtain dynamic hysteresis curves and loss under hybrid excitation. U 1m is the amplitude of fundamental harmonic voltage, U km is the amplitude of k-th harmonic voltage. φ k is the phase difference between the k-th harmonic and fundamental The ratio η is fixed to be 30% in experiment.
Eventually the flux density B(t) is calculated by integrating the measured data of u 2 (t). Here the peak value of B(t) is defined as B acm .
where u 2 (t) is the induced voltage of the secondary winding.

B. STATIC HYSTERESIS MODEL AND ITS PARAMETER IDENTIFICATION
To investigate one and two dimensional (1&2-D) hysteretic and loss properties of ferromagnetic materials, researchers has proposed various hysteresis models such as Jile-Atherton (J-A), Energetic, E&S, and Preisach models [14]- [16], [25]- [29]. Concerning the physical mechanisms of magnetization, the Preisach model has been presented and developed from the perspective of a stringent mathematical definition to accurately model the hysteresis nonlinearity and predict the magnetic characteristics [30], [31]. It is known that the hysteretic and loss characteristic depends closely on the specific forms of magnetizations. The concentric hysteresis loops exhibits the symmetric feature and there is only one major loop under sinusoidal excitation. The existence of DC magnetic field in the hybrid magnetization will make the hysteresis loops extremely asymmetric, and the high-order harmonic will result in minor loops during dynamic magnetization process. Therefore, further work is required for the improvement of hysteresis model and its identification strategies for accurate estimation of asymmetric hysteresis curves and iron loss under AC-DC hybrid magnetizations.
The Everett function is critical to the static hysteresis simulation by using the Preisach model, and it can be established through the first order reversal curves (FORCs) of ferromagnetic materials. Considering of the difficulties in obtaining FORCs by experiment, an alternative method to create FORCs was proposed in [32]. This method is based on the property of return-point memory, which is a most important and distinctive feature of history-dependent hysteresis behaviour in magnetization of soft magnetic material. Considering the asymmetric characteristic of static hysteresis loops in this article, a more satisfactory result of Everett function should be obtained by the FORCs generated from asymmetric major loops measured in DC-biased excitation rather than symmetric ones in sinusoidal cases.
It is difficult to obtain ideal static hysteresis curves through measurement due to the limitation of experimental condition. Therefore, in experiment the static hysteresis loops are measured at low frequencies (f = 3Hz). The measured static DC-biasing major hysteresis loop can be used for numerical generation of FORCs and the consequent identification of hysteresis model. Fig. 2 depicts the asymmetric major hysteresis loop under DC-biased magnetization with B acm = 1.9T and H dc = 100A/m. Two sets of FORCs are numerically generated by the aforementioned algorithm. It should be noted that the major hysteresis loop will be different when DC offset H dc is changed. In this case, the Everett function needs to be reconstructed. Compared with the FORCs under sinusoidal magnetization, it is obviously that the FORCs with reversal points starting from the ascending branch shown in Fig. 2(a) and that with reversal points starting from the descending branch shown in Fig. 2(b) are not centrosymmetric any more.
As a result, the inverted Everett function E a (b u , b v ) and E d (b u , b v ) can be constructed from the asymmetric FORCs with reversal points starting from the ascending branch and descending branch of major loop, respectively.
where b u and b v (or b v and b u ) are the magnetic flux density values of the starting points like A and C (or points to be solved like B and D) on the FORCs in Fig. 3, and their corresponding field intensity are denoted as well. According to the inverted Everett function (see Fig. 4) in the B-based model, the quasi-static magnetic field intensity can be calculated when B is increasing with time (dB/dt >0) as (20) and when B is decreasing (dB/dt <0) as where B Mk and B mk (k = 1, 2, . . . , n) are the sequences of local maximum values and minimum values, respectively. Besides, B M0 and B m0 represent the initial state of γ -operators, and are the coordinates of the right angel of the limiting triangle [33].
In order to testify the superiority of proposed method, the static hysteresis loops in Fig. 5 provide a comparison between measurement and simulation. We identify the Preisach model by using two different strategies. Firstly only one Everett function is established for identification of the Preisach model to simulate the asymmetric hysteresis curve with minor loops (denoted by Simulated I in Fig. 5), and then two different Everett functions are established by the proposed method for new simulation (denoted by Simulated II in Fig. 5). Obviously, when only one Everett function generated by symmetrical major loop is employed to identify the Preisach model, the simulated result is not satisfying, especially in the unsaturated region of ascending branch and the saturated region of descending branch. However, based on the ascending branch and descending branch of asymmetric DC-biasing major loop, the utilization of two Everett functions improves the fitting accuracy greatly, which decreases the relative error of static hysteresis loss P s,hy by more than 18%. Moreover, the static hysteresis loops under sinusoidal excitations are given in Fig. 6, which present a comparison between the proposed approach and sinusoidal measurements to show the accuracy of the model without any DC offsets. The relative errors of static hysteresis loss P s,hy are 3.28% and 0.15% under the sinusoidal conditions of B acm = 0.6T and B acm = 1.4T respectively. This further proves the applicability and accuracy of proposed method in the sinusoidal conditions.

IV. PARAMETER EXTRACTION OF DYNAMIC LOSS MODEL AND HYSTERESIS SIMULATION
Under the AC-DC hybrid excitation, the classical loss, depends on the derivative of flux density with respect to time dB/dt instead of DC magnetic field H dc , thus the magnetic field intensity H s,cl can be directly calculated by (12). According to (13), the parameter extraction of the dynamic field H s,ex is crucial for the excess loss modeling. It is pointed out that the dynamic magnetic field H s,ex as well as the statistical parameter V 0 are related to the peak AC flux density B acm and DC magnetic field strength H dc in a DC-biased condition [19]. Therefore the influence of the harmonic order k on V 0 should be considered under the superimposed excitation of DC bias and high-order harmonic.

A. PARAMETER EXTRACTION OF EXCESS LOSS
Experiments under various AC-DC hybrid magnetizations specified in Table 2 are carried out to provide enough data to fit the function V 0 (k, B acm , H dc ). As shown in Table 2, the DC magnetization, applied in incremental proportions in experiment, are represented by H dc,i . The AC excitation consisting fundamental wave and k-th harmonic is also applied in different situations, which is expressed by the step-increased peak value of alternating flux density B acm,j . The values of parameter V 0 extracted from the measured data of excess loss are shown in Fig. 7. It can be seen from VOLUME 8, 2020 Fig. 7(a) that V 0 increases with the incremental DC magnetic field H dc , whereas decreases with the enlarged AC flux density B acm and tends to stay unchanged in saturation region of the magnetic material. In Fig. 7(b), it can be seen that V 0 is approximately linear with the harmonic order k under different hybrid excitations. According to the variation tendency of V 0 with the aforementioned three factors, the following formulas based on exponential function and linear function are presented.
where g i is the linear function with respect to harmonic order. p mi , p ni and q i (i = 1,2,3,4) are the fitting coefficients to be determined by experimental data. The comparison between the fitted excess loss and the measured ones under two different magnetizations is shown in Fig. 8 and shows good agreement, which indicates that the proposed expression of parameter function V 0 in (22) can meet the desired requirements on accurate estimation of magnetic field and the corresponding excess loss.

B. SIMULATION OF DYNAMIC HYSTERESIS AND LOSS PROPERTIES
Gathering all the contributions, eventually we can obtain the dynamic hysteresis model H s,t (B) as given in (7) based on the field separation technique. Fig. 9 depicts the fitted and experimental dynamic hysteresis curves under different AC-DC hybrid excitation. The results indicate that in addition to the asymmetry led by DC bias, harmonic order can affect the number and shape of minor loops, thus aggravating the complexity of hysteresis and loss properties. By using the proposed dynamic hysteresis model, the simulation of hysteretic curves of the GO silicon steel sheet can be achieved for estimation of the iron loss.
Furthermore, the iron loss can be obtained by the enclosed area of dynamic hysteresis loops. The total iron loss P s,t under two different excitations (I: H dc = 25A/m, k = 3 and II: H dc = 75A/m, k = 9) calculated by the method proposed in this article (P s,cal−1 ) and in [32] (P s,cal−2 ) are compared with the experimental data (P s,mea ) in Table 3 and Table 4, respectively. It can be concluded that the proposed method in this article has better performance in dynamic iron loss simulation under AC-DC hybrid excitations, and the relative errors of calculated results are within 10%. From the comparison and analysis above, a conclusion that the proposed dynamic Preisach hysteresis model capable of estimating the hysteresis loss, classical loss, and excess losses simultaneously can be used for accurate simulation and prediction of hysteretic behavior and loss characteristic under AC-DC hybrid magnetizations, is drawn.

V. ANALYSIS ON IMPACT OF DIFFERENT FACTORS UPON LOSS PROPERTIES
According to the measured results of dynamic experiment, total iron loss can be analyzed from the perspective of DC  [32] and proposed method in this article (Excitation I).

TABLE 4.
Comparison between measured iron losses and calculated ones obtained by method in [32] and proposed method in this article (Excitation II).
magnetic field strength H dc , harmonic order k and phase shift angle φ k . As depicted in Fig. 9, the presence of DC bias and high-order harmonics in the excitation will lead to asymmetric hysteresis loops and possible minor loops which inevitably has great influence on the iron loss. It is necessary to observe how the iron loss varies with DC and harmonic excitations.
In Fig. 10, the loss curves under the hybrid excitations that includes the fundamental and seventh harmonic superimposed with different DC bias are demonstrated. The iron loss increases significantly when the DC magnetic field strength H dc increases from 0A/m to 25A/m; however, the growth rate is getting slower with the increased DC excitation since the iron has been highly saturated. In Fig. 11, the loss curves under the hybrid excitations including the fundamental harmonic and specific DC VOLUME 8, 2020 magnetic field (H dc = 100A/m) superimposed with different high-order harmonics, are depicted. It can be concluded that due to the enlarged eddy current in the silicon steel sheet at higher frequency, the higher the harmonic order is, the faster the loss increases. The loss curves are depicted in Fig. 12 when H dc and k change simultaneously. The loss curves under hybrid excitation with k = 3, H dc = 25A/m and k = 5, H dc = 0A/m intersect at point A. The DC bias affects the static hysteresis loss and excess loss, and will be limited by the saturation of core magnetic flux density. The harmonic order of hybrid excitation affects classical loss and excess loss, which is not restricted by the saturation characteristic however. This explains why the effect of harmonics order on the total loss is dominant when the peaked AC flux density B acm is high. In order to analyze the significant impacts of phase shift between the harmonic field and the fundamental field on the resulting iron loss, the proposed model is employed under the hybrid excitations (k = 3, H dc = 25A/m, B acm = 1.6T) with various phase shift angles. The same good agreement between measurement and calculation is observed in Table 5 when the phase shift angle φ k varies from 0 • to 360 • . The calculated results of loss separation are shown in Fig. 13.  As shown in Fig. 13, the hysteresis loss P s,hy , classical loss P s,cl and excess loss P s,ex increase firstly and then decrease with the increase of phase shift angle φ k . Overall, these losses are all approximately symmetrical about φ k = 180 • , resulting in the same symmetric feature of total iron loss P s,t .

VI. CONCLUSION
Magnetic and loss properties under AC-DC hybrid magnetizations are measured by an enhanced experimental platform and analyzed extensively to observe the effects of DC magnetic field and harmonic flux density on static and dynamic hysteretic characteristics. The asymmetric FORCs can be generated numerically based on the asymmetric major hysteresis loop from measurement under the DC-biased magnetization and it can be further utilized in the construction of asymmetric Everett function as well as the consequent improved static Preisach model. Furthermore, by extracting the statistical parameter of magnetic field intensity H s,ex , a new dynamic hysteresis model is presented with the influences of DC bias and high-order harmonic on the excess loss taken into account. The proposed model achieves high precision in simulation of asymmetric hysteretic curves with minor loops, and thus proves to be an appropriate choice in iron core loss estimation and prediction when the high-order harmonic excitation and DC bias are applied simultaneously.
Different from the loss property of GO silicon steel sheet subjected to sinusoidal excitation, inclusion of either the DC bias or the high-order harmonic in excitation will give rise to the increase of total iron loss. The DC-biased magnetization leads to the growth of iron loss by uplifting the magnetic flux in magnetic core to saturate the iron core significantly, while the high-order harmonic excitation increases the iron loss due to the enlarged eddy current and minor loops at higher frequency. It is noticed that the high-order harmonic may play a dominant role in variation of iron loss. In addition, iron losses will increase at first and then decrease with the increment of phase shift angle φ k , and are approximately symmetrical about φ k = 180 • .
RUI WANG was born in Tangshan, Hebei, China, in 1994. She received the B.S. degree in electrical engineering and automation from North China Electric Power University, in 2018, where she is currently pursuing the master's degree with the Department of Electrical Engineering.
Her research interests include operational safety in power equipment and numerical analysis of electromagnetic field. He has been engaged in basic research on transformers, electrical design, and key common technology for a long time. His key research interests include key components and material properties of transformers, key technology research on vibration and noise of transformers, reactors, and electromagnetic thermal simulation technology.
LICHUN DONG was born in China, in 1983. He received the B.S. degree in mechanical engineering from North China Electric Power University, in 2006.
His research interests include operational safety in power equipment and numerical analysis of electromagnetic field and material properties of transformers.
LANRONG LIU was born in 1978. She received the B.S. degree in electrical engineering from the Shenyang University of Technology, in 2000, and the M.S. degree in electric engineering from North China Electric Power University.
She has been engaged in basic research on transformers and key common technology for a long time. Her key research interests include key components and material properties of transformers, key technology research on vibration and noise of transformers, reactors, and electromagnetic thermal simulation technology.