Fractional model of magnetic field penetration into a toroidal soft ferromagnetic sample

We propose an original approach to solve the coupled problem of alternative magnetic field penetration inside a toroidal soft ferromagnetic sample and frequency dependent magnetic hysteresis. Local repartition of ferromagnetic losses depends on the instantaneous material properties and on the frequency of the excitation field waveform. A correct solution to the model, with respect to this repartition, implies a higher resolution in two dimensions of the diffusion equation including local dynamic hysteresis consideration. The resulting model gives precious local information but requires complex parameter setting, high computational capacity and long simulation time. Due to the toroidal shape, a single dimension algorithm solving of the diffusion equation is clearly insufficient and it would lead to inaccurate simulation results. Consequently, a large number of discretization nodes and extended simulation time must be considered in two dimensional configurations. In our alternative solution, starting from a lumped model, we add a fractional time derivative of the dynamic hysteresis losses. It leads to an accurate formulation of the problem with a reduction in complexity and simulation times.


Introduction
Local ferromagnetic losses through the cross section of a toroidal magnetic core cannot be physically measured.
The understanding and modeling of such losses is of a major interest as soon as we need accurate simulations of complete electromagnetic systems. Commonly used differential breakers could be examples of this case. The simulation of the magnetic field evolution and propagation should be considered with proper geometrical constraints and with presence of environing electronics. Fairly high response levels (amplitude, frequency) can be reached in transient phases. A correct treatment should also take into account material properties. In particular, the magnetic material law is necessary to obtain correct simulations of the whole system [1][2][3][4].
A number of research teams have already implemented successfully a hysteresis dynamic model in the diffusion equation of the magnetic excitation field [5][6][7][8][9]. The resolution in one dimension problem of such coupled consideration leads to accurate simulations. Unfortunately, it implies work on samples of specific geometry. For instance, ferromagnetic sheets which exhibit much lower thickness comparing to both remaining dimensions. In this article, a similar approach will be used in a two dimensional consideration which is necessary as soon as we work on toroidal geometry. The instantaneous resolution of the diffusion equation and of the dynamic hysteresis model provides very interesting information such as the repartition of the magnetic hysteresis losses through the cross section of the tested sample. It gives an accurate losses cartography of the tested toroidal sample. This sensitive and accurate simulation provides good results on a large frequency bandwidth. Unfortunately, such coupled techniques have some disadvantages. It is not modular. It lacks flexibility requires an important space in the memory allocation and usually excessive simulation times. An alternative option should be using a lumped model, without space consideration and focusing on the dynamic hysteresis behavior. Note that some parameters should be fit to geometrical constraints.In terms of frequency dependence or time derivative, the model can be a first order model (product of a constant to the time derivation of the induction field) or/and something a little more complicated such as taking into account of the excess losses by another contribution varying proportionally to the square root of the frequency. However, as demonstrated frequently these simple hysteresis models give correct behavior on a relatively limited frequency bandwidth [9].
In this article, we study frequently used configurations of toroidal soft magnet with applied oscillating magnetic field H. A schematics of the experimental configuration is shown in Fig. 1 To reveal the dynamical hysteretic behaviour we propose to use a lumped model where the space consideration is replaced by a fractional derivation taking into account the dynamic losses. The fractional order, using the system memory to parametrize the delayed response of the deeper regions, gives an alternative freedom in the simulation. We can perfectly fit or simulate the hysteresis curve area versus frequency curve on an extremely large frequency bandwidth by adjusting it [9][10][11].
To describe this alternative approach and to draw the main conclusions, this article has the following structure. After the current introduction (Sect. 1), a theoretical presentation of the two simulations scheme used during this study is provided (Sect. 2). Our first approach is constituted by the coupled two dimension discretization technique in finite difference/dynamic hysteresis model. The second consideration is a lumped model with a fractional derivative hysteresis. In the second part of this article, a large number of comparison simulation/results are given. Advantages and disadvantages of both simulations are also discussed. In Sect. 3, experimental results are presented and compared to simulations. Finnaly, Sect. 4 finish the paper with conclusions.

2D finite differences scheme for the simulation of the magnetic field diffusion
The diffusion equation is solved through the square crosssection of a toroidal magnetic flux. We assume unidirectional and inhomogeneous surface excitation field [according to the diameter of the magnetic core (see Fig. 1)]. In addition, we assume constant and homogenous electrical conductivity. By considering the dimensions of the studied magnetic core (namely width and thickness) a 2 dimensions study is carried out. Symmetrical considerations allow reducing the studied area to half the width of the magnetic flux. This consideration constitutes a significant reduction of the calculation time. The magnetic field diffusion results from the Maxwells equations [12]: where H is an applied magnetic field, B is the induction, and σ is electrical conductivity, As the magnetic field, in our configuation, is always perpendicular to a cross section of the toroidal sample (see Fig. 2), div (H) = 0, using Cartesian coordinates Eq. 2.1 becomes By using a single dimension finite differences temporal discretization of the temporal derivation of the magnetic induction in Eq. 2.2, the following expression appears: By considering a linear relation between B and H, one can obtain the analytical solution of Eq. 2.3. Unfortunately, this simple consideration gives simulation results far from the experimental ones. As soon as the frequency of the excitation field exceeds a few hertz, frequency effect appears. If we want to obtain correct simulation results under such external conditions, it implies taking into account the frequency dependence and the hysteresis phenomenon due to wall domain motions through the tested sample. A simple consideration of the magnetic dynamic hysteresis gives: where H stat (B) is a fictitious contribution derived from a quasi-static consideration of hysteresis. Namely, a quasistatic hysteresis model has been used to provide this contri-bution. A lot of quasi-static hysteresis models exist, however, the specific one we used is described in the next subsection. The parameter β depends on material and temperature and in our case is fitted from the experiment. It is independent of the geometry and of the excitation waveform. Former experimental validations have shown that this dynamic material law constitutes a good description of dynamic (damping) effects due to domain Blochs wall motions [6][7][8][9][10][11]. This material law represents a statistical behavior of the wall motions. It can then be considered as isotropic and characteristic of the material. Equation 2.4 compiles intwo hysteresis contributions, a quasi-static contribution described hereafter and a dynamic contribution product of the constant β and the time derivation of the induction field B.

Quasi-static hysteresis model
The quasi-static hysteresis model we used is based on the following assumption: under quasi-static excitation, we assume that Blochs domain wall movements behave like mechanical dry frictions. A static (frequency-independent) equation based on this assumption is proposed for the numerical simulation of these movements. A major hysteresis loop is obtained by translating a nonhysteretic curve. The sign of this translation is equal to the sign of the time derivative of the induction B and its amplitude which is equal to the coercive field, H c .
where, f (H ) (and reciprocally f −1 (B)) is a nonhysteretic and saturated function determined using an experimental major hysteresis loop (H max > H c ). B = B z and H = H z are defined on at a single point along z direction. The function f is characterized by the material parameters σ 1 and γ by comparison (simulation/measure) to a quasi-static major hysteresis loop: Note that, Eqs. 2.5 and 2.6 give a description of an analogue of a single parameter dry friction phenomenon. It is obviously not enough to describe correctly the whole range of tested sample behaviors related to a set of much larger contributions. However, this model can be considerably improved by taking into account a set of similar dry frictions with coeficients distributed in some specific domains. In our case each componnent is characterized by its own coercive field H ci and its own weight, depending on the position, in the final reconstitution of the magnetization and induction fields. More realistic cycles, including minor loops, are obtained by introducing a distribution of a basic element (spectrum), characterized by the corresponding coercive field and weight.
where k is number of componnents and Spectrum(i) represents the distribution of various elementary dry frictions.  (B(x, y, t))) (2.8)

Two dimension finite differences resolution
Thanks to the simplicity of a two dimensional case and the symmetries of the proposed problem (half width of the magnetic core as illustrated on Fig. 3), a formulation by finite differences method of the diffusion equation is carried out. Finite differences method is in that case accurate enough to get correct numerical results. A two dimension discretization is done in the Oy and Ox directions. In the case of a 25 nodes discretization, the matrix system equations resulting from the finite difference resolution is: where H i composes locally defined magnetic field and the matrices M, S 1 , S 2 , are defined as follows:

Fractional model
Numerical resolutions (finite differences approximation) lead to precise information (local distribution of the magnetic losses). Unfortunately, if no restrictive assumptions are proposed, whatever the solving numerical method adopted, the non-linear behavior of the system always implies long time calculation and uncertain convergences [13]. To overcome such inconvenient, alternative solution is proposed hereafter. Under weak frequency conditions, quasi-static  [14][15][16], Jiles-Atherton model [17][18][19] or a fractional model [12,20] (quasi-static contribution previously described [12]); provide accurate results of the evolution of the average magnetic induction B versus magnetic excitation field H . Such particular extern conditions mean homogeneous distributions of the induction through the cross section of the sample and consequently homogeneous distributions of the magnetic losses. Unfortunately for those simple models as soon as the quasi-static external conditions expire, huge differences appear. Small improvements can be obtained by adding to this lump model a simple dynamic contribution,product of a damping constant ρ to the time domain derivation of the induction field B. This product is homogeneous to an equivalent excitation field H. Here again, even if this adjunction provides a relative improvement, correct simulation results are obtained on an unfortunately too narrow frequency bandwidth. It seems that a simple viscous losses term ρd P/dt leads to an overestimation in the high frequency part of the magnetization hysteresis loop area versus frequency curve. Another correction of the lump model must be done to reach correct simulation results on a large frequency bandwidth. A mathematical operator dealing with the low frequency and the high frequency component in a different way than a straight time derivative is required. Such operators can be found in the framework of fractional calculus; they are the so-called non entire derivatives or fractional derivatives. Fractional derivation generalizes the concept of derivative to complex and non-integer orders. Fractional time derivative d n B/dt n can be added in our lump model thanks to Grunwald-Letnikov or Riemman-Liouville definitions [21][22][23][24][25]. Both of them are particular cases of a general fractional order operator namely, the first one represents the n order derivative, while the other represents the n fold integral. In this sense, the class of functions described by the Riemman-Liouville definition is broader (function must be integrable) than the one defined by Grunwald and Letnikov. However, for a function from the Grunwald-Letnikov class, both definitions are equivalent. In the present paper, we use the RiemmanLiouville form for n ∈ [0, 1] where is the Euler gamma function. According to the above definition (Eq. 2.9), the fractional derivative of a function f (t) can also be considered as the convolution of a f (t) function and t n / (1 − n). In that case n is the order of the fractional derivation. The additional time derivative present in the formula coincides with the occurrence of positive argument of the gamma function, (.), leading to its convergence to a finite value. It is obvious that fractional derivative includes memory of the previous states. From a spectral point of view, an interesting consequence of fractional derivative consideration is the correction of the frequency spectrum f (ω) of the time domain function f (t) which will multiplied by jω n instead of jω for a first order usual derivation. As a consequence, the fractional derivative provides an interesting alternative able to content the balance requirements, previously discussed between the low and high frequency component of f (ω) [12,26]. Fractional derivative is introduced in our lumped model through a dynamic contribution. We replace the term βdB/dt by ρd n B/dt n , and add this contribution to the quasi-static contribution, Eq. 2.5. Both contributions are included in the last version of the dynamic hysteresis lump model (2.12)

Experimental validation 3.1 Simulation results
Specific experimental set up has been carried out in order to compare both simulation models. An illustration of this experimental set up is given in Fig. 1 The sample tested is a low cobalt iron electrical alloy referenced SV142b. Both thickness and width are 5 mm; its conductivity is 1.4 × 10 7 ( m) −1 . Due to the high electrical conductivity of the tested sample, even for relatively weak variation of excitation frequencies, large differences appear in the losses distribution (Fig. 4). The magnetic excitation amplitude decreases in the radius direction (Ampre theorem), as a consequence a weak non symmetrical distribution of the local magnetic losses is obtained in the same direction. Figure 5a-d show the superimposition of simulated and measured averaged dynamic loops when the surface excitation field is a 1, 50, 200 and 500 Hz frequency sine wave.

Experimental results
All measurments have been done following the international standard on toroidal samples, i.e., with a sinus induction field B imposed of maximum value of 1.3 T. The lumped model simulation is calculated using both contributions. Static contribution provided by a quasi-static model (developed elsewhere), and the fractional model including both wall movement dynamic consideration provide the dynamic contribution. Finite differences model is performed using 10 × 20 nodes discretization. Simulation dynamic parameters are displayed in Table 1. Macroscopic temporal induction field is determined by calculating the average value of the local induction field of each node.
In Fig. 6, we check the accuracy of the fractional lumped model under non sinus type waveform; the magnetic excitation field is composed by a fundamental sinus 50 and 100 Hz of a harmonic 3 sinus excitation. The good correlations between simulations and experimental results give a validation to the approach developed here and confirm that the coupled diffusion equation/wall movement consideration can be described by the lumped fractional derivative model. Note that in the above simulations, the fractional order n = 0.52 gives consistent results for any thickness of the sample with respect to the numerical solution of the corresponding diffusion equation, as well as the experimental data.

Conclusions
In this paper, we succeed in replacing the complex 2D coupling model of magnetic hysteresis and diffusion problem by a lumped model including fractional derivative contribution. Local repartition of ferromagnetic losses depends on the instantaneous material properties and on the frequency of the excitation waveform field. A correct solution to model precisely this repartition implies the resolution in two dimensions of the diffusion equation including local dynamic hysteresis consideration. The resulting model gives precious information but requires complex parameter setting, high computational capacity and long simulation time. Due to the toroidal shape, a 1 dimension solving of the diffusion equation is clearly insufficient and will lead to inaccurate simulation results. A 2 dimensions resolution must be considered and consequently, a large number of discretization nodes is necessary, resulting in prohibitive simulation time and complex memory solicitation. These toroidal shape magnetic cores are mainly industrially employed in current sensors or differential breaker applications. In such application the excitation frequency is clearly unpredictable and high frequency levels are common. Precise model including frequency dependence hysteresis are required. Indeed as soon as the frequency overpassed the usual 50 Hz, the dynamic hysteresis losses increased exponentially. Lumped model including entire fractional derivation provides correct behavior on very large frequency bandwidth; a large improvement is obtained by replacing this entire consideration to fractional one. The fractional lumped model gives precise macroscopic behavior B(H ), it leads to a very accurate formulation of the problem with a high reduction of the complexity and of the simulation times. In future work, our fractional model will be tested to other diffusive and hysteresis phenomenal such as those of ferroelectric materials or other physical domain such as thermodynamic [27].