Detailed numerical simulations of a single stage of rotatory active magnetic regenerators: Influence of the pin geometry

International Journal of Thermal Sciences 149 (2020) 106198

Contents lists available at ScienceDirect

International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts

Detailed numerical simulations of a single stage of rotatory active magnetic regenerators: Influence of the pin geometry Ibai Mugica a ,∗, Sébastien Poncet a , Jonathan Bouchard b a b

Mechanical Engineering Department, Université de Sherbrooke, 2500 boulevard de l’Université, J1K 2R1, Sherbrooke, Canada IREQ-Laboratoire des Technologies de l’Énergie (LTE), 600 avenue de la Montagne, G9N 7N5, Shawinigan, Canada

ARTICLE

INFO

Keywords: Magnetic refrigeration 3D Numerical simulation Active magnetocaloric regenerator Demagnetization Tortuosity Multiphysics

ABSTRACT The following paper examines different configurations of pin geometries for rotatory magnetocaloric regenerators, through detailed multiphysics simulations. The analysis is first done by characterizing the performance of each geometry against the multiple physical phenomena that intervene in a magnetocaloric regenerator: magnetic, fluid and thermal fields. Then, full 3D AMR cycles are simulated for a selected group of pin geometries, and conclusions are drawn between the temperature span of the regenerators and the separately analyzed physical phenomena. The main novelty of this study is the analysis of the tortuosity of the fluid flow, a characteristic parameter of porous media. One of the main conclusions is that an increased tortuosity might be a valid alternative strategy to increase the power of magnetocaloric regenerators.

1. Introduction Current commercial refrigeration systems exploit a liquid–gas phase transition to transport heat from one place to another (by latent heat more specifically). These ‘‘vapor-compression’’ cycles have thus far provided the refrigeration needs for the general population worldwide. Within the research field of refrigeration, new alternative technologies seek to overcome the drawbacks of current vapor-compression technology by exploiting a different thermodynamic phase transition. The use of alternative refrigeration technologies could result in a reduction in polluting refrigerants (halogen fluid blends eg. HFCs, CFCs, HCFCs), reduced noise levels due to the absence of compression/expansion stages, and an increase in energy efficiency due to a more isentropic phase transition amongst other outcomes. One of the most researched and promising alternative refrigeration technologies is magnetocaloric refrigeration (see [1]). Magnetocaloric cycles use the magnetic phase transition of solidstate refrigerants to transport heat. Much like the liquid–gas phase transition, the ferro-paramagnetic phase transition provides a sudden change of the thermodynamic properties of the refrigerant (see [2]), which can be exploited into a cycle to transport heat. In the vicinity of its transition temperature (Curie temperature 𝑇𝐶 ), a paramagnetic material can be externally magnetized to become ferromagnetic again. In the ferromagnetic state, an increased amount of molecular magnetic dipole moments are aligned with the external magnetic field, so the magnetic entropy of the solid refrigerant decreases (there are some exceptions, see [3]). The excess entropy is released in the form of heat

(or thermal entropy, see Eq. (1)). This happens vice versa, such that if the Magneto-Caloric Material (MCM) is demagnetized adiabatically around 𝑇𝐶 , it would cool down. 𝑠𝑡𝑜𝑡𝑎𝑙 (𝐻, 𝑇 ) = 𝑠𝑚𝑎𝑔𝑛𝑒𝑡𝑖𝑐 (𝐻, 𝑇 ) + 𝑠𝑡ℎ𝑒𝑟𝑚𝑎𝑙 (𝑇 )

(1)

Magnetocaloric refrigerators are composed of at least the following parts: a magnet (permanent or electrically powered), a porous regenerator made of MCM, a system that pumps the Heat Transfer Fluid (HTF), and heat exchangers at the thermal sink and source (see Fig. 1). Most published prototypes try to follow a Brayton-like thermodynamic cycle (see Fig. 1): 1. 2. 3. 4.

Fast (quasi-adiabatic) magnetization of the regenerator. Surge of the HTF from the cold source to the hot sink. Fast demagnetization of the regenerator. Surge of the HTF from the hot sink to the cold source.

The reader is referred to [4] for a more exhaustive explanation of the magnetocaloric cycles. The first room temperature magnetocaloric refrigerators used an electromagnet to magnetize the regenerator and they operated at low frequencies. Later on, the electromagnets were substituted by permanent magnets that do not require as much input power to induce the (de)magnetization steps. Two main types of magnetic refrigerators have been built, linear and rotatory (see the exhaustive reviews of Yu et al. [5] and Greco et al. [6] for more information). Linear refrigerators

∗ Corresponding author. E-mail address: [email protected] (I. Mugica).

https://doi.org/10.1016/j.ijthermalsci.2019.106198 Received 6 May 2019; Received in revised form 5 November 2019; Accepted 20 November 2019 Available online 27 November 2019 1290-0729/© 2019 Elsevier Masson SAS. All rights reserved.

International Journal of Thermal Sciences 149 (2020) 106198

I. Mugica et al.

Nomenclature 𝐴 𝐵 𝑐 𝐹 𝑓 𝑁 𝑁𝑢 𝑀 𝑄 𝑅𝑒 𝑆 𝑇 𝑡 𝑡𝑟 𝑉 𝑊

Area/Surface [m2 ] Magnetic field flux [𝑇 ] Heat capacity [J kg−1 K −1 ] Force [N] Frequency [Hz] Demagnetization factor [−] Nusselt number [−] Magnetization [A m−1 ] Heat power [W] Reynolds number [−] Specific entropy [J m−3 K −1 ] Temperature [K] Time [s] Trace of the fluid particles [m] Volume [m3 ] Work power [W]

Fig. 1. (a) Schematic representation of a magnetic refrigerator. (b) Brayton-like thermodynamic cycle of an Active Magnetic Regenerator (AMR).

Acronyms AMR CFC COP GSA HCFC HFC HTF MCE MCM NTU

Active Magnetic Regenerator Chloro Fluoro Carbons Coefficient Of Performance Grid Sensitivity Analysis Hydro Chloro Fluoro Carbons Hydro Fluoro Carbons Heat Transport Fluid Magneto Caloric Effect Magneto Caloric Material Number of Transfer Units

cycle. The inertia of the wheel helps to maintain a constant relative movement of the magnetic field. Examples of this kind of machine are the prototypes published by Zimm et al. [7] and Eriksen et al. [8]. Note that the magnetic field of these last two machines (and of the overwhelming majority, see [5]) is perpendicular to the flow of the HTF. This configuration reduces the gap between the poles of the magnet, but it distributes the magnitude of the MCE in a different spatial dimension (perpendicular to the temperature span of the regenerator and the fluid flow). As checked by Mugica et al. [9], and Mugica et al. [10], simplifying the distribution of the magnetic field along the extra spatial dimension might cause an over prediction of the steady state temperature span of regenerators. The most common regenerator geometries are made of parallel plates or packed particle beds. Parallel plates offer lower pressure drops than packed beds but have a lower cooling power because of the slower heat transfer between the MCM and the HTF. In recent years, several studies have investigated the influence of the regenerator geometry on the performance of the AMR cycle. The following paragraphs sum up the most prominent findings: ∙ Engelbrecht et al. [11] ran passive regenerator cycles (no applied magnetic field) with controlled heat sink and source to measure steady state temperature spans. The regenerators were made of flat, corrugated and dimpled plates. Corrugated plates were made by pressing flat plates onto cylinders. In a similar way, the dimpled plates were fabricated by pressing a flat plate onto spheres. Both dimpled and corrugated plates presented higher temperature spans than flat plates. ∙ Similarly, Šarlah et al. [12] tested various passive regenerator geometries. They compared the pressure drop and heat transfer coefficients of various types of corrugated plates, flat plates and a packed bed. The latter showed higher heat transfer and pressure losses than the rest, but the lowest COP. ∙ Aprea et al. [13] continued investigating packed bed AMR cycles through 1D numerical modeling, more specifically, the aspect ratio of the whole regenerator volume. They concluded that longer and narrower regenerators suffer less from conduction losses, and more from viscous losses. ∙ Nielsen et al. [14] extensively mapped the temperature span and available cooling power of plate regenerators with the help of a 2D numerical model. For all the tested frequencies and utilization ratios (𝑓 = 0.14 − −4 Hz, 𝜑 = 0.14 − 6.4 see Eq. (11)), and for the cases where the mass of MCM was constant, the maximum temperature span and available cooling power decreased with an increasing porosity (𝜀 = 𝑉𝑓 𝑙𝑢𝑖𝑑 ∕𝑉𝑡𝑜𝑡𝑎𝑙 ). The maximum available cooling power always increased with an increasing frequency. This demonstrated that, the available

Greek symbols 𝜀 𝜑 𝜙 𝜆 𝜇𝑓 𝜇𝑟 𝜌 𝜏

porosity [−] Utilization ratio [−] Magnetic scalar potential [A] Thermal conductivity [W m−1 K −1 ] Fluid viscosity [Pa s] Relative magnetic permeability [−] Density [kg m−3 ] Tortuosity [−]

Subscripts 0 app C cond dip f HT int p reg s

Vacuum Applied Curie Conduction Dipolar Fluid Heat Transfer Internal Particle Regenerator Solid

are simpler to build, but are less efficient at higher frequencies of operation due to the considerable inertia of the moving parts. Rotatory configurations work by spinning either the permanent magnet or a wheel of regenerators along the same rotational axis to create an undulatory magnetic field that crosses the regenerator volume every 2

International Journal of Thermal Sciences 149 (2020) 106198

I. Mugica et al.

cooling power favored lower porosities and higher frequencies of operation. A higher porosity was linked to a lower predicted NTU. Higher frequencies of operation increase the velocity of the fluid surges (for a constant utilization ratio) which increases the convective heat transfer between phases (generally the Nusselt number (𝑁𝑢) increases with the Reynolds number (𝑅𝑒), see [15]). However, it leaves less time for the heat transfer to take place. Even though the NTU is lower at higher frequencies, the available cooling power increases by reducing the time of the AMR cycles. ∙ Tušek et al. [16] analyzed the consequences of varying the particle diameter of packed bed regenerators through a 1D numerical model. Decreasing the particle diameter increases the heat transfer surface of the packed be. So at a constant frequency, the maximum cooling load increased with a decreasing particle diameter. More heat could be transferred to the liquid phase at the same surge time. Similar to [14], they found that higher frequencies of operation are linked to a higher optimal mass flow rate. Nonetheless, contrary to [14], the cooling power decreased for higher frequencies of operation, after hitting a maximum. This means that the convective heat transfer rate did not increase enough at higher mass flow rates to compensate for the reduced surge time. This effect was more pronounced at higher temperature spans. Indeed, the temperature span is representative of a cooling load that takes place within the regenerator. The equilibrium between the longitudinal conduction and the cooling power of the machine results in the temperature span at ‘‘zero’’ cooling load. Therefore, at higher temperature spans, the regenerator has less available power to handle a cooling load at the cold heat exchanger. The COP showed a maximum at lower particle diameters and higher frequencies of operation because of the increased viscous losses. The pumping work increased at higher frequencies of operation, due to an increase of the mass flow rate. At the same time, the viscous friction term in the energy equation increased, thus increasing the rejected heat. So, the optimum particle diameter increased at higher operating 𝑓 . ∙ Vuarnoz and Kawanami [17] calculated the performance of active regenerators made of wires that are parallelly stacked to the direction of the flow. The calculations were done by a 1D numerical model and the investigated geometrical parameter was the diameter of the wires. Similar conclusions to past publications were drawn: the reduction of the wire diameter increases the cooling load and the viscous losses. This kind of regenerator showed slightly better COP and cooling power when compared to packed bed regenerators. Comparable to [16], the cooling load hit a maximum when increasing the operating frequency. In this case, it was most certainly due to a Nusselt number that was independent of the fluid velocity. ∙ Tušek et al. [18] calculated the optimal mass flow rate for the maximum COP and cooling power (𝑄𝑐𝑜𝑜𝑙 ) for various geometrical parameters of plate and packed bed regenerators, using a 1D numerical code. The simulations included 2 frequencies of operation (0.5 and 3 Hz). The studied geometric parameters consisted of: the regenerator length, particle diameter, plate thickness, and porosity. Similar to past findings, plate regenerator geometries achieved a higher COP and a lower 𝑄𝑚𝑎𝑥 when compared to packed-bed geometries. In accordance 𝑐𝑜𝑜𝑙 with past authors, they concluded that thinner plates and smaller particles increase 𝑄𝑚𝑎𝑥 to the detriment of increasing the viscous losses. 𝑐𝑜𝑜𝑙 ∙ Tušek et al. [19] tested the performance of different AMR geometries experimentally. The geometries included 3 packed bed regenerators with different particle diameters and shapes, and another 3 plate regenerators with different plate thicknesses, spacings, and different orientations in regard to the applied magnetic field. In terms of the maximum temperature span, the optimum frequency and mass flow rate were related to the overall convective heat transfer properties of the regenerator. A higher heat transfer surface led to higher optimum frequencies and mass flow rates. They concluded that the efficiency of parallel plate regenerators was superior, and that the research towards more efficient regenerators should be focused on ordered structures.

∙ Lei et al. [20] calculated the performance of the following regenerator geometries: packed spheres, packed screens, parallel plates, and micro channels through a 1D numerical model. The mass flow rate was tuned to get a cooling capacity of 100W/kg at a temperature span of 280–300 K. They mapped the COP and the generated entropy for two variables: operational √ frequency and the aspect ratio of the regenerator volume (𝑅𝑎 = 𝐿∕ 𝐴𝑐𝑟𝑜𝑠𝑠 ). The COP did not include the pumping power (𝐶𝑂𝑃 = 𝑄𝑐𝑜𝑙𝑑 ∕(𝑄ℎ𝑜𝑡 − 𝑄𝑐𝑜𝑙𝑑 )), and the volume of the regenerator was always kept constant. Packed screens, parallel plates, and micro channels were found to have a higher optimum COP when compared to packed sphere beds. The following approximation exemplifies the inverse relationship between the hydraulic diameter and the heat transfer area (𝐴𝐻𝑇 ) for the case of packed beds [13]: 𝐴𝐻𝑇 =

6 𝑉 (1 − 𝜀) 𝐷𝑝

(2)

Lei et al. [20] described how the optimum COP (in terms of frequency and 𝑅𝑎 ) for each of the tested geometries at constant porosity, had a maximum for the hydraulic diameter (or equivalently for 𝐴𝐻𝑇 ). As 𝐴𝐻𝑇 increases, the COP increases too due to an increase in the cooling power. As calculated by Tušek et al. [16], past a maximum, the heat generated by friction increases, which increases the heat rejected at the hot heat exchanger and lowers the COP. Geometries with lower 𝐴𝐻𝑇 favored lower frequencies of operation and larger 𝑅𝑎 . As [16] and [19] detailed, the reduced heat transfer (by reducing 𝐴𝐻𝑇 ) demands more surge time and lower conduction losses to achieve the optimum COP. The highest COP values seemed to favor higher porosities (around 0.64) for all the studied geometries, while the porosity of packed spheres is usually around 0.33 (see [21]). ∙ Trevizoli et al. [22] measured the cooling power and COP of packed bed, parallel plate and square pin regenerator geometries. They kept similar Gd mass, 𝐷ℎ , and porosity in the same regenerator housing. Thus, the results were characteristic to the specific geometries. The matrix of pins gave a close performance to the packed-bed regenerator. It had a slightly lower cooling capacity and better efficiency due to the lesser viscous losses. Congruent with previous authors [12,18,19, 22] concluded that packed bed geometries trade heating power for efficiency. A few other authors have investigated the influence of geometry on the perceived magnetic field of a porous structure. In an AMR, the geometry not only determines the heat transfer properties and viscous losses of the HTF, but it also determines how the external magnetic field is funneled into the MCM, which is the fundamental phenomenon that drives the MCE. Typically, the analytic methods to calculate the perceived magnetic field are limited to isolated ellipsoids of constant magnetization under a uniform applied magnetic field [23]. However, averaged solutions have been found for isolated rectangular prisms and cylinders [24], and [25]. Martínez-Huerta et al. [26] experimentally validated an analytical method to approximate the demagnetizing field of a conglomerate of identical 2D particles of relatively high porosity (𝜀 > 0.9). They detailed how the total demagnetizing factor is a sum of the self-demagnetizing factor and the dipolar interaction factor: 𝑁𝑇 ≈ 𝑁𝑠𝑒𝑙𝑓 + 𝑁𝑑𝑖𝑝 ≈ 𝑁𝑠𝑒𝑙𝑓 + (𝑁𝑚𝑎𝑐𝑟𝑜 − 𝑁𝑠𝑒𝑙𝑓 )(1 − 𝜀)

(3)

The demagnetization due to dipolar interaction (𝑁𝑑𝑖𝑝 ) was proposed as a function of the demagnetization of the macroscopic conglomerate 𝑁𝑚𝑎𝑐𝑟𝑜 , of the isolated pins 𝑁𝑠𝑒𝑙𝑓 , and 𝜀. Bjørk and Bahl [27] calculated the total demagnetization of rectangular prismatic packed beds made of spherical particles using a 3D magnetostatic model. In terms of the particle size, the packed beds had monodispersed, normal and lognormal distributions. They also included 3 different aspect ratios for the macroscopic shape of the packed beds: 1, 1.5, and 2. The dimension that changed to accommodate the aspect ratios was parallel to the external uniform field. For low relative magnetic permeabilities, (𝜇𝑟 ≈ 2, characteristic of Gd at 1T) Eq. (3) fit the simulations well. However, 3

International Journal of Thermal Sciences 149 (2020) 106198

I. Mugica et al.

when 𝜇𝑟 was increased, the demagnetization decreased asymptotically. Thus, the total demagnetization 𝑁𝑇 must be a function of 𝜇𝑟 as well. Prior to studies by Bjørk and Bahl [27] and Martínez-Huerta et al. [26], Christensen et al. [28] studied the demagnetization of parallelly stacked plates. They studied different combinations of porosity, plate thickness, total size of the stack, and orientation of the external magnetic field. Their findings followed the logic of Eq. (3). Trevizoli et al. [29] summarized the main challenges that magnetic refrigeration faces. The total MCM mass and regenerator volume should be kept to a minimum, while the regenerator delivers a useful cooling load at a reasonable temperature span. Therefore, the optimal geometry will depend on how three physical phenomena interact: the magnetic field, fluid flow and thermal energy. So far, the literature suggests that larger porosities seem to reduce the demagnetization and friction losses, but increase the magnetized volume and reduce the heat transfer by convection. Those larger porosities are translated into lower optimum frequencies of operation, a higher COP, and lower 𝑄𝑐𝑜𝑜𝑙 . Also, most of these results have been obtained with 1D numerical models. These models employ averaged correlations to model the heat transfer (𝑁𝑢), demagnetization (𝑁𝑇 ), and viscous losses (friction factor). This approach limits the scientific research to geometries where the dimensions of the regenerator are much larger than the characteristic length of the porous media. This might be sufficient when correlations are available and when only the macroscopic behavior of the regenerator is of interest. However, correlations cannot model arbitrary geometries, and do not provide insight on how exactly the physical phenomena develop within the regenerator. The latter is important to understand and engineer the geometry of future regenerators. In order to breach that knowledge gap, a few 2D and 3D models have been published. Petersen et al. [30] and Aprea et al. [31] presented numerical models that can directly solve the convective heat transfer of 2D regenerator geometries. Bouchard et al. [32] was the first to develop a 3D AMR model that would solve the coupled equations of magnetostatic, thermal and fluid fields at every time step. The reader can find a more general review on AMR numerical models by Nielsen et al. [33]. The novelty of the following study is threefold:

[8]). An optimally engineered AMR has an MCM with a 𝑇𝐶 in the vicinity of the thermodynamic cycle that the regenerator is undergoing. However, as Fig. 1 depicts, the thermodynamic cycles are continuously spread throughout the temperature span of the machine. This has led to multi-stage (or multi-layer) designs such as the ones published by Govindappa et al. [34] and Navickaitä et al. [35], where various MCM are chained with an increasing 𝑇𝐶 . That is why only a single stage of MCM is considered in this study. This results in temperature spans that are not as large as the ones reported for a whole regenerator. Monomaterial regenerators are unlikely to become more popular due to their reduced MCE, and the conclusions reached in this study about the geometry of regenerators are extrapolatable to regenerators with multiple stages. 2. Numerical solver Instead drawing upon correlations to model the multiple physical phenomena, the 3D solver described by Mugica et al. [10] directly solves the equations that describe the physics of an AMR cycle: (i) The magnetic field was computed from the magnetostatic scalar potential equation: ∇ ⋅ (𝜇𝑟 ∇𝜙) = 0

(5)

The moving magnet was represented by a Neumann boundary condition 𝜕𝜙 = 𝐻𝑚𝑎𝑔𝑛𝑒𝑡 , and the magnetic permeability was 𝜕𝑛 considered constant within the different domains 𝜇𝑟𝐺𝑑 = 2, 𝜇𝑟𝑎𝑖𝑟 = 1. The accuracy of the last hypothesis was studied by Mugica et al. [10]. (ii) The fluid surges described by the Navier–Stokes equations were precalculated prior to computing the AMR cycle. The fluid properties were chosen to be the same as the water–ethanol mixture (90% − 10%) used by Bahl et al. [36] 𝜌𝑓 = 981 kg∕m3 , 𝑐𝑓 = 4330 J∕(kg K), 𝜆𝑓 = 0.52 W∕(m K), and 𝜇𝑓 = 1.6 mPa s. The influence of temperature dependent properties on the velocity field was analyzed and deemed negligible by Mugica et al. [10]. ( ) In addition, the thermal diffusivity 𝜆𝑓 ∕(𝜌𝑓 𝑐𝑓 ) of this mixture changes less than 2% between 295–300 K (see [37]). (iii) Once the solutions for the magnetic and fluid fields were obtained, the calculation of the fluid and solid thermal fields could be performed by reading the fluid (𝐮) and magnetic (𝐻 = ∇𝜑) fields from a file, and solving Eqs. (6) and (7). All the physical parameters were kept constant, except for the heat capacity of the solid phase 𝑐𝑠 = 𝑐𝑠 (𝐻, 𝑇 ) and the magnetization 𝑀 = 𝑀(𝐻, 𝑇 ), which were extracted from [38] and [39], respectively. The thermal conductivity and density of Gd were the same as those employed by Mugica et al. [9]: 𝜆𝑠 = 10.5 W∕(m K), 𝜌𝑠 = 7900 kg∕m3 .

(i) The presented data has been directly computed employing the thoroughly validated 3D AMR solver presented by Mugica et al. [10]. This solver brings an unprecedented level of accuracy to the mathematical solution of the AMR thermodynamic cycle, which is key to understanding the dynamics of the interacting physical phenomena present in the 3 spatial dimensions. (ii) The influence of a characteristic parameter of porous geometries has been studied: the tortuosity of the fluid surges. This is defined as the fraction of average length that a fluid particle must travel in order to reach the other end of the regenerator. The average is calculated by the mean length of 𝑛 streamlines: 1 ∑𝑛 𝑖=1 ∫ 𝑡𝑟(𝑥)𝑑𝑥 −1 (4) 𝜏= 𝑛 𝐿 (iii) The cases studied here are centered on regenerators made of pins. These pin geometries are characterized by having their largest dimension (easy axis of magnetization) parallel to the applied magnetic field. Another advantage is the reduced conduction losses due to the lack of contact between the pins.

𝜆𝑠 𝜕𝑀(𝐻, 𝑇 ) || 𝜕𝐻 𝜕𝑇 2 + 𝜇0 𝑇 | 𝜕𝑡 = 𝜌 𝑐 (𝐻, 𝑇 ) ∇ 𝑇 𝜕𝑡 𝜕𝑇 |𝐻 𝑠 𝑠

(6)

𝜆𝑓 𝜕𝑇 + 𝐮∇𝑇 = ∇2 𝑇 𝜕𝑡 𝜌𝑓 𝑐 𝑓

(7)

As detailed in [10], the walls of the regenerator are modeled adiabatic, to the exception of the interface between MCM and HTF. At the thermal interface the heat flux is conserved through a 1st order discretization scheme of the temperature gradients. Also, the thermal inflow and outflow boundary conditions are modeled by heat exchangers that assume a perfect mixture of the fluid and hold exactly the volume of the fluid surges. The validation of the 3D AMR solver [10] had two steps. First a validation of the resolution of the separate physical phenomena. Second, a comparison between the AMR cycle simulation and detailed parallel plate experiment by Bahl et al. [36] (see [10]). As it is depicted in Fig. 2, the detailed resolution of the internal magnetic field provides a much better estimation of the regenerator temperature span. By

First, the proposed geometries are analyzed according to the three separate physical phenomena mentioned earlier. Then the full AMR operation of those geometries is directly simulated. This methodology has the objective of elucidating the relationship between the three separate phenomena (magnetic field, fluid flow and thermal energy) and the performance of the AMR. The AMR simulations are centered on a single MCM stage of Gd and have a cyclic magnetic field characteristic of rotatory configurations (much like the prototype by Eriksen et al. 4

International Journal of Thermal Sciences 149 (2020) 106198

I. Mugica et al.

Table 1 Values for the geometrical parameters of Fig. 3. All pin cross-section areas are 1 mm2 . Parameter

Value

r s 𝑟𝑓 j a b A e L

0.5642 mm 1 mm 1.12 mm 106.26◦ 0.7979 mm 1.5958 mm 0.6204 mm 0.25 mm 4 mm

Fig. 2. Validation of the 3D AMR solver of Mugica et al. [10] against the experiments and 2D simulations of Bahl et al. [36]. The temperature span of the plate regenerator (𝛥𝑇𝑟𝑒𝑔 ) was obtained for zero cooling load conditions.

Fig. 3. Description of the simulated pin cross-sections. From left to right and top to bottom: circle, square, foil, ellipse, hexagon, and plate.

Fig. 4. Tested pin regenerator layouts: (a) Staggered, (b) Straight.

solving the magnetic field in 3D, the user obtains a detailed distribution of the MCE and an accurate description of the shape demagnetization and parasitic reluctances of the magnetic circuit (airgaps between the magnet poles and the regenerator). In contrast, the calculus of Bahl et al. [36] employed an instantaneous and homogeneous magnetic field in the solid phase without including any parasitic reluctances of the magnetic circuit. Furthermore, the overwhelming majority of prototypes have their magnetic field perpendicular to the fluid flow and temperature span of the regenerator. Therefore, the detailed simulation of the regenerator becomes a 3D problem. Both Mugica et al. [10] and the present study, show Grid Sensitivity Analyses (GSA) for every case. All the discretized differential terms were represented by second order accurate numerical schemes. The solver employed in this study was developed within the OpenFOAM framework (Open Field Operation And Manipulation).

3.1. Demagnetization In this subsection, the validity of Eq. (3) is tested for the geometries presented in Fig. 3, and the two arrangements of Fig. 4. This is done by simulating the perturbation that the pins exert on a uniform magnetic field in 3D. Since the magnetostatic hypothesis is embraced, the electromagnetic laws of Maxwell can be reduced to a scalar potential for the magnetic field. Hence, 𝐻 = 𝜕𝜙 and the law of continuity for 𝜕𝑛 the magnetic flux reduces to Eq. (5). The different magnetization of the domains (distinct 𝜇𝑟 ) is the source of the difference between the applied and internal magnetic field of the samples. This is usually approximated by a constant value 𝑁 (see [23]): 𝐻𝑖𝑛𝑡 = 𝐻𝑎𝑝𝑝 − 𝑁𝑀

3. Phenomenological study

(8)

It is yet unclear how the arbitrary outlines of the pins affect 𝑁𝑑𝑖𝑝 and 𝑁𝑠𝑒𝑙𝑓 . The height of all the tested geometries was 10 mm, so the √ factor 𝐻∕ 𝐴𝑐𝑟𝑜𝑠𝑠 , and the volume of the pins stays constant through all the cases. The porosity 𝜀 defines the space between pins, which is always equal in both Cartesian directions. The symmetries of the problem were exploited so that the computational domain could be simplified (see Fig. 5). The sides of the domain were set to a zero gradient of the magnetic scalar potential. The applied magnetic field at the top and bottom boundary conditions imposed a uniform magnetic field of 1T parallel to the longitude of the pins

In this section the properties of the pin geometries of Fig. 3 and Table 1 are studied. All cross-sections of Fig. 3 have an area of 1 mm2 . The simulations of the magnetic, fluid, and thermal fields were done by solving the equations presented in the last section. Two types of pin clusters are studied: straight and staggered (see Fig. 4). Straight clusters have the center of the pins aligned to each other, and staggered clusters have alternate columns of pins aligned with the interstice of the neighboring pins. 5

International Journal of Thermal Sciences 149 (2020) 106198

I. Mugica et al.

Fig. 5. Example of the calculated magnetic field on the staggered circular cross-section of 𝜀 = 0.5. The yellow colored boundary cells had a zero gradient condition for the magnetic scalar potential, and the gray colored boundaries represent the interface between pins and the surrounding paramagnetic domain.

Table 2 Standalone demagnetization for the pins described in Fig. 3 and a height of 10 mm. The side far field boundary conditions were chosen so far, there was only a 10−3 T deviation in respect to the applied magnetic field (1T in the Z direction). Type of pin

𝑁𝑠𝑒𝑙𝑓 ⋅ 103

Circle Square Foil Ellipse Hexagon Plate

59.165 59.55 57.72 58.44 59.126 45.337

Fig. 6. Magnified vision of the magnetic field at the top of the staggered plate pins as the cell size decreases (𝜀 = 0.5). The origin of Y axis is 0.125 mm away from the wide side of the plate and its direction is coherent with Fig. 4.

(Z dimension). These boundaries were chosen to be far enough so as not to influence the solution near the pins. The magnetic permeability of the pin was chosen to be 𝜇𝑟 = 2, similar to the one of Gd around 0-1T (see [27] and [10]). Fig. 6 shows the GSA for the case with the finest geometrical features: staggered plates of 𝜀 = 0.5. A minimum cell size of 6 μm at the surface of the pins was necessary to obtain a magnetostatic solution that is independent of the discretization size. The baseline GSA case of a minimum cell size of 24 μm (see Fig. 7) was created by uniformly discretizing the X and Y dimensions of the pin and air domains (the air domain started with 48 μm and a single surface refinement). The lower cell sizes on Fig. 6 correspond to subsequent refinements at the interface by halving the cells. In fact, as the graph in Fig. 6 shows, the discontinuity of the magenetostatic solution is located at the interface between air (𝜇𝑟 = 1) and pin domains. The mesh in the Z dimension had a similar cell size of the X and Y dimensions at the top and bottom faces of the pins (see Fig. 7), and it expanded from there at a rate of 1.25. These elongated cells lowered the total cell count of the mesh, while still giving an accurate solution of the magnetic field. Indeed, the magnetic field vectors are mostly straight in the Z direction, especially away from the interface, where the top and bottom boundary conditions enforce a uniform magnetic field aligned with the Z direction. Table 2 shows the self-demagnetization factor calculated with the magnetostatic solution and Eq. (8) [23]. To calculate 𝑁𝑠𝑒𝑙𝑓 there is no need for a symmetry boundary condition, so the far-field in X and Y directions was chosen to be far enough away so as not to disturb the solution around the pins. There is only a small variation of the demagnetization factor between the different cases of Table 2. Thus, it can be concluded that the outline of the pins does not affect 𝑁𝑠𝑒𝑙𝑓 substantially, as long as

Fig. 7. Details of the mesh employed in the GSA of Fig. 6.

√ the form factor 𝐻∕ 𝐴𝑐𝑟𝑜𝑠𝑠 is kept constant. Next, 𝑁𝑇 is calculated for straight and staggered clusters of pins (Fig. 8). 6

International Journal of Thermal Sciences 149 (2020) 106198

I. Mugica et al.

Fig. 9. Reference for the dimensions specified in Tables 3 and 4.

Table 3 Dimensions of the tested regenerator geometries (see Fig. 9 and Table 4). A similar width (W) was kept to have a comparable magnetic field change between cases. All geometries have 24 pins, and H = 4 mm for all cases. W [mm]

L [mm]

𝜀

0.5

0.625

0.75

0.5

0.625

0.75

Circle and Square Foil Ellipse Hexagon Straight plate Staggered plate

11.314 12.426 12.846 10.668 12 12.67

13.064 14.945 15.385 12.416 16 13.084

16 19.22 19.686 15.349 24 13.849

4.243 3.863 3.737 4.5 4 3.788

4.899 4.282 4.16 5.155 4 4.891

6 4.995 4.877 6.255 4 6.932

Table 4 Number of pins along the X and Y dimensions for the regenerator geometries described in Fig. 9 and Table 3. Fig. 8. Variation of 𝑁𝑇 with the porosity (𝜀 = 0.5, 0.675, 0.75). 𝜀 = 1 is attributed to isolated pins, where 𝑉𝑓 𝑙𝑢𝑖𝑑 ≈ 𝑉𝑡𝑜𝑡𝑎𝑙 . The approximation corresponds to Eq. (3).

Circle and Square Foil Ellipse Hexagon Straight plate Staggered plate

𝑁𝑚𝑎𝑐𝑟𝑜 in the approximation of Eq. (3) was set towards 1, because in the limit case of 𝜀 = 0 the demagnetization of an infinite slab tends to 1. Since the magnetic field is a solenoidal field (∇ ⋅ 𝐵 = 0) when all the applied magnetic field crosses to the MCM: 𝜇𝐼 𝐻𝐼 = 𝜇𝐼𝐼 𝐻𝐼𝐼

𝜇0 𝐻𝑎𝑝𝑝 = 𝜇0 (𝐻𝑖𝑛𝑡 + 𝑀)

X

Y

3 2 2 3 24 3

8 12 12 8 1 8

Fig. 9). The latter constrain ensures similar maximum and minimum bounds of the applied magnetic field. The difference between the minimum and maximum values of the magnetic field drives the MCE, and as [40] investigated, the magnet can be optimized to increase this difference. Note that the regenerators of Table 3 keep a constant volume (𝑉𝑟𝑒𝑔 = 𝑊 ⋅ 𝐿 ⋅ 𝐻) for the cases of equal porosity. This means that, for a constant 𝜀, the magnetized volume and the mass of MCM are constant, which establishes a good baseline for a comparative study. The GSA of the regenerator geometries was done for the case that presents the highest temperature and velocity gradients in the solution: staggered plates of 𝜀 = 0.5. For the GSA, the fluid inlet velocity was set to match a conservative operational point that will not be surpassed: 𝑓 = 1 Hz, and 𝜑 = 10. 𝑈 𝐴 𝜒 (11) 𝜑 = 𝑖𝑛𝑙𝑒𝑡 𝑖𝑛𝑙𝑒𝑡 4𝑓 𝑉𝑠

(9)

(10)

Eq. (10) is an expression analogous to Eq. (8) if 𝑁 = 1. Fig. 8 shows how Eq. (3) overestimates the total demagnetization of straight clusters by almost 0.1 in the worst case (𝜀 = 0.5). Staggered clusters present a slightly higher demagnetization. The tendency seems to be quite linear for both kinds of clusters. Other authors [25,27] have pointed out that the total demagnetization reaches an asymptotic minimum as 𝜇𝑟 increases. So, for higher values of 𝜇𝑟 the discrepancy between the approximation and the computed 𝑁𝑇 should increase. Eq. (8) might be a good option to approximate 𝑁𝑇 for clusters of high porosity (𝜀 > 0.5), and low 𝜇𝑟 difference with the reference case (𝑁𝑠𝑒𝑙𝑓 ).

𝜒=

3.2. Pressure losses and heat transfer In this subsection, the pressure head and heat transfer of the geometries of Fig. 3 are studied for the specific regenerator geometries of Tables 3 and 4. These regenerator geometries were designed to accommodate the presented pins, and to keep a similar width (W, see

𝑓=

𝜌𝑓 𝑐𝑓 𝜌𝑠 𝑐̄𝑠 1 4𝑡𝑓

(12)

(13)

The computed cycle has a magnetic field that drifts along the Y dimension at a constant velocity, and so, the cycle period is divided 7

International Journal of Thermal Sciences 149 (2020) 106198

I. Mugica et al.

Fig. 10. Velocity and temperature profiles of the GSA case at 𝑥 = 1.184 mm away from the cold heat exchanger (see depiction of the regenerator at the bottom).

Fig. 11. (a) Final mesh of the fluid domain in X and Y dimensions with a uniform cell size of 9.75 μm. (b) Velocity and temperature fields of the GSA case (Fig. 10), around the point 𝑥 = 1.184 mm, 𝑦 = 0.

in four equal time periods Eq. (13). Eq. (11) and 𝑐̄𝑠 = 282.64𝐽 ∕(𝑘𝑔 𝐾) yield an inlet velocity of 𝑈𝑖𝑛𝑙𝑒𝑡 = 39.8267 mm∕s for 𝑓 = 1 Hz, and 𝜑 = 10. Top and bottom boundary conditions of the regenerator (limits in Z dimension, see Fig. 9) were modeled as walls, but the sides (limits in Y dimension) were modeled by a symmetrical boundary condition. That way, the influence of the regenerator width on the fluid domain is kept out of the scope of this study. With a conservative approach in mind, the GSA case had the initial temperature of the walls set to 310 K, and the inlet and the fluid initial temperatures to 290 K. Fig. 10 shows the velocity and temperature profiles of a steady state simulation of the fluid domain at the narrow passage located at 𝑥 = 1.184 mm, 𝑦 = 0. This past figure shows how a uniform cell size of 12.5 μm in the X and Y directions is sufficiently small to solve the velocity and temperature fields at the point with the highest local velocity. The final mesh (Fig. 11) was chosen to have a uniform cell size of 9.75 μm with a minimum accuracy of 0.5𝐾∕(310𝐾 − 290𝐾) ≈ 0.025 at 𝑥 = 1.184 mm, 𝑦 = 0. A similar analysis was undergone in the Z dimension to check the solution of the velocity field at the top and bottom walls. Fig. 13 shows the temperature distribution along the X axis (𝑦 = 0 mm) after 0.025s of simulation for different cell sizes in the solid

domain. After the first refinement (50 μm), the solution is not very sensitive to further refinements at the fluid–solid interface (see Fig. 12). However, to minimize the error in the calculation of the heat flux at the interface, a cell size of 12.5 μm at the surface of the pins was finally employed. This size is closer to the discretization size of the fluid domain (9.75 μm) and reduces the difference between the boundary heat flux computed in the solid and fluid mesh. Discretization of the Z dimension of the solid mesh is further explored in Section 4, because the same solid mesh is employed to compute the internal magnetic field. Fig. 14 shows the influence of the time step on the previously analyzed case. An accurate solution could be obtained for 𝛥𝑡 = 1 ms using an implicit second order numerical scheme. 3.2.1. Pressure losses To compare the pressure head of the regenerators of Table 3, 2D steady state flow was directly simulated in the X and Y dimensions for the operational point of 𝑓 = 1 Hz, and 𝜑 = 0.5. By keeping 𝜑 constant, the mass flow is also kept constant. Different regenerator widths provide different inlet velocities (the smaller the width, the higher the inlet velocity, see Table 3). The width (W) of the straight 8

International Journal of Thermal Sciences 149 (2020) 106198

I. Mugica et al.

Fig. 12. Slice at 𝑧 = 0 of the solid mesh refinements at the surface of the pin, near the point located at 𝑥 = 1.184 mm and 𝑦 = 0.

Fig. 15. Pressure head of the regenerator geometries described in Table 3. Lines were drawn for a better visualization, continuous lines refer to staggered cases and discontinuous to straight cases.

Fig. 13. Temperature profiles of the GSA of the solid mesh along the X axis, at 𝑦 = 𝑧 = 0 and 𝑡 = 0.025 s. Vertical dotted lines represent the limits of the solid domain. See Fig. 12 for more information on the cell distribution.

is the 𝛥𝑃 between both cases. This difference is most apparent for the lower 𝜀 cases. The tortuosity of each geometry was calculated for an average of 103 streamlines that were longer than the length of the regenerator (L, see Eq. (4)). Therefore, recirculation zones are out of the scope of 𝜏. Straight clusters had a negligible tortuosity (𝜏 < 0.02). Since the Reynolds number 𝑅𝑒 of these simulations is rather low (𝑅𝑒 ≈ 15), the increase in pressure head of the staggered cases is attributed mostly to the increased shear stress at the walls of the pins. This increase in shear stress comes from an increase of the local velocity of the staggered cases. A higher tortuosity means that the fluid travels a longer distance to end up on the other side of the regenerator. Thus, if the mass flow is constant between straight and staggered cases, the local velocity must increase when the tortuosity increases. Furthermore, in a staggered configuration, the mainstream flow (away from the walls) has a stagnation point at the next line of pins. This exposes the mainstream flow to the shear stress of the wall, thereby increasing the viscous pressure losses. The flow field is always symmetrical around the pins. So if the flow was irrotational, the pressure head of the pins would theoretically be zero. Because of that, the pressure head attributed to the shape of the pins is probably negligible.

Fig. 14. Temperature profiles along the X axis for different time steps, at 𝑦 = 𝑧 = 0 and 𝑡 = 0.025 s. Vertical dotted lines represent the limits of the solid domain.

3.2.2. Heat transfer In order to analyze the heat transfer properties of each of the geometries in Table 3, a series of 2D thermal simulations were done in the X and Y dimensions. They all started with a 2 K linear temperature gradient on both solid and fluid domains. This was approximately the

𝑣𝑒𝑟𝑠𝑢𝑠 staggered cases is the same, so more direct comparisons can be made between these cases. It can be concluded from Figs. 15 and 16 that the higher the tortuosity change between the straight and staggered cases, the higher 9

International Journal of Thermal Sciences 149 (2020) 106198

I. Mugica et al.

Fig. 17. Evolution of the average temperature gradient in the fluid domain at the surface of the pins (𝜑 = 0.5). Lines were drawn for a better visualization, continuous lines refer to staggered cases and discontinuous to straight cases.

Fig. 16. Tortuosity of the staggered regenerator geometries described in Table 3. Lines were drawn for a better visualization.

equivalent maximum temperature gradient seen on [36]. On the first set of simulations, an additional 1.5 K was added to the solid domain to mimic the MCE. This is equivalent to the 𝛥𝑇𝑎𝑑 of Gd for an internal magnetic field of 0.5-1T [38]. Then, the temperature field was simulated for 0.25s without the HTF flow, and from that solution, another 0.25s were run with the fluid flow calculated in the last Section 3.2.1. These time frames are in concordance with the choice of 𝑓 = 1 Hz and Eq. (13). The no-flow/flow alternation brings these experiments a bit closer to the AMR cycle. Figs. 17 and 18 show the spatial average of the normal temperature gradient at the pin boundaries, for the second half of the simulation (when the HTF is flowing). Fig. 17 shows how the regenerators made up of hexagonal, square, and circular pins achieve a slightly higher heat transfer for the staggered cases. This is replicated for all the other porosities in Fig. 18. The decline of the surface temperature gradient of Fig. 18, with an increasing 𝜀, is a consequence of a diminished inlet velocity due to an increase of the regenerator width for the higher 𝜀 cases (mass flow is constant for all cases). This is why, as seen in the literature, the maximum performance of the AMRs with higher 𝜀 is found at higher 𝜑. However, for 𝜑 = 0.5, the staggered plates do not seem to follow this small improvement seen for hexagonal, square and circular pins. An increase in 𝜏 diverts the flow in a perpendicular direction to the initial temperature gradient (it increasingly flows in the Y dimension). Thus, the cold flow does not reach as far down in the X dimension as it does for the straight plates with the same 𝜑. Consequently, for lower 𝜑 cases of high 𝜏 geometries, some stream-down plates do not come

Fig. 18. Overall average temperature gradient in the fluid domain at the surface of the pins for different porosities (𝜑 = 0.5). Lines were drawn for a better visualization, continuous lines refer to staggered cases and discontinuous to straight cases.

in contact with the cold inlet flow, reducing the average 𝑁𝑢 number. Fig. 19 shows the same calculations, but for the case of 𝜑 = 1.5. There, the staggered plates show a much higher heat transfer, when compared to the straight plates working in the same conditions. 10

International Journal of Thermal Sciences 149 (2020) 106198

I. Mugica et al.

4. AMR Simulation In this section the AMR cycle of some of the geometries of Table 3 is simulated. First of all the GSA of the magnetic field solution is checked. The magnetic field is set up to represent the conditions of a rotatory regenerator. The poles of the magnet represented by a rectangular = 1𝑇 (parallel to Z) swipe cyclically at boundary condition of 𝜇0 𝜕𝜙 𝜕𝑛 a constant speed, always moving negative to Y (see Figs. 21 and 22). The speed of the magnet is determined by the cycle frequency (1 Hz). As implied by Eq. (13), (de)magnetization time is equal to the fluid surge time (𝑡𝑚𝑎𝑔 = 𝑡𝑓 ). In order to model the influence of neighboring regenerators on the internal magnetic field, half a regenerator is placed on both sides of the investigated geometry (green plates on Fig. 21). The neighboring half regenerators are linked to each other through a cyclic Arbitrary Mesh Interface (cyclicAMI). This boundary condition transmits the information between boundary pairs (black contours in Fig. 21 at the end of the green regenerators). Effects of curvature were neglected, so the circumference of the machine located in the Z dimension is hypothesized to be substantially larger than the height of the pins (4 mm). The air gap was 1 mm. The modeled regenerators were designed to represent a single MCM stage in a chain of multimaterial regenerators, with each stage working near their 𝑇𝐶 . So the boundary conditions of the cold and hot sides mirror the magnetic field solution = 0), as if another stage was present next to the cold (𝑥 = 0) and ( 𝜕𝜙 𝜕𝑛 hot (𝑥 = 𝐿) boundaries. Fig. 21 shows the magnetic field GSA on the regenerator with the finest geometry details (staggered plates of 𝜀 = 0.5). The center of the magnet poles in the Y dimension was set at the interface of the gray regenerator (the one used in the calculation of the AMR cycle) and the neighboring green regenerator. This configuration was set up to test the proper transmission of information though the interface between the gray and green regenerators, as well as through the cyclic boundary condition. The size of 12.5 μm in the X and Y dimensions employed previously, for the MCM domain (see Figs. 12 and 13) seem to be fine enough for calculating the magnetic field. The cells in the Z direction were approximately half the size of X and Y dimensions at the top of the pins, where the solution of the magnetic field changes the most. From there, the cells expanded approximately at a rate of 1.25 in the Z dimension. To ensure an accurate representation of the boundary conditions, the boundary surfaces had an extra level of refinement relative to the surface of the pins (see Fig. 22). Sections 3.2.1 and 3.2.2 show that 𝜏 needs to be larger than 2 to have a glaring effect on the regenerator properties. The computational cost of simulating an AMR cycle is relatively high (see [10]), so only the following regenerator geometries were simulated: staggered plates (𝜀 = 0.5, 0.625) and circles (𝜀 = 0.5, 0.625, 0.75), straight plates (𝜀 = 0.5, 0.625, 0.75). The last two were chosen so that the reader can make comparisons with other published regenerators. Fig. 24 shows the temperature span of the regenerators at zero cooling load, with a fixed temperature at the hot heat exchanger 300 K. A steady state solution was reached when 𝛥𝑇𝑟𝑒𝑔 changed less than 0.01 K between cycles. The best results are observed for staggered plates. It seems that the increased heat transfer and the reduced conduction losses that come from increasing 𝜏 enable the staggered circle and plate regenerators to work at noticeably higher temperature spans. Staggered circles and straight plates show a decreasing 𝛥𝑇𝑟𝑒𝑔 , when 𝜀 increases. As Figs. 8, 18 and 20 demonstrate, by increasing 𝜀, demagnetization decreases, 𝑁𝑢 decreases (at constant 𝜑), and the conduction losses are slightly reduced. Thus in the case of staggered circles and straight plates, the reduction of 𝑁𝑢 overwhelms the two other positive effects of increasing 𝜀. Also, the maximum 𝛥𝑇𝑟𝑒𝑔 comes at higher 𝜑, for higher 𝜀 cases. This maintains coherence with past findings, as more mass flow (𝜑) is needed for the higher 𝜀 cases to reach the same 𝑁𝑢 value. The staggered plates show an inverse behavior: 𝛥𝑇𝑟𝑒𝑔 increases when 𝜀 increases. This means that the benefit of the decreased demagnetization and conduction losses outweighs the reduced 𝑁𝑢 of higher

Fig. 19. Evolution of the average temperature gradient in the fluid domain at the surface of the plate pins (𝜑 = 1.5). Lines were drawn for a better visualization, the continuous line refers to the staggered case and the discontinuous to the straight case.

Fig. 20. Generated entropy by conduction in the X dimension after 0.25s, for an initial 𝛥𝑇 = 2 K and with no fluid flow. Lines were drawn for a better visualization, continuous lines refer to staggered cases and discontinuous to straight cases.

From Figs. 17 and 19 it can be concluded that, if 𝜑 is high enough, the heat transfer between phases increases when 𝜏 increases. So, 𝜏 acts as a multiplier of the inlet velocity. It increases the local velocity within the regenerator, and thus the convective heat transfer. A second set of simulations was performed to evaluate the longitudinal heat diffusion. For that, no additional temperature was added to the solid phase, and the calculations were run without the fluid flow for 0.25s. Fig. 20 shows the volume average of the total generated entropy by conduction along the X dimension (see Eq. (14)). ( ) ( ) 𝜆𝑓 𝜕𝑇𝑓 2 𝜆𝑠 𝜕𝑇𝑠 2 𝑆𝑐𝑜𝑛𝑑 = + 𝑑𝑡 (14) 2 2 ∫ 𝑇 𝜕𝑥 𝜕𝑥 𝑇𝑠 𝑓

The difference between staggered and straight cases is minimal, with the exception of the plate pins (see Fig. 20). In the case of the straight plates, the solid phase holds the entire initial temperature gradient (2 K/4 mm), which is translated in an increase of the conduction losses due to 𝜆𝑠 ≫ 𝜆𝑓 . For the staggered plates, the initial gradient is split between the fluid and solid domains many times. So the thickness of the staggered plates holds a much smaller temperature gradient. The diminishing trend of 𝑆𝑐𝑜𝑛𝑑 with 𝜀 in Fig. 20 is caused in part, by a lower initial temperature gradient due to the larger length (L) of the higher 𝜀 cases, and in part by the increased proportion of fluid with lower 𝜆. The latter lowers the volumetric average of 𝑆𝑐𝑜𝑛𝑑 . 11

International Journal of Thermal Sciences 149 (2020) 106198

I. Mugica et al.

Fig. 23. Magnetic and Velocity fields of the AMR cycle simulations, at the middle of the hot surge, for the case where 𝜑 = 0.6 and 𝜀 = 0.625. Top to bottom: staggered plate, staggered circle, straight plate. Velocity magnitude contours show how the higher tortuosity cases have a higher local velocity.

Fig. 21. (a) Magnetic field at the top of the pins of the GSA case at 𝑥 = 1.184 mm (vertical dotted lines represent the limits of the gray regenerator employed in the calculation of the AMR cycle). (b) Origin of the plotted Y axis, and a slice of the computational mesh at 𝑥 = 1.184 mm.

Fig. 22. 3D view of the case plotted in Fig. 21 with the rectangular boundary condition representing the poles of the magnet ( 𝜕𝜙 = 1𝑇 ). 𝜕𝑛

porosities. The benefits of reducing the conduction losses can be clearly seen for 𝜑 = 0.5. Figs. 8 and 17 show that 𝑁𝑇 and 𝑁𝑢 for staggered and straight plates of 𝜀 = 0.5 and 𝜑 = 0.5 are similar. That provides a good reference to understand the reduced conduction losses of the staggered 𝑣𝑒𝑟𝑠𝑢𝑠 straight cases in 𝛥𝑇𝑟𝑒𝑔 (see Fig. 20). The way in which the MCM has to be arranged to increase 𝜏, reduces the conduction losses. So, 𝛥𝑇𝑟𝑒𝑔 increases in Fig. 24 between the staggered and straight plate cases of 𝜀 = 0.5 and 𝜑 = 0.5. A steep increase of 𝛥𝑇𝑟𝑒𝑔 is seen for staggered plates of 𝜀 = 0.625 between 𝜑 = 0.6 − 1. As explained in Section 3.2.2 the increased 𝜏 diverts the flow in a perpendicular direction to the temperature gradient of the regenerator. This reduces the transport of temperature in the X dimension for lower 𝜑 cases. This effect becomes more prominent at higher 𝜀 as demonstrated by Fig. 20. Once 𝜑 reaches a higher value, the regenerator is able to pump heat in the direction of 𝛥𝑇𝑟𝑒𝑔 . A drawback of the staggered plates is the increased pressure head. Fig. 25 shows the pressure head of the regenerators calculated in

Fig. 24. Temperature span of the staggered plates (red), staggered circles (green), and straight plates (blue) for a zero cooling load and a constant 300 K in the hot heat exchanger. Symbols represent different values of 𝜀.

Fig. 24. As commented earlier, tortuosity increases the pressure head of the regenerators. However, the minimum ratio between the rejected 𝑄

heat and the pumping work was 𝑊 ℎ𝑜𝑡 = 368.4 (staggered plates of 𝑝𝑢𝑚𝑝 𝜀 = 0.5 at 𝜑 = 1.8). Thus, the magnetic work is still significantly higher than the pumping work for these geometries. The magnetic work, as expressed by Kitanovski and Egolf [4], can be calculated the following 12

International Journal of Thermal Sciences 149 (2020) 106198

I. Mugica et al.

5. Conclusion A detailed phenomenological analysis of AMR pin regenerator geometries was undergone. Demagnetization, pressure head and heat transfer properties were analyzed separately before computing the temperature span of their AMR cycle. All regenerators had the same amount of MCM distributed equally between the pins. The outer shape of the regenerators was kept as similar as possible, in order to maintain a similar magnetic field change and flow velocity. The introduction of a larger tortuosity (by changing the arrangement and shape of the pins) resulted in the following: slight increase of the total demagnetization (less than 0.1), increase of the pressure head, increase of the convective heat transfer, and decrease of the conduction losses. The increase of the pressure head and convective heat transfer is caused by an increase of the local velocity within the regenerator, due to the longer trajectory the HTF must travel at a constant mass flow. Simulations of the AMR cycles indicate that the main contribution to the increased temperature span of the higher 𝜏 cases is due to the reduction of the conduction losses. Furthermore, 𝛥𝑇𝑟𝑒𝑔 seems to increase with 𝜀 for the highest tortuosity cases. This means that the decrease of conduction and demagnetization losses of a higher 𝜏 outweighs the reduced inlet velocity for higher 𝜀 cases. The increase of tortuosity could be an alternative solution to increase the power of magnetocaloric regenerators. Thus far, the strategy to increase the cooling power of packed beds has been to increase the heat transfer surface by employing smaller particles. However, particles in packed beds are in contact with each other, which increases the conduction losses as 𝜀 decreases. Staggered arrangements of elongated pin cross sections, have inherently less conduction losses, and could in principle provide a similar increase in convective heat transfer due to the local increase of the fluid velocity. Furthermore, if most of the heat transfer is controlled by the fluid surges (due to a very different 𝑁𝑢 number between (de)mangetization and surge processes), the thermodynamic cycle will be closer to a Brayton-like cycle. Also, to achieve a higher tortuosity, it is necessary to develop channels for the HTF that are perpendicular to the temperature span of the regenerator. In many cases, this perpendicular direction is the same as the swiping motion of the magnet (Y dimension, see the HTF channels of staggered plates in Fig. 23). Since there is no practical reason for magnetizing the HTF, and higher 𝜏 geometries seem to perform better at higher 𝜀, it might be possible to reduce the material of the permanent magnet that sandwiches the HTF channels that are parallel to Y.

Fig. 25. Pressure head of the regenerators simulated in Fig. 24.

way: 𝑊𝑚𝑎𝑔 =

∫

𝜇0 𝑀(𝐻, 𝑇 )𝑑𝐻

(15)

Hence, knowing 𝑀(𝐻, 𝑇 ), the temperature field, and the internal magnetic field should be sufficient to calculate 𝑊𝑚𝑎𝑔 . However, the calculation of 𝑊𝑚𝑎𝑔 did not match the values of 𝑄ℎ𝑜𝑡 as it should when applying the first law at zero cooling load. Therefore, the authors did not consider these results to be trustworthy enough for publication. There are a couple of reasons for the observed discrepancy:

Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

(i) The internal magnetic field was precalculated as an independent variable from the temperature field (Eq. (5) with 𝜇𝑟𝐺𝑑 = 2). A more precise internal magnetic field could have been obtained by solving the following nonlinear magnetostatic equation at runtime:

Acknowledgments This work is part of the NSERC chair on industrial energy efficiency established at Université de Sherbrooke in 2014 with the support of Hydro-Québec, Canada, Ressources Naturelles Canada (CanmetENERGIE) and Rio Tinto Alcan, Canada, who are gratefully acknowledged for their financial support. The calculations have been performed on the Mammouth Parallèle 2 cluster of Compute Canada, which is also gratefully acknowledged.

∇ ⋅ (𝜇𝑟 (𝐻, 𝑇 )∇𝜙) = 0 (16) 𝑀(𝐻, 𝑇 ) 𝜇𝑟 (𝐻, 𝑇 ) = 1 + (17) 𝐻 The hypothesis of a constant 𝜇𝑟 to calculate the temperature field was thoroughly checked by Mugica et al. [10]. Solving Eq. (16) at runtime is certainly possible and [10] presented a solver that would operate in this way. However, that algorithm runs at a substantially larger computational cost, that would have made the present study intractable. (ii) The literature is sparse on the physical properties of MCMs. Therefore, 𝑐𝑠 (𝐻, 𝑇 ) and 𝑀(𝐻, 𝑇 ) of Gd had to be obtained from different sources. This probably introduces some kind of inconsistency between the solved temperature field and 𝑀(𝐻, 𝑇 ).

References [1] W. Goetzler, Energy Savings Potential and R&D Opportunities for Non-VaporCompression HVAC Technologies, Tech. Rep., U.S. Department of Energy, 2014. [2] J. Smart, Effective Field Theories of Magnetism, Saunders, Philadelphia, 1966. [3] T. Krenke, E. Duman, M. Acet, E. Wassermann, X. Moya, L. Mañosa, A. Planes, Inverse magnetocaloric effect in ferromagnetic Ni–Mn–Sn alloys, Nature Mater. 4 (2005) 450–454. 13

International Journal of Thermal Sciences 149 (2020) 106198

I. Mugica et al.

[23] A.H. Morrish, The Physical Principles of Magnetism, John Wiley & Sons, Inc., Massachusetts, USA, 2001. [24] A. Aharoni, Demagnetizing factors for rectangular ferromagnetic prisms, J. Appl. Phys. 83 (6) (1998) 3432–3434. [25] R. Prozorov, V. Kogan, Effective demagnetizing factors of diamagnetic samples of various shapes, Phys. Rev. A 10 (2018) 014030. [26] J.M. Martínez-Huerta, J.D.L.T. Medina, L. Piraux, A. Encinas, Configuration dependent demagnetizing field in assemblies of interacting magnetic particles, J. Phys.: Condens. Matter 25 (22) (2013) 226003. [27] R. Bjørk, C.R.H. Bahl, Demagnetization factor for a powder of randomly packed spherical particles, Appl. Phys. Lett. 103 (10) (2013) 102403. [28] D.V. Christensen, K.K. Nielsen, C.R.H. Bahl, A. Smith, Demagnetizing effects in stacked rectangular prisms, J. Phys. D: Appl. Phys. 44 (21) (2011) 215004. [29] P.V. Trevizoli, T.V. Christiaanse, P. Govindappa, I. Niknia, R. Teyber, J. Barbosa, A. Rowe, Magnetic heat pumps: An overview of design principles and challenges, Sci. Technol. Built Environ. 22 (5) (2016) 507–519. [30] T.F. Petersen, K. Engelbrecht, C.R. Bahl, B. Elmegaard, N. Pryds, A. Smith, Comparison between a 1d and a 2d numerical model of an active magnetic regenerative refrigerator, J. Phys. D: Appl. Phys. 41 (10) (2008) 105002. [31] C. Aprea, A. Greco, A. Maiorino, C. Masselli, A comparison between rare earth and transition metals working as magnetic materials in an AMR refrigerator in the room temperature range, Appl. Therm. Eng. 91 (2015) 767–777. [32] J. Bouchard, H. Nesreddine, N. Galanis, Model of a porous regenerator used for magnetic refrigeration at room temperature, Int. J. Heat Mass Transfer 52 (5–6) (2009) 1223–1229. [33] K. Nielsen, J. Tusek, K. Engelbrecht, S. Schopfer, A. Kitanovski, C. Bahl, A. Smith, N. Pryds, A. Poredos, Review on numerical modeling of active magnetic regenerators for room temperature applications, Int. J. Refrig. 34 (3) (2011) 603–616. [34] P. Govindappa, P.V. Trevizoli, O. Campbell, I. Niknia, T.V. Christiaanse, R. Teyber, S. Misra, M.A. Schwind, D. van Asten, L. Zhang, A. Rowe, Experimental investigation of MnFeP1-x As x multilayer active magnetic regenerators, J. Phys. D: Appl. Phys. 50 (31) (2017) 315001. [35] K. Navickaitä, H.N. Bez, T. Lei, A. Barcza, H. Vieyra, C.R. Bahl, K. Engelbrecht, Experimental and numerical comparison of multi-layered la(fe,si,mn)13hy active magnetic regenerators, Int. J. Refrig. 86 (2018) 322–330. [36] C.R.H. Bahl, T.F. Petersen, N. Pryds, A. Smith, A versatile magnetic refrigeration test device, Rev. Sci. Instrum. 79 (9) (2008) 093906. [37] I.H. Bell, J. Wronski, S. Quoilin, V. Lemort, Pure and pseudo-pure fluid thermophysical property evaluation and the open-source thermophysical property library coolprop, Ind. Eng. Chem. Res. 53 (6) (2014) 2498–2508. [38] M. Risser, C. Vasile, B. Keith, T. Engel, C. Muller, Construction of consistent magnetocaloric materials data for modelling magnetic refrigerators, Int. J. Refrig. 35 (2) (2012) 459–467. [39] J.S. Lee, Evaluation of the magnetocaloric effect from magnetization and heat capacity data, Phys. Status Solidi b 241 (7) (2004) 1765–1768. [40] R. Bjørk, C. Bahl, A. Smith, D. Christensen, N. Pryds, An optimized magnet for magnetic refrigeration, J. Magn. Magn. Mater. 322 (21) (2010) 3324–3328.

[4] A. Kitanovski, P.W. Egolf, Thermodynamics of magnetic refrigeration, Int. J. Refrig. 29 (1) (2006) 3–21. [5] B. Yu, M. Liu, P.W. Egolf, A. Kitanovski, A review of magnetic refrigerator and heat pump prototypes built before the year 2010, Int. J. Refrig. 33 (6) (2010) 1029–1060. [6] A. Greco, C. Aprea, A. Maiorino, C. Masselli, A review of the state of the art of solid-state caloric cooling processes at room-temperature before 2019, Int. J. Refrig. 106 (2019) 66–88. [7] C. Zimm, A. Boeder, J. Chell, A. Sternberg, A. Fujita, S. Fujieda, K. Fukamichi, Design and performance of a permanent-magnet rotary refrigerator, Int. J. Refrig. 29 (8) (2006) 1302–1306. [8] D. Eriksen, K. Engelbrecht, C. Bahl, R. Bjørk, K. Nielsen, A. Insinga, N. Pryds, Design and experimental tests of a rotary active magnetic regenerator prototype, Int. J. Refrig. 58 (2015) 14–21. [9] I. Mugica, S. Poncet, J. Bouchard, Entropy generation in a parallel-plate active magnetic regenerator with insulator layers, J. Appl. Phys. 121 (7) (2017) 074901. [10] I. Mugica, S. Poncet, J. Bouchard, An open source DNS solver for the simulation of active magnetocaloric regenerative cycles, Appl. Therm. Eng. 141 (2018) 600–616. [11] K. Engelbrecht, K. Nielsen, N. Pryds, An experimental study of passive regenerator geometries, Int. J. Refrig. 34 (8) (2011) 1817–1822. [12] A. Šarlah, J. Tušek, A. Poredoš, Comparison of thermo-hydraulic properties of heat regenerators applicable to active magnetic refrigerators, Strojniški Vestnik J. Mech. Eng. 58 (1) (2012) 16–22. [13] C. Aprea, A. Greco, A. Maiorino, Modelling an active magnetic refrigeration system: A comparison with different models of incompressible flow through a packed bed, Appl. Therm. Eng. 36 (2012) 296–306. [14] K. Nielsen, C. Bahl, A. Smith, N. Pryds, J. Hattel, A comprehensive parameter study of an active magnetic regenerator using a 2d numerical model, Int. J. Refrig. 33 (4) (2010) 753–764. [15] N. Wakao, S. Kaguei, Heat and Mass Transfer in Packed Beds, Gordon and Breach Science Publishers, 1983. [16] J. Tušek, A. Kitanovski, I. Prebil, A. Poredoš, Dynamic operation of an active magnetic regenerator (AMR): numerical optimization of a packed-bed AMR, Int. J. Refrig. 34 (6) (2011) 1507–1517. [17] D. Vuarnoz, T. Kawanami, Numerical analysis of a reciprocating active magnetic regenerator made of gadolinium wires, Appl. Therm. Eng. 37 (2012) 388–395. [18] J. Tušek, A. Kitanovski, A. Poredoš, Geometrical optimization of packed-bed and parallel-plate active magnetic regenerators, Int. J. Refrig. 36 (5) (2013) 1456–1464. [19] J. Tušek, A. Kitanovski, S. Zupan, I. Prebil, A. Poredoš, A comprehensive experimental analysis of gadolinium active magnetic regenerators, Appl. Therm. Eng. 53 (1) (2013) 57–66. [20] T. Lei, K. Engelbrecht, K.K. Nielsen, C.T. Veje, Study of geometries of active magnetic regenerators for room temperature magnetocaloric refrigeration, Appl. Therm. Eng. 111 (2017) 1232–1243. [21] A. Kitanovski, J. Tušek, U. Tomc, U. Plaznik, M. Ozbolt, A. Poredoš, Magnetocaloric Energy Conversion, Springer, Switzerland, 2015. [22] P.V. Trevizoli, A.T. Nakashima, G.F. Peixer, J.R. Barbosa, Performance assessment of different porous matrix geometries for active magnetic regenerators, Appl. Energy 187 (2017) 847–861.

14