A numerical analysis of a magnetocaloric refrigerator with a 16-layer regenerator

A numerical analysis was conducted to study a room temperature magnetocaloric refrigerator with a 16-layer parallel plates active magnetic regenerator (AMR). Sixteen layers of LaFeMnSiH having different Curie temperatures were employed as magnetocaloric material (MCM) in the regenerator. Measured properties data was used. A transient one dimensional (1D) model was employed, in which a unique numerical method was developed to significantly accelerate the simulation speed of the multi-layer AMR system. As a result, the computation speed of a multi-layer AMR case was very close to the single-layer configuration. The performance of the 16-layer AMR system in different frequencies and utilizations has been investigated using this model. To optimize the layer length distribution of the 16-layer MCMs in the regenerator, a set of 137 simulations with different MCM distributions based on the Design of Experiments (DoE) method was conducted and the results were analyzed. The results show that the 16-layer AMR system can operate up to 84% of Carnot cycle COP at a temperature span of 41 K, which cannot be obtained using an AMR with fewer layers. The DoE results indicate that for a 16-layer AMR system, the uniform distribution is very close to the optimized design.


Numerical Method and Verification
The numerical simulations of the 16-layer magnetocaloric parallel plates regenerator are based on layers LaFeMnSiH with different Curie temperatures and a water/glycol mixture. The magnetic refrigerator is comprised of parallel plates with 200 mm as length L, 0.75 mm as height H, and a unit length width as shown in Fig. 1(a).
x and y are the direction of fluid flow and the direction perpendicular to the fluid flow, respectively. Figure 1(b) shows the entropy change with temperature for each of the 16 MCMs used in the model when the applied magnetic field is varied from 0 to 1.6 T, which also reveals the different Curie temperatures for the MCM.
The effect of plate orientation with respect to external magnetic field on the magnetization is not considered for simplicity. The porosity ε in present work is 1/3. The cold end temperature T C and hot end temperature T H are fixed at 266 K (−7.15 °C) and 307 K (33.85 °C), based on the Curie temperature of MCM materials at the two ends, so the temperature span in present work is fixed at 41 K. The maximum applied magnetic field B max is 1.5 T. The density and thermal conductivity of MCM are 7000 kg/m 3 and 10 W/(m⋅K), respectively, as reported by the MCM provider. Therefore, the total mass of MCM in one plate is 1.05 kg. The properties of the water/glycol mixture are obtained from 36 as a function of temperature except the density, since the variance of fluid density has a trivial effect on the results. Similar to the single-layer AMR refrigeration system, the refrigeration cycle with 16-layer parallel plates regenerator also consists of four steps, which are magnetization, cold-to-hot blow, demagnetization, and hot-to-cold blow. The details of the cycle can be found in the literature 37 and will not described here. Figure 1(c) shows the variation of the fluid mass flow and magnetic field during the refrigeration cycle.
In order to analyze and optimize the magnetic refrigerator system with a multi-layer regenerator, a 1D numerical model was employed, which applies the first law of Thermodynamics to the MCM plates and the heat transfer fluid. The governing equations that describe the heat transfer in MCM plates and heat transfer fluid are shown in Equations (1) and (2), respectively, which are coupled by a convective heat transfer term. This mathematical model has been widely used and validated by different investigators 30,[37][38][39] .
In these two equations, T, k e , ρ, and c denote temperature, effective thermal conductivity, density, and specific heat, respectively, while subscripts s and f indicate MCM and fluid respectively. a, ε, s, B, p and u are volume-specific surface area, porosity, MCM entropy, applied magnetic field, fluid pressure, and fluid velocity, respectively. The k e,s and k e,f are defined as Engelbrecht suggested 37 , based on the axial dispersion for a fluid flowing between two infinite plates from Beard 40 , Scientific RepORts | 7: 13962 | DOI:10.1038/s41598-017-14406-9 where k s and k f are the thermal conductivities of MCM and fluid, respectively. Pr and Re are the Prandtl Number and Reynolds Number, respectively. k e,s relates the actual rate of axial conduction through the composite MCM/ fluid matrix in the absence of fluid flow to the rate of conduction heat transfer that would be experienced by a comparable solid piece of material. k e,f considers an augmentation of the thermal conductivity in the fluid due to the eddy mixing of the fluid that occurs when fluid flows 37 .
In Equations (1) and (2), the MCM entropy and specific heat vary with temperature as well as applied magnetic field. However, these two parameters are not independent, i.e. the specific heat of MCM can be calculated from entropy variation with temperature, c s = T s (∂s/∂T s ). Therefore, for the 16 MCMs with different Curies temperatures, 16 different sets of s property data have to be utilized, which are from the experimental measurements provided by the material provider.
The convection heat transfer coefficient, h sf , is calculated from correlations for the flow between parallel plates 41 ,

Nu Gz
where Gz is the Graetz number in the regenerator, which is defined as The friction factor, which includes the viscous developing region of the flow, is shown in Equation (10).
The boundary conditions are shown in Table 1. The cooling capacity, rejected heat, COP, and utilization, are calculated using the following equations, c c r In the above equations, τ is the cycle time, which is the multiplicative inverse of frequency f. τ f is the time for cold flow and hot flow, since they are identical in present work.
 m f is the mass flow rate.  W pump is the pumping power and calculated using the pressure drop and the flow rate. Utilization is the ratio between the heat capacity of the fluid and that of the magnetocaloric material during a single fluid flow period 43 . The cooling power density (Q c ) is defined as the generated cooling power by unit mass of MCM.
The 1D model was implemented in MATLAB ® . The MATLAB ® code is based on 37 , which starts from an initial temperature distribution and terminates after a cyclical steady state has been achieved. Steady state is defined as when the relative change in the energy (fluid energy + regenerator energy) of the regenerator from cycle to cycle is less than a specified tolerance. To verify the developed model, simulation of a single-layer AMR made of Gd using water as the heat transfer fluid with published results from Petersen et al. 44 (numerical results) were compared, since Petersen et al. 44 provide the comprehensive parameters in their paper to benefit the comparison. The comparisons in Fig. 2(a) show a good agreement between each other. Moreover, a comparison was also conducted against experimental data from Engelbrecht et al. 45 . As shown in Fig. 2(b), although, there is an over-prediction of Q c from our numerical model due to the lack of external losses.The single-layer code needs to be extended for multi-layer AMRs. Extension of the original single-layer model to multi-layer suffered from excessively long computation time. The speed issue may be insignificant if only four or six layers are considered. However, for a 16-layer AMR, the increase in computational time makes the model undesirable, e.g. the time required for a 16-layer case using this model running 10 cycles was 10809 seconds, while a one-layer case under the same conditions needed only 1804 seconds. In the simulation, the MCM properties (entropy and specific heat) are preprocessed as matrices in MATLAB that are used as lookup tables, so for multi-layer AMR, every MCM layer has two property matrices. The property matrices need to be interpolated to obtain the entropy and specific heat for a specific temperature and magnetic field. It was found that the most time consuming process in the simulation is the interpolation in properties matrices. The time consumption increases significantly with the number of layers as 2*N layer property matrices have to be interpolated every time step in the simulation. It is also found that the calculation speed of the native MATLAB interpolation function is almost independent of the size of matrix, i.e. the time consumption to conduct an interpolation in a large matrix is almost as same as in a small one. This finding was exploited to reduce the simulation computation time. Instead of using multiple matrices for every layer separately, the multi-layer properties were integrated into two large matrices for entropy and specific heat, individually. A mathematical transformation was employed to change the temperature index of different layers, i.e. the temperature index for layer N will add N*100 K, so the properties can be identified for each layer in one matrix. As a result, the speed of 16-layer cases has been significantly increases by a factor of 6. As a result, the 16-layer case solved by the new speed improvement model only needs 1867 seconds to run 10 cycles, which is very close to the time needed for single-layer model. Figure 2(c) to (f) shows the results of the mesh and time independent study, indicating that COP is sensitive to the time step, so a sufficiently small time step is required in order to accurately capture the dynamics of the 16-layer AMR. Additionally, Fig. 2(e) and (f) show that mesh size does not greatly impact the COP and Q c from N x = 256 to 512. According to the mesh and time step study, N t = 40000 (Δt = 1/4000 s) and N x = 256 (Δx = 0.78 mm) are chosen as the parameters for an optimal balance between accuracy and simulation speed.

Results and Discussions
Performance Analysis of the 16-layer AMR. Q c and COP for different utilization and frequencies for single layer AMR have been well investigated by the previous work experimentally, e.g. Trevizoli et al. and Trevizoli et al. 46,47 , so the discussions in this section will focus on the performance of the 16-layer AMR.
With the configuration shown in Fig. 1, a series of simulations were conducted to investigate the effect of utilization U on the performance of the refrigerator. Figure 3(a) depicts the cooling power Q c change with utilization at different operating frequencies. For a fixed frequency, the cooling power increases with U initially and then decreases after reaching a maximum value. When the utilization is low, the heat transfer fluid moves slowly and small quantity of low temperature fluid reaches the cold side. Therefore, only little cooling power is produced. On the other hand, when the fluid moves too quickly, the temperature profile of the AMR is disturbed, resulting in a loss of cooling power as well 48 . Note that the cooling power is negative, and not shown in this plot, when U > 0.255 and f = 0.5 Hz, indicating that the disturbed temperature profile leads to T f (x = L) > T c in the cooling processing and the refrigerator does not produce cooling effect under this condition. Consequently, for a 16-layer AMR running under a certain frequency, there exists a U, as shown in Fig. 3(a), that produces maximum cooling power. The corresponding U of the maximum cooling power increases as frequency decreases. It is because a higher U (e.g. fluid speed) is required to disturb the temperature profile in a lower frequency AMR. Figure 3(a) also indicates that the maximum cooling power increases with frequency. This occurs because increasing the operation frequency not only increases the rate of fluid motion, but also enhances the rate of entropy change (see Equation (1)) of the MCM to generate more cooling power. Figure 3(b) illustrates the COP change with U at different operating frequencies.
The COP at f = 0.5 and U > 0.255 is not shown here because it is less than zero due to the negative cooling power discussed above. Results indicate that when f = 0.5, the COP always decreases with U from 0.05 to 0.25. The trends of the COP of other frequencies are similar to the cooling power, i.e. there exists a utilization that results in a maximum COP. Figure 3(b) indicates that the maximum COP at f = 0.167, 0.1, and 0.05 Hz are 4.7 and 5.5 and 5.0 occurring at U = 0.1, 0.15, and 0.2 respectively, so the COP of the 16-layer AMR system reaches its maximum at frequency around 0.1 Hz. This is attributed to the competition between Q c and  W mag , both of which increase with as frequency increases. When the frequency is higher than 0.1 Hz,  W mag increases faster than Q c , leading to reduction of COP max with increasing frequency. When frequency is less than 0.1 Hz, the increase of Q c overcomes  W mag . As a result, COP max increases with increasing frequency. The corresponding U of the COP max increases as the frequency decreases, although it is smaller than that of maximum cooling power at the same frequency. This is due to the mismatch between the maximum values of Q c and  W mag . Because of this mismatch, it is impossible to obtain the highest Q c and COP in one configuration simultaneously. Therefore, after considering the priority of COP and comparing the 3 optimized point candidates, which are f = 0.167, f = 0.1 and f = 0.05 Hz at U = 0.1, U = 0.15, and  Fig. 4. Note that all investigated AMRs had the same length, height, width, and porosity, regardless of their number of layers. The only difference among them was the length of each MCM layer. In Fig. 4, the layer # represents the layer number of the MCM in the 16-layer regenerator, e.g. in the 4-layer regenerator, Layer # 1, 6, 11 and 16 indicate that the materials are the same as the MCM Layers 1, 6, 11 and 16 in the 16-layer regenerator. Simulations of AMRs with different N layer were conducted under the same conditions to compare their performances. Figure 5(a) depicts the Q c change with U for different AMRs with different N layer . As described above, for all the studied regenerators, Q c increases to a maximum then decreases with U, generating a maximum cooling power (Q c,max ) at about U = 0.4. As expected, the 4-layer AMR demonstrates the lowest cooling power, i.e the highest cooling power is less than 6 W from the 4-layer AMR. With the N layer increasing from 6 to 10, the cooling power is improved significantly, which reaches as high as 14.3 W for the 10-layer AMR. To clearly illustrate the relationship between Q c,max and number of layers in AMR, Fig. 5(c) depicts the Q c,max as a function of the numbers of layers. Figure 5(c) shows that Q c,max of the 6-layer AMR is less than the 8-layer one, while the 10-layer AMR provides the greatest Q c,max . Since the Curie temperatures of the MCM layers of AMRs fit the temperature distribution along the length of the regenerator, the shorter layer length leads to a narrower temperature span in one single layer of 10-layers AMR than AMR with less layers. This allows AMRs to take full advantage of the cooling ability of the individual MCM layers, leading to a higher cooling power. However, when N layer = 16, the maximum cooling power drops. As shown in Fig. 5 (c), the Q c,max of the 16-layer AMR is less than the 6-layer one. This is because going from 10-layer to 16-layer AMR, the arrangement of the MCM layers becomes denser. This means that the adjacent layers have only 2.5 or 2.25 K temperature difference for the 16-layer AMR, which may lead to a mismatch in the regenerator operating temperature. This mismatch in the regenerator operating temperature range  in 16-layer AMR may not provide the optimal total entropy change comparing to the AMR with fewer layers (e.g. 10-layer), which has been reported by Govindappa et al. 28 . This phenomenon was also observed by Lei et al. 34 .
Their results indicate that the cooling power does not monotonically increase with N layer . Since the properties of different layers in their simulation is identical (only shifting by Curies temperatures), it only occurs when the temperature span is low. However, due to the difference of entropy curve for every layer, the phenomenon observed in present work happens in a greater temperature span. Moreover, it implies that an optimal N layer exists (here is N layer = 10) for the AMRs in present work, that could provide the highest cooling power among the candidate options. Figure 5(b) shows the COP change with U for AMRs with different N layer , which illustrates that the COP from the 4-layer regenerator is always less than 1.7 for any utilization. When N layer is increased from 4 to 6, 8 and 10, improvement in COP in these regenerators is observed (COP max from 1.6 for 4-layer to 3 for 10-layer) due to the cooling power enhancement of the AMR with more layers. However, the COP max in these three regenerators never exceeds 3.1 as shown in Fig. 5(b). From the previous utilization discussions, it was concluded that a desynchronization exists between Q c,max and COP max as a function of U, i.e. a certain value of utilization cannot provide both the highest Q c and the highest COP simultaneously. This trade-off between COP max and Q c,max also can be observed for the performance variance as a function of N layer . Figure 5(b) shows that although N layer = 10 is the optimal number of layers for cooling power, it is not the optimum choice for COP; the16-layer AMR offers the highest COP (about 5.5) at U = 0.15. As discussed above, due to the competition between Q c and  W mag , the U for the COP max is always less than the U for the Q c,max , so COP max occurs at U = 0.15 or 0.2 for all the AMRs instead of U = 0.4, corresponding to Q c,max . It can be observed that in the range of U = 0.15 to 0.2, except the 4-layer AMR, the cooling powers of the different AMRs are very close. Therefore, the reason for the high COP max of the 16-layer AMR can be attributed to the low  W mag in the 16-layer AMR due to the significant variation of MCM properties (e.g. entropy and specific heat) with temperature. It can be concluded that the combined effects of the mismatch of Q c and  W mag as a function of U and the low  W mag in the 16-layer AMR lead to the highest COP max from 16-layer AMR. Also, it is not feasible to obtain a highest Q c and COP from an AMR simultaneously.  It can be concluded that significant COP enhancement was achieved by utilizing the 16-layer AMR. As mentioned above, the COP of the 16-layer AMR can reach as high as 5.5, which amounts to 84.6% of the Carnot cycle COP.
Optimization of the 16-layer parallel plates regenerator. Concerns over optimal distribution of MCM layers in multi-layer AMRs have stimulated studies to optimize the MCM layers 19,26,32 . It is essential to investigate the optimal distribution of the 16-layer parallel plates regenerator since it could provide guidance to achieve higher cooling power and COP. To accomplish this task, a DoE was employed to guide the parameters in numerical experiments and condense the essential number of numerical simulations. After the simulations were conducted following the DoE, a plot was generated to pick the cooling power and COP. Since the total length of the MCM layers is L, which is a constraint on the component proportions, a constrained mixture design was used to assess the influence of the length of every MCM layer on the overall regenerator performance 49 . Three levels of    L. Following the guidance of the DoE, a set of 137 simulations were conducted for f = 0.1 Hz and U = 0.15 obtained in Section 1, since the two parameters are the optimized working conditions for COP of all of the 137 runs. The cooling power and COP results from these simulations are compiled in Fig. 6(a), in which the squares and circles denote the results of COP and cooling power, respectively. The results are in the form of relative difference defined as ( where subscript i and ref indicate the ith run from DoE and the reference case, respectively, and R represents either Q c or COP. Figure 6(a) reveals that most of Q c from the DoE cases are close to the reference case as shown by the negative relative difference with small absolute values. Only a few runs return a higher cooling power, which are more than a 10% increase from the reference case. By examining the DoE parameters, it was discovered that every high cooling power result corresponds to the cases with a longer layer# 16, e.g. case 16,31,45,135 and 136 in Table 2. Note that layer# 16 is the last layer of MCM, which is located just next to the cold end of the regenerator, as shown in Fig. 1(a). Figure 6(b) depicts the temperature profiles along the MCM at the end of cooling cycle for case 16 and the reference case. It shows the MCM temperatures next to the hot end are very close between the two cases. However, next to the cold end, the temperature in the case 16 is much lower than the reference case, due to the longer layer# 16. Therefore, in the cooling cycle, a longer layer# 16 generates a lower temperature range next to the cold end as compared to the reference case, which is the reason for the 10% increase in the cooling power.
Moreover, Fig. 6(a) also presents the COP distribution of the DoE cases. It reveals that the values of COP of most of the DoE cases are about 10% less than the reference case. Considering that most of these cases yield smaller cooling power compared to the reference case, it can be concluded that the non-uniform layer lengths have a stronger effect on the magnetic work than on the cooling power, leading to higher drop of COP than cooling power. Figure 6(a) also reveals that when the non-uniform layer lengths are utilized in the 16-layer AMR, only trivial COP benefit is gained, i.e. the greatest COP enhancement from DoE cases is only 1.19% compared to the reference case. However, it can be identified from Fig. 6(a) that all of the cases with a longer layer# 16 result in a COP recession, e.g. the COP in case 16 is about 20% less than the reference case, while in case 31, a smaller increase in the length of layer 16 leads to an increase in Q c that is as great as the decrease in COP (about 10% for both). This is due to the competition between  W mag and Q c : the longer layer# 16 not only enhances Q c , but also requires more power to overcome the increasing  W mag . In the cases with a longer layer# 16, the increasing  W mag is more than Q c enhancement, resulting in a lower COP output. Therefore, although the longer layer# 16 provides a higher cooling power, the dramatic COP drop eliminates them from the list of optimization.
The DoE results can be summarized as: (1) Most of non-uniform layer cases lead to a lower Q c and a much lower COP compared to the reference case; (2) Although some cooling power benefits could be obtained from some cases with a longer layer# 16, all these cases result in a significant COP drop; (3) only a few combinations of non-uniform layer length are capable of providing COP benefits, which are no more than 1.2% compared to the reference case. Therefore, some optimization work may be required for the previous work with fewer layers MCM; however, this work suggests limited benefit of optimizing the length of every MCM layer for present 16-layer AMR case. This is because the Curie temperatures of the adjacent MCM layers of the 16-layer regenerator are very close to each other. As a result, there is no space to improve the COP.

Conclusions
A room temperature magnetic refrigerator with a 16-layer parallel plates regenerator has been numerically investigated in present work. The property data of these 16 MCMs with different Curie temperatures was derived from experimental measurements provided by the material manufacturer. The Curie temperatures were from 268.4 K to 305.4 K with temperature steps 2.25 or 2.5 K. A 1D model was employed based on the previous single layer AMR model, in which a new method was developed to accelerate the model to make it executable for the multi-layer AMR model. Results indicate that the computational time of a 16-layer case is very close to the single-layer case using the new method, while, without the new method, the time required is almost 6 times that of the single-layer case.
Using the new model, performance analysis of the 16-layer AMR was conducted between 41 K temperature span and a maximum applied magnetic field of 1.5 T. It was found that for every frequency, there exists a maximum cooling power for the 16-layer AMR. The Q c, max decreases as frequency increases. The same maximum COP with frequency variance is also observed except for the highest frequency case. The results reveal that around f = 0.1 Hz, the 16-layer AMR reaches its highest COP.
Additionally, a performance comparison between AMRs with different layers is provided in the present work. The results suggest that the AMR with a 10-layer regenerator offers the greatest cooling power. However, the AMR with a 16-layer regenerator can significantly enhance the system COP max , which could reach 84% COP of the Carnot cycle. These results guided an optimization investigation of layer length of the 16-layer AMR. A constrained mixture design, as a kind of Design of Experiments, was employed to generate a series of numerical experiments for the optimization. The design includes 137 simulations. According to the results, it is not essential to use a 16-layer AMR with non-uniform layer lengths, since it cannot find a significant COP enhancement by adjusting the layer lengths in DoE. This is explained because the Curie temperature between the adjacent layers is very close (2.25 or 2.5 K), and there is no room to improve COP by adjusting the space length of the layers.