Eddy-Current Loss Model for Soft Magnetic Composite Materials Considering Particle Size Distribution

Dynamic magnetization curves of soft magnetic materials are often written in terms of magnetic flux density <inline-formula> <tex-math notation="LaTeX">$b$ </tex-math></inline-formula> and magnetic field strength <inline-formula> <tex-math notation="LaTeX">$h$ </tex-math></inline-formula> as <inline-formula> <tex-math notation="LaTeX">$h = \mathcal {H}(b) + c {\text {d} b}/{\text {d} t}$ </tex-math></inline-formula>, where <inline-formula> <tex-math notation="LaTeX">$\mathcal {H}$ </tex-math></inline-formula> is a static magnetization model and <inline-formula> <tex-math notation="LaTeX">$c$ </tex-math></inline-formula> is a real number describing eddy-current effects. In this article, an analytical derivation for <inline-formula> <tex-math notation="LaTeX">$c$ </tex-math></inline-formula> is presented for soft magnetic composite (SMC) materials. The parameter <inline-formula> <tex-math notation="LaTeX">$c$ </tex-math></inline-formula> will depend explicitly on the conductivity of the material particles as well as the geometry of the particles, described by mean particle volume, variance of the particle volumes, volume fraction of the material, and insulation thicknesses. No experimental or empirical parameters appear in <inline-formula> <tex-math notation="LaTeX">$c$ </tex-math></inline-formula> in our treatment.


I. INTRODUCTION
S OFT magnetic composite (SMC) materials consist of electrically insulated compacted ferromagnetic powder [1]. Compared with electrical steels, the particle scale structure of an SMC provides a reduction to classical eddy currents. The magnetic behavior of SMCs is often isotropic. As a drawback, the permeabilities of SMC materials tend to be low. Fig. 1 depicts an SMC material.
There are different approaches for modeling SMC materials. Let us discuss three of them. First, there are models that include some geometric imitation of an SMC particle structure into a computational model, usually based on the finite-element method (FEM). Such studies have been conducted by Maruo and Igarashi [2], Sato et al. [3], Sato and Igarashi [4], and Vesa and Rasilo [5], [6], [7]. In the works of Ren et al. [8], [9], [10], a similar approach in the general level appears. In such computational studies, the particles are usually treated as isotropic and homogeneous continua. The obvious advantage of this paradigm is that classical eddy currents and magnetic behavior, that are affected by the particle structure, can be modeled in a very detailed manner. Furthermore, it is possible to analyze how the particle scale geometry affects the electromagnetic properties of the material. As an obvious disadvantage, we must note that the approach is computationally expensive. Modeling a magnetic component, such as an inductor, with these methods could be heavy in terms of computation time.
Second, there is an approach more suited to modeling actual components. It is to ignore the particle scale structure and replace the ignored scale of the material by a material model. Such a paradigm is efficient when analyzing magnetic components made out of SMC materials. As an example, Guo et al. [11] used this approach to analyze a permanent magnet machine with an SMC core. One might argue that the second approach suffers from the lack of explicit information about the particle scale structures. If the materials were to be improved, the approach provides no insight into the material details that could be tuned. This means that there is a call for yet another approach that combines the best bits of the mentioned ones. In the works of Sato et al. [12], a classical eddy-current loss model for SMC materials was presented. The calculation was based on solving eddy currents in a spherical particle analytically and integrating the loss contributions of several particles defined by some particle diameter distribution. Later, a Cauer network was utilized to convert a complex permeability of a spherical particle for time-domain computations [13].
It seems that there is room for further developments in the approach of treating classical eddy-current effects of SMC particles analytically and combining the results into some macroscale material model. A fully analytical closed-form material model for SMC materials taking into account the geometry of the particles seems to be missing. Laminated cores are often treated using a material model of the form h = H(b) + c db/dt, where H is a static magnetization model and c = σ d 2 /12, where σ and d are the electrical conductivity and the thickness of the laminations, respectively [14]. In this article, we aim at deriving a constant similar to c for SMC materials. But, for SMC materials, the constant c is going to contain the conductivity of the material particles and a description of the geometry of the particles. The geometry of the particles will be described by mean particle volume, variance of the particle volumes, volume fraction of the material, and particle insulation thicknesses. All the dependencies of c are going to be physical or geometric. Fitting parameters that do not have such an interpretation will not be introduced.

II. ANALYSIS OF ONE PARTICLE IN FREQUENCY DOMAIN
In this section, we derive a magnetic model for one magnetically linear material particle in frequency domain with sinusoidal excitations. The derivation is based on assuming the SMC particles to be cylinders, much like the one depicted in Fig. 2. The particles are connected to one another at the ends of the cylinders through insulation layers.

A. Fields in 2-D
Referring to Fig. 2, we denote the radius of the cylindrical particle as R. We treat the magnetic field as rotationally symmetric around the symmetry axis of the cylinder and translationally symmetric in the direction of the symmetry axis. We write the magnetic field h in terms of the radial coordinate r . Next, we investigate how the magnetic flux through the bottom disk of the cylinder and the magnetic field strength on the lateral boundary surface of the cylinder, h s , depend on each other.
Following the treatment presented by Lammeraner and Štafl [15], the equations for the magnetic field strength h are where k 2 = −µσ jω. Furthermore, µ and σ are the permeability and the conductivity of the particle, and j and ω are the imaginary unit and the angular frequency of the excitation, respectively. This is a Bessel-type equation whose solution is given as follows: where J 0 is the zero-order Bessel function of the first kind. Substituting b(r ) = µh(r ) and integrating the flux density over the spherical surface yield the magnetic flux through the surface, which is given by where J 1 is the first-order Bessel function of the first kind. Let us assume that the spherical cross section of the cylindrical particle is surrounded by some insulating material. We assume that the effective cross-sectional area that the magnetic flux passes through is given by π R 2 /η, where 0 < η ≤ 1 is a coefficient that models the volume fraction of the material. The effective flux density is now given by in the frequency domain.

B. Reluctance Network in 1-D
Let us now connect the particle into another particle through a thin insulating layer. We use well-known methods of reluctance networks.
Referring to Fig. 2, the darker region on the right-hand side of the cylinder represents a thin insulation. The reluctance of the insulating layer is given by where µ ins is the permeability of the insulation, l = 2Rκ is the length of the particle, where κ describes the length relative to the diameter 2R of the particle, and τl is the thickness of the insulation, where the constant 0 < τ < < 1 determines the relative thickness of the insulation. Assuming that the flux through the particle passes also the insulation, the magnetic field strength in the insulating region is related to the effective flux density of the particle by The volume fraction coefficient η appears in both (4) and (6). The reason for such a convention is that all the flux is assumed to pass the cylindrical particle and the associated insulation, but there are insulating surroundings outside the lateral surfaces of the particle and the insulation that are taken into account by the coefficient η.
The effective magnetic field strength h eff,l over the particle of length l and the insulation of thickness τl is given by Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. Combining (4), (6), and (7) yields where the function g is farsightedly characterized by the argument k R = R(−µσ jω) 1/2 .

C. Taylor Expansion
We now aim at turning the single particle model (8) into time domain. We also want to perform an integration over a distribution of particle volumes. Let us approximate the function g as its N th-order Taylor expansion, and let us substitute the volume of the particle V = 2κπ R 3 . We get (9) where only the even coefficients γ 2i = 1, −(1/8), . . appear, as explained in the literature [16]. Since only even exponents are nonzero in the expansion, formula (9) expresses the magnetic field strength in terms of rising integer exponents of jω, which is straightforward to transform into time domain. Furthermore, the formula is expressed in terms of exponents of the volume V , which provides a neat access to integrations over V .

III. MULTIPARTICLE MODEL
Let us assume that the particle volumes of the material more or less obey the gamma distribution, given by the density function where the parameters α and β are related to the expected value E and variance Var by the formulas Roughly speaking, the density function expresses "how many particles have the volume V ." The function Γ is given by Referring to Fig. 3, we assume that we have a sequence of particles whose volumes obey the gamma distribution. Each particle has the same effective magnetic flux density b eff . Fig. 3 illustrates the consequences of the particles having the same effective magnetic flux density. If a larger particle is connected to a smaller particle, some flux diverges from the sequence to the surroundings, and if a smaller particle is connected to a larger particle, some flux is introduced to the larger particle from the surroundings. Even though skin effect is present in model (8), we assume that the magnetic coupling between the particles by equating b eff is not affected by the flux density distribution in individual particles.
The effective magnetic field strength over the sequence of particles is obtained by computing the total magnetomotive force over the length of the sequence. We get where the multiset S contains all the lengths of the particles in the sequence. Next, we aim to approximate the sums of (13) by integrals involving the gamma distribution (10). Since (13) involves the lengths of the particles, but the density function (10) is expressed in terms of particle volumes, it is worth denoting the length l of a particle in relation to the volume V of the particle. Given that V = π R 2 l and l = 2Rκ, we have Taking the leap from (13) to the continuous domain, we get Substituting (9) into (15), we get where the rather technical computations of the integrals are available in the Appendix. We manipulate (16) further by substituting (11) and replacing jω by time differentiation. We get (17) where inside the big parentheses the static contribution of the insulations, the static contribution of the particles, the firstorder dynamic contribution of the particles, and the higher order dynamic contributions of the particles have been separated. Truncation of (17) from an appropriate order, usually one, yields a simple dynamic hysteresis model. The material model (17) is already in the desired form, but the model is linear and scalar. Generalization of the model into 3-D is straightforward, provided that linearity holds. Due to superposition, the scalar model is applicable componentwise. Furthermore, truncating (17) from the first-order dynamic term, we get where the magnetic model H is assumed to be linear through the relations H ins (B) = 1/µ ins B and H part (B) = 1/µB. The novelty of this article is the derivation of the coefficient c ed . The decomposition of H is just a consequence of the wellknown reluctance network approach. A summary of the symbols and their possible numerical values is given later in Table I. Whether (18) is applicable for nonlinear static models H, it is beyond the scope of this article, but it seems that for laminated cores, the corresponding equation is valid in the nonlinear case [14]. This is based on the assumption that skin effect of the magnetic flux density is negligible in the sheets, which is valid when the excitation frequencies are relatively low. Then, the contributions of the nonlinear static model and the dynamic part may be separated by a simple superposition. If the same argument was applied to (18), it would mean that the operators H, H ins , and H part could be nonlinear static magnetic models.

IV. NUMERICAL EXAMPLE
In this section, we provide an experimental study to demonstrate the applicability and the limitations of the model. Magnetic properties of a sample are studied over a wide frequency range with low-amplitude magnetic excitations.

A. Measurements
The SMC material depicted in Fig. 1 serves as a case study sample. Multiple microscope images, much like the one appearing in Fig. 1, were taken, and the particle surface areas were computed from the images. The computed particle surface areas were raised to the power of 3/2 to estimate the volumes of the particles. Fig. 4(a) depicts the measured distribution of the particle volumes. Integration over an interval of volumes tells us how many particles belong to that volume range relative to the total number of the particles. Fig. 4(b) is explained later.
The SMC material was characterized using a two-coil setup with an LC R meter. Frequency-dependent complex permeability µ meas was computed from the measured impedance between the sinusoidal low amplitude primary current and the sinusoidal low amplitude secondary voltage. Root-meansquared (rms) secondary voltage was set to a fixed value V 2,rms = 1 V. By the Faraday law, the rms magnetic flux density was b rms = V 2,rms /(2π A f ) ≈ 1500 T · Hz/ f . However, at lower frequencies, the current limit 20-mA rms of the device is expected to limit the flux density. The Ampère law yields a bound b rms = |µ meas |I 1,rms /l ≈ 0.2 A/m · |µ meas |. Measurements were repeated for multiple values of the secondary voltage within a range of V 2,rms = 0.01, . . . , 1 V to ensure that the frequency-dependent µ meas is independent of the amplitude of the flux density. This means that the measurements were carried out in the magnetically linear region of the material.
where w meas corresponds to energy loss density over one cycle of excitation. The terms of (19) are arranged in such a way that the right-hand side contains no information other than what can be determined by the measured complex effective permeability. By linearity, the loss density scales to the second Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. power of the flux density, and hence, the left-hand side is in a normalized form. It should be noted that (19) accounts for all the normalized losses that have been written into µ meas . In case of measurements, this includes all the losses that are visible to the LC R meter through the measured impedances. The effective frequency-dependent complex permeabilities obtained from our computational models contain information about the eddycurrent losses.

B. Models
We compare two models with the measurements. Both models provide a relation between complex b and complex h and, thus, provide a complex effective permeability. Hence, also, losses are available by (19).
On the one hand, we use model (18) in its complex scalar form, and we call this model "Computed." On the other hand, a model with less approximations is obtained by integrating (8) over the estimated gamma distribution numerically by formula (15). We call this model "Computed, Bessel," since it involves numerical evaluations of Bessel functions. There are two reasons for the use of "Computed, Bessel" in the comparison. First, our model (18) was derived from the truncated Taylor polynomial of (8), whereas the "Computed, Bessel" uses (8) directly. Hence, we can take a look at how the truncation of the Taylor polynomial affects the numerical results. Second, we need to identify local material parameters from macroscale measurements, and we will see that the "Computed, Bessel" model provides a good forward model for the estimation problem.

C. Identifications
Model (18) involves quite many parameters that should be identified from the material depicted in Fig. 1. Table I provides a summary of all of them. Next, we discuss how to provide reasoning for each of the identified parameters.
In Section IV-B, we discussed how the measured particle volume distribution was obtained for Fig. 4(a). The expected particle volume E and the variance of the particle volumes Var are just the expected value and the variance of the measured distribution. Numerical values of these parameters are available in Table I. The dashed red curve "Estimated" in Fig. 4(a) represents the gamma distribution (10) that is assumed to approximate the measured distribution. The gamma distribution is fully determined by E and Var. Fig. 4(a) does not really reveal the weaknesses of the use of the gamma distribution. Hence, Fig. 4(b) is provided. It contains transformed versions of the distributions that appear in Fig. 4(a). The transformations include multiplications by particle volumes, after which normalizations were carried out. Finally, transformation of the x-axis was conducted. From Fig. 4(b), it can be seen that the measured distribution is more skewed toward smaller particle sizes than the estimated distribution even though expected values and variances match. This is a sign that the actual particle distribution does not exactly follow the gamma distribution. But, in the end, some reasonably realistic distribution has to be used for an analytical derivation of the eddy-current coefficient c ed .
The permeability of the insulations µ ins was chosen to be the permeability of free space µ 0 for non-magnetic behavior. The conductivity σ and the permeability µ of the particles as well as the insulation thickness factor τ of the material cannot be directly measured from the material sample, but they have to be estimated from the macroscale magnetic behavior of the sample. Let us denote the computed effective frequencydependent complex permeability obtained using the model "Computed, Bessel," by µ B . We define an objective function where the summation is defined over the set F containing all the measurement frequencies. The parameters µ, σ , and τ were found by minimizing the objective function e under the constraint τ = |µ meas,dc |/µ − η η − |µ meas,dc |/µ ins where µ meas,dc is the measured permeability at the dc limit. This constraint is to ensure that the measured and computed permeabilities match at the dc limit. The concrete computations were carried out using the reflective trust region optimization algorithm that is available in the Python Scipy library [17]. Numerical values for µ, σ , and τ are available in Table I. It was assumed that the particles are more or less round, and hence, the relative length κ of the particles was chosen to be 1. Furthermore, the volume fraction η of the material was set to be 90%. This is close to the values discussed earlier in the literature [5], [18].

D. Numerical Results
We use the two models, "Computed, Bessel" and "Computed," with the terminology summarized in Section IV-B. The models use the parameters that are summarized in Table I. In Fig. 5, we see that the measured permeability curve agrees well with the "Computed, Bessel" model. This is due to two factors. First, model (8) that is integrated numerically over the particle volume distribution is based on reasonable assumptions. It seems that classical eddy currents, including skin effect, explain the behavior of the material for such a frequency sweep well. Second, since not all of the parameters listed in Table I can be measured directly from the material sample, an optimization algorithm was utilized to find some of the parameters. This was discussed in Section IV-C. Such fitting may also compensate some inaccuracies related to, for example, the estimation of the material particle volume distribution.
From the losses depicted in Fig. 6, we can state that the overall shape of the measured losses agree with the losses obtained using the "Computed, Bessel" model. There is a 4.4% relative difference of the losses at 1 MHz and a 19% relative difference at 10 MHz. It is unclear why there is this discrepancy. On one hand, we could consider if the used model contains some unnecessarily restrictive assumptions. One assumption we made is that the particles in the "Computed, Bessel" model are not capacitively coupled, and only eddy currents inside the particles contribute to the losses. However, the measurements contain all the losses that are visible to the LC R meter by the measured impedances. Considering capacitively conducting currents in the model could yield an increase for the losses in the high-frequency range. On the other hand, it may very well be that the optimization task that is based on minimizing (20) does not give enough weight for the imaginary part of the "Computed, Bessel" model, and hence, there is a discrepancy in the high-frequency losses.
The model "Computed" is based on a first-order Taylor expansion of the more complicated model. Hence, the agreement between the model "Computed" and the measurements is reasonable only in the low-frequency region, where the Taylor polynomial was developed. The permeability trend in Fig. 5 is similar to the measured, but no overall quantitative agreement is expected. At 400 kHz, the relative difference between the measured permeability and the "Computed" permeability is 9.6%. The losses of the "Computed" model capture only the linear part of the loss curve. At 400 kHz, the relative difference of the measured and the "Computed" losses is 7.8%. We conclude that the truncated model (18) is applicable only in the low-frequency region.

V. CONCLUSION
The contribution of this article was to derive the parameter c ed appearing in the dynamic hysteresis model (18) for SMC materials. The coefficient c ed depends explicitly on the geometric and physical parameters of the material, and no additional experimental coefficients were brought into the consideration. A list of the physical parameters with numerical values for one specific material was given in Table I.
It was shown experimentally that the underlying more complicated eddy-current model for SMC materials captures the behavior of the case study material in the sense of permeability and losses. The eddy-current coefficient c ed was derived from the first-order Taylor coefficient of a more complicated model. Hence, the simplified dynamic hysteresis model (18) with the coefficient c ed should only be used for low-frequency applications. For the case study SMC material, it was found that permeabilities and losses were at 10% relative vicinity of the measured quantities when the excitation frequency was less than 400 kHz. This 400 kHz can be taken as "low frequency." APPENDIX In Section III, the derivation of (16) was omitted. Here, we provide the missing pieces.
First, we show one integral formula. Let α, β > 0 and q ≥ 0. Using the substitution V = β −1 W between the variables V and W , we have by definition (12) of the function Γ . Second, we compute a fraction of integrals farsightedly. Let i ≥ 0. We use definitions (10) and (14) for the gamma distribution function G and the particle length l. We notice that Substituting (23) and using formula (22), we find the integral formula R + V Third, we substitute (9) into (15) and use linearity of integration. Computing the integral with (24), we have which is equal to (16).