Low computational cost semi-analytical magnetostatic 1 model for magnetocaloric refrigeration systems 2 3

The analysis of the active magnetic refrigeration (AMR) cycle for different wave9 forms of both the magnetic field and the velocity of the heat transfer fluid is an essential chal10 lenge in designing and implementing heating and cooling systems based on the magnetocaloric 11 effect. One of the most important issue is the correct modelling of the magnetic and thermal 12 behavior of the active magnetocaloric materials (MCM) in order to estimate precisely cooling 13 capacity of the magnetocaloric system. As the multiphysics coupling implies successive calls 14 for both the thermal and the magnetic modelling subroutines, the execution time of these sub15 routines has to be as short as possible. 16 For this purpose, a new magnetostatic model based on reluctance network has been performed 17 to calculate the internal magnetic field and the internal magnetic flux density of the active mag18 netocaloric material (gadolinium, Gd) inside the air gap of the magnetic circuit. Compared to a 19 3D Finite Element Model (FEM), our magnetostatic semi-analytical model leads to a sharp drop 20 of the computation time, while offering a similar precision for all magnetic quantities in the 21 whole magnetocaloric system. 22


Introduction
The magnetic refrigeration and heat pumping is an emerging technology which offers environmental benefits compared to conventional vapor compression machines.
The magnetic refrigeration is based on the magnetocaloric effect (MCE) exhibited by some materials at room temperature. 1The magnetocaloric effect occurs during the critical transition paramagnetic/ferromagnetic of these ferromagnetic materials: any change in external magnetic field around the Curie temperature induces a reversible change in the correlated electronic spin entropy, directly related to strong specific thermal power density production/absorption.When these magnetization changes occur in an adiabatic way, they produce an adiabatic temperature change ΔTad, which is a characteristic property of the magnetic material and depends on both the magnetic field evolution and the initial temperature.Since the best ΔTad is in the range of only a few kelvins per tesla, an active magnetic refrigeration (AMR) cycle has to be imposed to magnetocaloric regenerators -micro-heat exchangers composed of magnetocaloric plates or spheres, etc.-to produce larger temperature gradients and significant thermal power for heating or cooling purposes.An AMR cycle imposes a fluid (coolant) to flow alternatively through the regenerator while synchronizing the magnetization-demagnetization of the regenerator. 2 several years, the Energy department of the FEMTO-ST Institute has been developing research on high efficiency magnetocaloric devices using magnetocaloric properties of some materials around ambient temperature. 3This technology is meant to provide large-scale ecological solutions for refrigeration and heat pumping, since the theoretical efficiency of the AMR cycles is much higher than conventional vapor compression technologies.On the other hand, its operation does not require any use of greenhouse gases, unlike conventional refrigeration machines.
However, the design of such efficient magnetocaloric devices still requires further research for faster precise simulation codes and multiphysics models.This paper aims at exposing the theoretical basis and numerical results of a new model based on reluctance network applied to the Magnetic Equivalent Circuit of our magnetocaloric device, and showing the large improvements provided in computational time for the simulation of large magnetocaloric systems, while offering a good compromise on its accuracy.Even if models using reluctance network have been widely developed in electrical machine modeling, there are very few reluctance network model in the magnetic refrigeration domain.
Dai and al. 4 evaluate the influence of the gadolinium stack on the magnetic flux inside the air gaps of a magnetic circuit including permanent magnet, using an average permeability of the gap based on a simplified reluctance network, with an averaged magnetic flux density in the air gap.The magnetic field inside the gadolinium plates and the permeability of gadolinium are simply calculated by the means of the average field Brillouin function.The authors show that the change in temperature of the gadolinium stack strongly modifies the magnetic field in the gap, so that iso-field phases cannot be used for an accurate modeling of Ericsson cycles.
However, nonlinear behavior of the ferromagnetic circuits is not considered.
More recently, Vuarnoz and al. 5 express the overall reluctance of the air gap as a seriesparallel assembly of mean reluctances of the magnetocaloric plates, fluid layers and air gaps, allowing for the calculation of the mean magnetic flux and the corresponding magnetic flux density inside the plates.Using this model to calculate a simplified Tad-based magnetocaloric heat-source term in a simple 1D thermo-fluidic model of an AMR cycle with imposed external magnetic field, the authors find rather underpredicted values of the plates and fluid temperatures, however with less discrepancy compared to the results of simple uniform magnetic field models.Their calculations as well as their experimental results confirm the inability of permanent magnets to impose a constant internal magnetic field during cold or hot blows in AMR cycles, as shown previously in Ref. 4.    In this paper, a new magnetostatic model based on reluctance network is applied to the global magnetocaloric system in order to calculate the internal magnetic field and the internal magnetic flux density of the active magnetocaloric material (gadolinium, Gd).The model takes into account the non-linear behavior and non-homogeneity of the ferromagnetic materials, while highly reducing the computation time compared to a previous FEM model.

Experimental prototype
In the magnetic refrigeration domain, different prototypes based on active magnetic regenerative refrigeration (AMRR) principles were built and tested for more than twenty years. 6pecific experimental prototype has been developed and created in our laboratory (Fig. 1), including a controlled power source for pulsed magnetic field (powerful electromagnet) and a controlled hydraulic cylinder specially designed to produce precise flow sequences through an active magnetocaloric regenerator block inserted between two micro-heat exchangers at both ends.The regenerator consists of 14 pure gadolinium rectangular parallel plates -12 central plates (13 x 1 x 45 mm 3 ) and 2 external plates (13 x 0.5 x 45 mm 3 ), all 0.5 mm apart from each other (Fig. 2) -; it is arranged inside the electromagnet 21 mm wide air-gap.The purpose of the bench is to characterize the thermofluidic behavior of the magnetocaloric regenerator and to optimize its refrigeration performances. 7,8IG.1: Femto-ST experimental prototype with gauss-meter to measure magnetic flux density FIG.2: Experimental magnetocaloric regenerator module

Magnetostatic modelling of the magnetocaloric system
Many numerical models of AMR cycles have been developed in the last ten years.Among the most recent models, some consider the magnetic field equal to the applied field; 9 some others consider the magnetic flux density to be uniform in the whole material. 10More precise models simply evaluate the internal magnetic field at each point of the regenerator by using averaged demagnetizing factors. 11n our previous work, a magnetostatic finite element model (FEM) was developed and implemented in a multiphysics analysis of the refrigeration test bench, and a significant magnetic field heterogeneity was observed. 12This heterogeneity can be explained by the strong dependence of the internal magnetization of the magnetocaloric material on the temperature close to its Curie point combined with the temperature gradient between cold and hot sides of the magnetocaloric regenerator.Another inhomogeneity factor is produced by 3D phenomena in the magnetic field distribution inside the air-gap of the magnetocaloric test bench. 13ever, even if the magnetostatic FEM model offers a very accurate resolution of Maxwell equations, its main drawback is the highly time-consuming calculation to solve these equations in the whole magnetic domain.The computation time of the magnetostatic model is crucial, since it has to be used in each loop of the iterative calculation process until convergence of the multiphysics analysis, which requires a huge number of calls of the magnetostatic FEM model for each magnetocaloric operating point. 8this work, a semi-analytical model is designed for both the magnetic flux density and the magnetic field calculation inside the whole magnetic system, allowing much faster calculation than the FEM model with comparable precision, and accounting for the heterogeneous distribution of the magnetic field in the AMR cycle.
Since the main objective of our test bench is to precisely characterize the AMR cycle, it is important to estimate the magnitude of the magnetic field distribution inside the regenerator, accounting for the non-linear behavior of both the external magnetic circuit and the active magnetic regenerator. 14

Magnetic device geometry and AMR discretization
The electromagnet working was first simulated with FEM.The ferromagnetic circuit of the experimental prototype has an "8-shaped" structure presented in Fig. 3.It is composed of soft Fe-Si sheets, and a cylindric hole is arranged along the center line to enable further optical measurements of the fluid velocity (micro-PIV measurements).The pulse-controlled magnetic field is produced by 4 large coils with 90 turns each (colored rectangles on the side view, Fig. 3).

FEM model
The internal magnetic flux density B and internal magnetic field H of each segment is calculated for electrical coil currents ranging from 0 to 50 A with 5 A steps.The results obtained with a 3D FEM using Flux© magnetostatic solver can be considered as a reference, since this computation code can take complex magnetic phenomena into account, such as 3D effects, nonlinear B(H) behavior, magnetic flux losses, demagnetization, etc.In particular, the demagnetizing field has a strong influence on the internal magnetic field, on the magnetic flux density and on the magnetization inside the AMR. 5 Therefore, oversimplified hypothesis on demagnetizing factors should be avoided when calculating the internal magnetic field in AMR materials.
FEM models lead to a very precise magnetic analysis but are highly time-consuming (about five minutes for each magnetostatic resolution point of our experimental prototype with a 2.80 GHz dual-core PC).Even if the FEM model is accurate, the computation cost seems to be prohibitive when this model is coupled to a multiphysics simulation loop.In this case indeed, a huge number of successive estimation points are required to insure convergence, leading to a very large computation time.The numerical results of FEM calculations are exhibited below and compared to those obtained with our semi-analytical model in section 4.

Semi-analytical model based on reluctance network
In the present work, a semi-analytical model has been designed to accurately model the distribution of the internal magnetic field at each point of the regenerator for each time step during the AMR cycle, allowing much faster calculation than FEM with comparable precision, while accounting for the heterogeneous distribution of the magnetic field in the AMR.

Mapping the non-linear behavior of the external ferromagnetic yoke
The first feature of our semi-analytical model is to take into account the non-linear behavior of the external ferromagnetic yoke of the test bench, by means of experimental measures of the magnetic flux density in the empty air-gap.The calculation scheme is displayed in Fig.

FIG. 5: Calculation scheme of mmf for the external ferromagnetic yoke with empty air-gap
During this step, the magnetomotive force (mmf) of the ferromagnetic circuit UmFe is calculated.The magnetic flux density of the central air gap B0 is precisely measured in the core of the empty air-gap with a gauss-meter (as shown in 0), for currents ranging from 0 to 55 A with a 5 A step.This leads to the magnetic flux in the central air gap 0 [Wb] depending on the magnetic flux density and the associated air-gap section S0 = 90 × 50 mm 2 (red part in Fig. 6 and Fig. 7): From there, the reluctance of the whole air-gap ( 2) is calculated considering the leakage flux between magnetic poles outside the central rectangular air-gap (external trapezoidal air gaps): where Lj and Sj are the length and section of the j th among n = 50 subdivisions of the symmetric trapezoidal leakage flux zone (an example of segmentation with n = 10 is presented in blue in Fig. 6).The magnetic flux density Bj in the externals trapezoidal air gap was measured and linearly estimated along the n = 50 subdivisions (the values for j = 1, n/2, n are displayed in Fig. 7).The magnetic flux in the whole empty air gap tot [Wb] expresses then as: In these conditions, the mmf values Um0 (equation 4) of the empty air gap are obtained using this experimental method and compared with those calculated by the FEM, leading to a relative difference smaller than 0.4%: An equivalent diagram of the magnetic circuit with an empty air gap is displayed in Fig. 8 showing the relation between the global magnetomotive force equal to NI, the magnetomotive force of the ferromagnetic circuit and the magnetomotive force of the empty air gap.

FIG. 8: Equivalent diagram for the calculation of the ferromagnetic circuit mmf
Therefore, the magnetomotive force UmFe across the ferromagnetic circuit can be deduced from the equivalent Kirchhoff voltage law for magnetic circuit (5), where I is the intensity of the electric current flowing in the coils and N is the number of coil turns (Fig. 9).
FIG. 9: Obtaining the magnetomotive force of the ferromagnetic circuit

Computing magnetic flux distribution in the AMR regenerator
In the second step, we analyse the magnetic flux distribution in the AMR regenerator plates for an inhomogeneous temperature distribution.The magnetic flux density distribution Bi,k for each element ei,k is computed in the discretized regenerator, for a given coil current at a given temperature distribution Ti,k in the regenerator.This step considers the magnetic behavior of the ferromagnetic yoke UmFe(Φtot) obtained previously and represented in Fig. 9. bounded by two R1 reluctances in series with its non-linear UmGd i,k(Φi,k) (as shown in Fig. 11).
Since the temperature Ti,k of each element ei,k is known, the curve Bi,k(Hi,k) is determined by the approximation of the experimental MCM characteristic curves B(H, T) for T = Ti,k.This allows to calculate the spatial distribution of the characteristic UmGd i,k(i,k), using ( 6) and (7), where Si,k is the element section and l the width of the element: The mmf UmGd i,k of each element in series with two R1 reluctances leads to the global mmf Um0 in the air gap, as a function of the flux Φi,k through the element ei,k :

FIG. 11: 3D view of the air-gap with some gadolinium plates
In order to identify the magnetic flux distribution in the air spacing of the air-gap, it is necessary to compute the reluctance Rair with (9): where Rk is the reluctance of the parallelipedic air space crossed by the magnetic flux flowing through the cross-section Sk of the gadolinium k plate.
The magnetic flux of the air spacing inside the air-gap is calculated by (10), considering the linear dependence of the magnetic flux in the air gap: This leads to the mapping function Φtot(Um0) of the global air gap flux noted tot function of the total mmf of the gap including the regenerator (Um0).For a given magnetomotive force Um0, the global magnetic flux is obtained by adding all element fluxes i,k with the magnetic flux air in the air gap (11): What remains to be done is to compute the cartography (tot) of the total magnetomotive force corresponding to a total magnetic flux, by adding the magnetomotive force in the global air gap, using the inverse interpolation of Φtot(Um0) obtained with (11).Thus, (N.I) is obtained by interpolation of Θ for a regular dataset of the magnetic flux between 0 and a maximum value, according to Kirchhoff's law for magnetic circuit (12): The proposed algorithm contains several sub-functions considering the nonlinearities in active magnetic materials (gadolinium plates, ferromagnetic yoke).In Table .I, we observe the complete scheme of the algorithm.The different steps are exposed in detail below, and two preliminary computations are necessary: -external ferromagnetic circuit behavior tot(UmFe) -Rair and R1 air reluctances

Graphical interpretation of the algorithm
Thereafter, in these conditions, a series of interpolations have to be performed in order to obtain the internal Bi,k and Hi,k values of each element of the magnetocaloric material: -the global magnetic flux through the gap Φtot is calculated from the current I in the coils (consequently NI) and the mapping tot(NI).Thus, the first curve fit method can be performed (step ⑥ in Table .I and interpolation ❶ in Fig. 12); -the second interpolation is performed to obtain the mmf of the gap Um0, which depends on the curve tot(Um0), using the values of the magnetic flux in the gap obtained with the previous interpolation (step ⑦ in Table .I and interpolation ❷ in Fig. 12); -the next interpolation aims at obtaining the magnetic flux in each segment using the mmf calculated previously.For this purpose, it is necessary to use the curve connecting the magnetic flux of the segment and the mmf of the gap since the dependence i,k(Um0) is known for each element ei,k (step ⑧ in Table .I and interpolation ❸ in Fig. 12); -the last curve fit method is achieved to obtain the mmf of each element, which depends on the curve i,k(UmGd i,k), using the values of the magnetic flux of the corresponding segment obtained with the previous interpolation (step ⑨ in Table .I and interpolation ❹ in Fig. 12).

FIG. 12: Successive interpolations for obtaining the magnetic flux and the magnetomotive force in each segment
The ultimate step is to get the internal magnetic flux density Bi,k and internal magnetic field Hi,k of each element, by dividing i,k and UmGd i,k respectively by the section Si,k and the length l (see example in Fig. 13).So, the characteristics Bi,k (Hi,k) of the magnetocaloric material (gadolinium in our case) are obtained after achievement of all the steps with our code developed in Python language.FIG.13: Obtaining the internal magnetic flux density and magnetic field of an element ei,k

Results and comparisons
The results obtained by the FEM (Flux3D©) for the magnetic flux density calculation in a ten-fold segmented plate of gadolinium at imposed uniform temperature T = 293 K are compared to the results of our semi-analytical model.The outputs for the central and the external elements of the central plate are presented in Fig. 14 and Fig. 15, respectively (the corresponding segments are highlighted in blue in Fig. 4).The maximum difference between these results is lower than 2%.The difference between the magnetic flux density curves of the two elements shows how much the magnetic flux density is depending on the position of the element inside the air gap, even if the temperature is uniform.Indeed, the external lines of the magnetic flux density concentrate at the ends of the regenerator plates, so that the magnetic flux density is higher (0.05 T) at the ends than at the center of the plates.More thoroughly, it has been shown in Ref. 13 that strong inhomogeneities in the magnetic field inside an empty air gap are induced when introducing a magnetocaloric plate into an applied magnetic field, leading to a higher magnetic field intensity in the air, along with a large drop of the latter inside the plate very close to the ends.This is due to both the magnetic permeability of gadolinium (µr ~ 1.5) and the demagnetizing field on the boundary surfaces, which increases the axial demagnetization factor close to the ends of the regenerator plate.Since the magnetic flux is conservative, the decrease of the magnetic field intensity at each end of the magnetocaloric plate leads to a shift of the magnetic curve toward lower values of H, as shown in Fig. 15.Moreover, both the magnetic flux density and the magnetic field depend strongly on temperature, more particularly around the Curie temperature (Tc  293 K for gadolinium).The comparison will further focus on the magnetic flux distribution in the central gadolinium plate, with imposed temperatures at each end (left cold side 290 K; right hot side 299 K).For this purpose, Fig. 16 shows the magnetic flux density distributions in the central gadolinium plate that have been obtained by the FEM and semi-analytical models respectively, with 50 A current intensity in the coils producing 1 T magnetic flux density in the empty air gap.

FIG. 16: Magnetic flux density distribution comparison
The accuracy of the semi-analytical model can be highlighted in the case of a linear temperature gradient from 290K to 299K imposed between the cold and the hot sides of the regenerator.The dependency of the internal magnetic flux density on the internal magnetic field in the gadolinium central plate divided into 10 equal elements is calculated with our semi-analytical model and with Flux3D © and displayed on Fig. 17   As can be seen in both figures, the magnetic flux density is much higher at the cold side of the plate, because of a higher magnetization due to the ferromagnetic state of the gadolinium under its Curie temperature, while being much lower at the hot side above the Curie temperature (paramagnetic state).
No significant difference can be detected between the two figures, which confirms the good accuracy of our semi-analytical model.No discrepancy higher than 3% was obtained between these two simulations.
As a final comparison, the two different methods for obtaining internal magnetic field and magnetic flux density of the magnetocaloric material exhibit the following advantages and drawbacks: -the finite element method (FEM) leads to very precise magnetic analysis (3D phenomena) but is highly time-consuming when using a fine spatial resolution grid.For instance, the calculation time needed for achieving the FEM numerical simulation applied to the evolution of internal magnetic field and magnetic flux density in the whole magnetic system with a 5A step [0 : 50, 5A] is around 35 min with a 2.80 GHz dual-core PC; -the semi-analytical model always requires much less time than the FEM resolution and offers good analysis precision with very similar results.The calculation time needed for the same analysis under same conditions drops to around 5 s with the same computer.

Conclusions
A semi-analytical model based on a reluctance network is proposed to calculate the internal magnetic field H and internal magnetic flux density B when introducing a magnetocaloric stack of parallel plates into a magnetic air gap.The computation time needed with our model for obtaining internal magnetic field H and magnetic flux density B at each point of the whole magnetic system is around four hundred time shorter, compared to a complete FEM solver resolution for the same system in the same conditions.
The semi-analytical model is directly applicable and can easily be extended to rotary magnetocaloric devices.A multiphysics analysis will be described in a further paper, in which the presented magnetostatic model is combined with a thermo-fluidic model, allowing to achieve the simulation of 5 AMR cycles applied to the whole system within only 5 minutes, instead of 4 hours with 3D FEM 12 (using a 2.80 GHz dual-core PC).Besides, the magnetic results obtained with our semi-analytical model are very close to those of the FEM solver Flux3D © , which proves the efficiency of our model.Moreover, our semi-analytical model was applied efficiently to an optimization process of an AMR refrigeration system.The model involves the multi-physics phenomena occurring during the magnetic refrigeration process and simulates both the magnetocaloric material behavior, the instantaneous heat exchanges and the evolution of the whole magnetocaloric system undergoing successive AMR cycles with a controlled periodic applied magnetic field.This allows to simulate the experimental behavior of our magnetocaloric test bench developed at the FEMTO-ST Institute.
More generally, our analytical reluctance-based model could easily be adapted to any other magnetocaloric regenerator device to save huge calculation time without losing precision and to allow easier multi-parametric optimization.

Fig. 10 FIG. 10 :
Fig. 10 displays the magnetic equivalent diagram of the air-gap with gadolinium plates.

①⑦
Inputs: Ti,k, I ② Mapping Φi,k(Um0), equations(6) to(8)   a) Obtaining magnetic behavior interpolation functions B(H)i,k for each element ei,k.b) Transforming B(H)i,k into a flux mmf interpolation function Φi,k(UmGd i,k) -Eq (6) and (7) c) Transforming UmGd variable in Um0 -Eq (8) ③ Considering the linear dependence of the magnetic flux Φair(Um0) in the air gap -Eq (10) ④ Mapping function of the global air gap flux Φtot(Um0) -Eq (11) ⑤ Computing the cartography Θ(Φ) -Eq (12) ⑥ Computing the total mmf Φtot for a given current I 16 Obtaining the gap mmf Um0 by interpolation of Φtot(Um0) ⑧ Obtaining the magnetic flux Φi,k through each element by interpolation of Φi,k(Um0) ⑨ Obtaining the mmf UmGd i,k of each element with the dependence Φi,k(UmGd i,k) ⑩ Obtaining the internal magnetic flux density Bi,k and internal magnetic field Hi,k of each element

FIG. 17 :
FIG. 17: B(H) curves at each point of the regenerator at corresponding temperature (semi-analytical)

TABLE . I
: Complete scheme of the air-gap with gadolinium plates and equivalences