Refrigeration Performance and Entropy Generation Analysis for Reciprocating Magnetic Refrigerator with Gd Plates

In the current work, a novel 2D numerical model of stationary grids was developed for reciprocating magnetic refrigerators, with Gd plates, in which the magneto-caloric properties, derived from the Weiss molecular field theory, were adopted for the built-in energy source of the magneto-caloric effect. The numerical simulation was conducted under the conditions of different structural and operational parameters, and the effects of the relative fluid displacement (φ) on the specific refrigeration capacity (qref) and the Coefficient of Performance (COP) were obtained. Besides the variations of entropy, the generation rate and number were studied and the contours of the local entropy generation rate are presented for discussion. From the current work, it is found that with an increase in φ, both the qref and COP followed the convex variation trend, while the entropy generation number (Ns) varied concavely. As for the current cases, the maximal qref and COP were equal to 151.2 kW/m3 and 9.11, respectively, while the lowest Ns was the value of 2.4 × 10−4 K−1. However, the optimal φ for the largest qref and COP, and for the lowest Ns, were inconsistent, thus, some compromises need be made in the optimization of magnetic refrigerators.


Introduction
Magnetic refrigeration near room temperature is a potential alternative to the traditional vapor compression counterpart [1][2][3]. An active magnetic regenerator (AMR), with the skeleton being manufactured by solid magneto-caloric materials (MCMs), is the key component of advanced room-temperature magnetic refrigerators. A typical magnetic refrigeration cycle consists of the four processes of magnetization, including heat rejection, demagnetization, and heat absorption, and its performance depends on a variety of factors. In addition to developing new MCMs [4,5], scholars and scientists have made great efforts to improve magnetic refrigerator performance from various aspects, such as efficient magnets [6][7][8], advanced refrigeration cycles [9,10], heat transfer enhancement [11][12][13][14][15], smart mechanical design [16][17][18][19][20], and optimal control [21,22]. On the other hand, in addition to building prototypes for experiments [23,24], various 1D and 2D numerical models, based on solving the energy conservation equations of fluid and solid MCMs, have been developed for the research of magnetic refrigeration [2,[10][11][12][13][14][15]25,26], and Nielsen et al. have made a review on the AMR modeling [27]. In these models the magnetic refrigeration fundamental was studied [10,12] and the effects of structural and operational parameters on performance aspects, such as temperature span, refrigeration capacity, and Coefficient of Performance (COP), were studied for optimization [28][29][30]. With the assumption of local thermal equilibrium between the fluid and solid phases, Teyber et al. developed a semi-analytic AMR model with two-layer MCMs [31], which could be extended for the optimization of AMRs with multi-layer MCMs.
As is well known, entropy generation takes place in an irreversible process [32], and the minimal entropy generation principle can be applied for the optimal design of various thermal equipment [33][34][35]. With the assumption of trivial fluid thermal capacity, Rowe and Barclay proposed an expression of temperature for optimal magneto-caloric effects, from the aspect of entropy generation minimization [36]. Li et al. [37] analyzed the irreversible losses in magnetic refrigerators and performed geometrical optimization for an AMR with Gd particles. Lei et al. [38] and Trevizoli et al. [39] conducted parameter optimizations for AMRs, with different geometries, based on the entropy generation minimization.
Compared with their counterpart of packed particles, the AMRs with MCM plates have the advantage of smaller flow resistance. A magnetic refrigeration apparatus, with an AMR of Gd plates, was built in the Risø DTU (Technical University of Denmark) [19] and a 2D numerical model of sliding grids was developed for it [26]. Both experiments and the numerical simulation of Risø DTU demonstrated that the refrigerator with MCM plates could generate a notable temperature span. To further improve the refrigerator performance of Risø DTU, thinner Gd plates, together with a smaller pitch, were adopted in our previous Computational Fluid Dynamics (CFD) simulation of a reciprocating refrigerator and the maximal temperature span was greatly improved. It is noted that the Navier-Stokes equations need be solved hundreds of thousands of times in the commercial CFD simulation of magnetic refrigerators and, thus, the computation is very time-consuming.
In the current investigation, a novel 2D numerical model based on stationary grids will be developed with Matlab for the reciprocating magnetic refrigerators with MCM plates. With the 2D model, numerical simulations will be performed to study the effects of structural and operational parameters on the refrigeration capacity and the COP of magnetic refrigerators. Moreover, the variations of the entropy generation rate and number with structural and operational parameters will be studied, and the contours of the local entropy generation rate will be presented for discussion. Figure 1a depicts the reciprocating magnetic refrigeration apparatus of Gd plates studied by the Risø DTU [19,26]. In their apparatus, the thickness (δ p ) and streamwise length (L) of the Gd plates were 40 mm of 0.9 mm, respectively, and the channel width, or the gap, between two adjacent plates (δ f ) was 0.8 mm. Five thermocouples were arranged to measure the fluid temperatures and the temperature span (∆T HC ) was determined by the subtraction of the 5th and 1st thermocouples under the no load condition. During the hot (or cold) blow, the two pistons on the two ends of the apparatus moved synchronously towards the right (or the left). The magnetic field (µ 0 H) at the center of the pole gap was~1.0 T, while the average field in the "out of field" position was~0.16 T due to the stray field of the electromagnet [19]. The current investigation is based on the above refrigeration apparatus. To enhance refrigeration performance, the magnetic field at the pole center had an intensity of 1.4 T and the AMR is demagnetized, with the field at approximately 0 T. Two thicknesses were adopted for the Gd plates, i.e., δ p = 0.4 or 0.8 mm, while the channel width (δ f ) was kept constant at 0.4 mm. Furthermore, the cycling period (τ) was 0.5, 1 or 4 s. Deionized water acted as the working fluid, and the hot and cold reservoirs were at the temperatures of T H (=300 K) and T C (=286 K), resulting in a temperature span of 14 K. The piston stroke is expressed by the relative fluid displacement (φ), which is calculated by:

Physical Model
Here, u p and τ/4 refer to the piston velocity and the duration time of the hot blow, respectively.

Governing Equations
During the magnetic refrigeration cycle, the solid MCMs were subjected to a magnetic field of varied intensity. With the consideration of the effect of the magnetic power, the refrigeration process could be expressed by the following fluid and solid energy conservation equations [26,27,35]: where T and t refer to the temperature and time, respectively, while the subscripts f and s represent the fluid and the solid, respectively. τ represents the viscous stress tensor and: τ = µ(∇u + ∇u T − 2 3 ∇ · uI).
As for the refrigerators manufactured by the MCM plates with a small pitch (see Figure 1a), the fluid flow between adjacent plates could be assumed to be 2D laminar and fully-developed, which indicates that the y-component velocity in a typical small channel equals zero, while the x-component counterpart is calculated by [26]: The q MCE of Equation (3) represents the energy source caused by the magneto-caloric effect of the MCMs. With the application of the thermodynamics equation of Maxwell, the energy source could be expressed by [11,26,27]: Here, the specific entropy (s) of the MCMs is calculated with the Weiss molecular field theory [25,40,41]. Lastly, the fluid and solid temperature fields in Equations (2) and (3) were coupled with the conservation of heat flux through the interface, i.e.,: where n represents the normal direction against the fluid-solid interface.

Computation Domain, Mesh Generation, and Solution Scheme
To minimize the computation load, the 2D unit structure, consisting of half a typical channel and half an MCM plate, was adopted as the spatial computation domain, and uniform grids were adopted for mesh generation, as depicted in Figure 1b. The spatial derivatives in Equations (2) and (3) were discretized with the central-difference scheme, while the scheme of fully-implicit forward differences was adopted for the unsteady terms. The software package of Matlab was adopted in the current computation. The solution scheme is similar to the scheme in You et al. [11], i.e., the discretized algebraic equations were solved with the sparse decomposition algorithm, and the whole solution procedure consisted of two iteration loops: The external one is for the periodical steady running, while the temperature fields of an entire cycle were obtained with the internal loop. The independence of the solution, on both spatial and temporal grids, were checked. With the compromise between the computation load and precision, the final computation of a typical case, with 0.8 mm plates, used 10 × 60 spatial grids together with 2400 time steps per cycle, and the iterative convergence criterion was set as 1.0 × 10 −6 .

Model Validation
The current specific entropy, specific heat, and adiabatic temperature increment of Gd, obtained by the Weiss molecular field theory, were compared with the experimental and numerical counterparts in the literature [19,26], and good agreement is observed among them.
The magnetic refrigerator in Ref. [19] was computed with the current 2D model for validation. The predicted maximal temperature spans (∆T max ) at various relative fluid displacements, obtained by the interpolations at zero refrigeration capacity, were consistent with their experimental counterparts in the literature. As the Weiss molecular field theory could overrate magnetic entropy variation and underrate specific heat near the Curie temperature [2], a mean overestimation of~2.17 K was observed in the current model validation.

Refrigeration Capacity and Coefficient of Performance
The refrigeration capacity and heat rejection per AMR volume, expressed by q ref and q rej , respectively, are calculated by [11]: The pressure drops of the hot and cold blows were calculated by empirical correlation. As the working fluid flowed into the smooth straight channel with a limited velocity, the pumping power consumption was trivial compared with the heat transfer rate and, thus, the Coefficient of Performance (COP), i.e., the ratio of the refrigeration capacity against the total power consumption could be calculated by:

Specific Entropy Generation Rate and Entropy Generation Number
As is well known, the entropy generation due to heat conduction can be calculated by s g,∆T = k|∇T| 2 /T 2 , thus, the entropy generation rate per AMR volume induced by the heat transfer in the fluid and solid MCMs, respectively, expressed by S f g,∆T and S s g,∆T , could be obtained by the integration of the local entropy generation rate over their corresponding computation domains. For example, the specific entropy generation rate of the fluid heat transfer is calculated by: Similarly, the local viscous entropy generation rate could be calculated by s g,∆p = µ where the tensors, u and x, have subscripts, i and j, and the viscous entropy generation rate per AMR volume could be obtained by: With the assumption that the irreversible losses of the magnetization and demagnetization were trivial, the specific entropy generation rate of the total AMR is equal to the sum of the fluid and solid counterparts, i.e., S AMR g = S f g,∆p + S f g,∆T + S s g,∆T . On the other hand, the total entropy generation rate could be calculated by the entropy balance equation, i.e.,: For a more rational comparison among different cases, the entropy generation number, normalized by refrigeration capacity, i.e., Equation (13), was adopted in the current study.

Variations of Refrigeration Performance with Relative Fluid Displacement
The variations of the specific refrigeration capacity and the Coefficient of Performance with relative fluid displacement under the cycling periods of 0.5, 1 and 4 s, respectively, are presented in Figure 2a-c, with the double vertical coordinates where the two plates' thickness (δ p = 0.4, 0.8 mm) were adopted in the AMRs.

Variations of Specific Refrigeration Capacity
From the left vertical axes of Figure 2a-c, it is observed that the two AMRs both generated the specific refrigeration capacity (q ref' ) to vary convexly with an increment of the relative fluid displacement (φ) for all the cycling periods (τ), which indicates that an optimal φ could be taken for the maximal q ref' . Moreover, the optimal φ increased with an increment of τ, and a moderate τ facilitated a larger peak of the q ref' . Specifically, under the conditions that τ equals 0.5, 1 and 4 s, the AMR, manufactured by the Gd plates of δ p = 0.4 mm, had the optimal φ of 0.125, 0.2 and 0.5, and the corresponding peak of the q ref' were equal to 120.8, 151.2 and 103.3 kW/m 3 , respectively. Moreover, it is found that the thinner Gd plates generated a larger q ref' , especially at a smaller τ. Specifically, when τ was 0.5, 1 and 4 s, the peak of the q ref generated by the AMR, with Gd plates of δ p = 0.8 mm, was about 43.9%, 69.1%, and 87.4% of the counterparts of the Gd plates of δ p = 0.4 mm, respectively.

Variations of the Coefficient of Performance
The right vertical axes in Figure 2a-c depict the variations of the Coefficient of Performance (COP) with relative fluid displacement (φ) under the cycling periods (τ) of 0.5, 1 and 4 s, respectively. Similar to the specific refrigeration capacity, the COP of the AMR follows the convex variation trend, with φ for all the τ. However, the largest peak of the COP was generated at the greatest τ. Furthermore, the AMR with thicker plates was found to go against a better COP, especially under the condition of a small τ. As for the AMR manufactured by the Gd plates of δ p = 0.4 mm, when the τ of 0.5, 1 and 4 s were adopted, the peak of the COP was 3.09, 5.67 and 9.11, respectively, while their optimal φ were equal to 0.075, 0.15 and 0.4, respectively. It is clear from Figure 2a-c that the optimal φ, corresponding to the maximal COP, was smaller than the counterpart of the maximal q ref .

Variations of Specific Entropy Generation Rates
The variations of specific entropy generation rates of the fluid heat transfer and the total AMR, expressed by S f g,∆T and S AMR g , respectively, with relative fluid displacement (φ) are depicted by the left vertical axes in Figure 3. Two thicknesses of the Gd plate (δ p = 0.4 and 0.8 mm) were adopted and Figure 3a-c shows the cycling periods (τ) of 0.5, 1 and 4 s, respectively. It is clear from Figure 3a-c that the S f g,∆T and S AMR g increased with the decreasing τ or rising φ. For example, for the AMR with δ p = 0.4 mm, the φ increased from 0.15 to 0.25, and the S f g,∆T and S AMR g at the τ of 1s varied from 31.5 and 38.4 to 55.5 and 71.2 W/(m 3 ·K), respectively. It is evident that these increments were related to the heat transfer enhancement due to a larger φ. Specifically, a larger φ could result in a larger heat transfer rate between the fluid and the solid and, thus, a greater entropy generation of the fluid heat transfer. Moreover, with the increment of φ, the irreversible loss of fresh water and the remanent fluid mixing in the channel may cause the viscous entropy generation rate to increase considerably. As for the solid Gd plates, the entropy generation rate and its increment with rising φ could be limited because it has a larger thermal conductivity. To explore the above phenomenon more clearly, the local mean entropy generation rates of the fluid and solid heat transfer (s f g,∆T and s s g,∆T ) were calculated by averaging the local entropy generation rates over their cross sections, and their variations over an entire cycle (τ = 1 s) are presented in Figure 4a,b and Figure 5a,b for the cases with the φ of 0.15 and 0.25, respectively. The log mean entropy generation rates based on two are adopted in Figures 4 and 5 for better distinction.
By scrutinizing the subfigures in Figure 4 or Figure 5, it is clearly seen that the fluid s f g,∆T was much larger than the solid s s g,∆T . Moreover, during the hot or cold blows, both s f g,∆T and s s g,∆T considerably increased with the increment of fluid displacement (φ). Contrastingly, as the fluid stayed stationary during the magnetization and demagnetization, the s f g,∆T and s s g,∆T in those processes were found to have no significant dependence on φ. All the contours are consistent with the curves of entropy generation rates in Figure 3.

Variation of the Entropy Generation Number
The variations of the entropy generation number (Ns), i.e., the total entropy generation rate normalized by the refrigeration capacity, with fluid displacement (φ) under the conditions with cycling periods (τ) of 0.5, 1 and 4 s are presented in Figure 3a-c, respectively, with the right vertical axes. It is clear from the right axes of Figure 3 that the Ns varied concavely with the increment of φ for all the τ values, and that thinner Gd plates facilitated a smaller pit Ns. Moreover, the smallest pit Ns were obtained by adopting the largest τ. Specifically, for the AMR manufactured by 0.4mm thick Gd plates, when the τ was 0.5, 1 and 4 s, the generated pit Ns were 9.1 × 10 −4 , 4.7 × 10 −4 and 2.4 × 10 −4 K −1 , respectively, and the corresponding optimal φ was 0.1, 0.15 and 0.3.
It is noted that with the increment of φ, the q ref and COP varied convexly, while the Ns changed concavely. Furthermore, the optimal φ for the largest q ref and COP, and the lowest Ns were inconsistent, and were 0.2, 0.4 and 0.3, respectively, in the current work. Thus, some compromises need be made in the optimization of AMR.

Conclusions
In the current work, a 2D numerical model of stationary grids was developed and numerical simulation was conducted for the reciprocating magnetic refrigerators with Gd plates under the conditions of various structural and operational parameters. By simulation, the effects of the relative fluid displacement (φ) on the specific refrigeration capacity (q ref ) and the Coefficient of Performance (COP) were obtained, and the variations of the entropy generation rate (S g ) and number (Ns) were researched. Furthermore, the contours of the local entropy generation rate were presented for discussion. From the current work, some conclusions were obtained: (1) The q ref and COP of the magnetic refrigerator followed the convex variation trend with increasing φ, and a smaller plate thickness (δ p ) facilitated a larger peak of the q ref and COP. Furthermore, the largest τ produced the largest COP (=9.11), while the greatest q ref (=151.2 kW/m 3 ) was generated at a moderate τ. (2) With increments of φ, the S g in the AMR rose monotonically, while the Ns varied concavely.
Moreover, a larger τ, or a smaller δ p, resulted in a smaller pit Ns, and the lowest Ns was 2.4 × 10 −4 K −1 . (3) The optimal φ for the largest q ref and COP and the lowest Ns was inconsistent, and were 0.2, 0.4 and 0.3, respectively, in the current work. Thus, some compromises need be made in the optimization of AMR.