Entropy Generation Minimization Analysis of Active Magnetic Regenerators *

A performance assessment of active magnetocaloric regenerators using entropy generation minimization is presented. The model consists of the Brinkman-Forchheimer equation to describe the fluid flow and coupled energy equations for the fluid and solid phases. Entropy generation contributions due to axial heat conduction, fluid friction and interstitial heat transfer are considered. Based on the velocity and temperature profiles, local rates of entropy generation per unit volume were integrated to give the cycle-average entropy generation in the regenerator, which is the objective function of the optimization procedure. The solidmatrix is a bed of gadolinium spherical particles and the working fluid is water. Performance evaluation criteria of fixed cross-section (face) area (FA) and variable geometry (VG) are incorporated into the optimization procedure to identify the most appropriate parameters and operating conditions under fixed constraints of specified temperature span and cooling capacity.

Regenerative cycles enable the achievement of large temperature spans and large cooling capacities.The Brayton cycle is the most common thermo-magnetic regenerative cooling cycle, and its T-S diagram is schematically presented in Fig. 1 (Rowe et al. 2005).The cycle starts with an adiabatic magnetization (1-2) of the magnetocaloric material, usually arranged as a porous matrix.In this process, the total entropy of the material remains constant, but the decrease in the magnetic contribution of the total entropy due to the magnetic field variation is compensated by an increase in the lattice and electronic entropy contributions.As a result, the temperature of the magnetic material increases, which is a manifestation of the MCE.The next step is a constant magnetic field cold blow (2)(3), where cold fluid (from the cold reservoir) flows through the matrix, cools down the solid phase (which releases heat to the fluid) and rejects heat to a hot reservoir.Next, an adiabatic demagnetization (3-4) of the material is performed and, similarly to the adiabatic magnetization, the demagnetization process reduces the temperature of the solid material, which undergoes an adiabatic temperature change.The final procedure is a constant magnetic field hot blow (4-1), where hot fluid (now from the hot reservoir) returns through the matrix, and the solid phase removes heat from the fluid phase, increasing the internal energy of the matrix.The fluid the absorbs heat from a cold reservoir, completing the cycle.The search for regenerator configurations and operating parameters capable of giving satisfactory values of temperature span and thermal efficiency at reasonable manufacturing costs are one of the main isuses of magnetic refrigeration research (Tura et al. 2012, Trevizoli et al. 2014a).
There is a direct relationship between the thermodynamic efficiency of a regenerative cooling cycle and the thermal effectiveness and viscous losses in the regenerator.The conflicting requirements of high heat transfer and low pumping losses can be equalized using thermal optimization methods based on Entropy Generation Minimization (EGM) (Bejan 1995).In AMRs, the structure and geometry of the solid matrix have to be optimized to reduce the thermal, viscous and magnetic losses and achieve the desired operating conditions of temperature span, cooling capacity and cycle efficiency.A concise review of applications of the method in the context of heat exchangers and storage systems was presented by Awad and Muzychka (2012).For passive regenerators (i.e., those without a magnetic-thermal interaction), several optimization procedures have been implemented (de Waele et al. 1997, Das and Sahoo 1999, Steijaert 1999, Nam and Jeong 2006, Trevizoli and Barbosa 2015).
In the context of AMR optimization, EGM analyses were conducted by Li et al. (2008) and Numazawa et al. (2012).Li et al. (2008) used a second-law analysis to find optimal regenerator dimensions and the diameter of the spheres in the packed bed for a real apparatus.They pointed out limitations regarding the regenerator housing diameter and the magnet circuit size.In order to keep the applied magnetic flux density constant and increase the housing diameter, the magnetic circuit volume has to increase exponentially.Thus, for a fixed housing diameter, an optimization of the housing aspect ratio was carried out, which resulted in an optimal coefficient of performance (COP) of 15.54 for a regenerator aspect ratio of 6.39 and 1.1-mm diameter spheres.Those simulations were carried out considering a hot reservoir temperature of 300 K and a 10-K temperature span.However, since Li et al. (2008) used the COP as the objective function, their analysis did not guarantee a fixed cycle-average heat transfer rate at the cold source for each combination of regenerator housing aspect ratio and sphere size.Numazawa et al. (2012) compared the performances of packed beds of spheres and parallel plates as regenerative geometries by changing geometric and operating parameters (plate thickness, sphere diameter, aspect ratio and frequency).However, as the cooling capacity was not kept fixed during the analysis, it was difficult to quantify the true effect of each parameter.
More recently, Aprea et al. (2013) carried out a detailed optimization analysis of a packed-bed AMR based on the second-law efficiency.They investigated the influence of the fluid properties, particle diameter, fluid blow time, mass flow rate, regenerator geometry and effect of axial thermal conduction.The mathematical model was written in terms of dimensionless parameters (e.g., number of transfer units, N T U and utilization) to facilitate the identification of the combination of variables that maximized the system performance.However, as in Li et al. (2008) and Numazawa et al. (2012), the cycle-average cooling capacity was not kept fixed for each combination of regenerator geometry and particle size, which may result in comparisons between systems that produce different values of cooling capacity.Ganjehsarabi et al. (2016) employed a genetic algorithm multi-objective optimization to determine the best design parameters of a cascade active magnetic refrigerator.The total cost rate of the system and the exergy (second-law) efficiency were used as objective functions.
In EGM methods, the total entropy generation can be used as an objective function in optimization procedures constrained by Performance Evaluation Criteria (PEC), such as those established by Webb and Kim (2005) for recuperators.For single-phase flow, a heat exchanger can be optimized according to fixed volume, fixed face area or variable geometry constraints.These constraints can be extended to the design of oscillating-flow regenerators for both passive and active applications.An explanation of the PEC employed in the present work is as follows: • Variable Geometry (VG): this PEC enables varying the regenerator housing cross-section area and length, while keeping the total housing volume constant.The VG PEC may be useful when there is freedom in terms of the final aspect ratio of the regenerator.For a given mass (volume) of magnetocaloric material, this PEC can be combined with the design optimization procedure of the magnet system in order to achieve an optimal configuration of the AMR device; • Fixed Face Area (FA): this PEC enables varying the regenerator housing length, while maintaining the cross-section area constant.The total mass of solid material is not constant in this PEC, since it changes with the housing length.In the context of magnetic cooling, this PEC can be useful in the optimization of the regenerator to be fitted in an existing magnet gap.In this case, the regenerator length can be adjusted to satisfy a specific design constraint or metric that may not consider the performance of the magnet itself, since its dimensions are already fixed.Conversely, the FA PEC can also be employed when neither the size of the magnet gap nor the mass of magnetocaloric material are constrained.In this case, both the magnet and the regenerator can be evaluated together in order to achieve the optimal applied magnetic flux density volume and the amount of raw magnetic material that satisfy a specific operating condition involving, say, the cooling capacity and the temperature span.
In order to achieve an optimal set of geometric parameters and operating conditions that return the minimum entropy generation in the regenerator, additional constraints must be established.For passive regenerators, as proposed by Trevizoli and Barbosa (2015), this additional constraint is imposed in terms of a target regenerator effectiveness, which must be met by the optimal geometry associated with each PEC.On the other hand, for AMRs, the proposed additional constraint is a fixed cooling capacity for a given temperature difference between the thermal reservoirs (temperature span), which are the desired output parameters of the actual cooling system.
This paper advances a performance analysis of AMRs based on the EGM theory.The mathematical model is composed of the 1-D Brinkman-Forchheimer equation for momentum transfer in porous media coupled with energy balance equations for the fluid and solid phases.The local instantaneous velocity and temperature fields are used in the calculation of the local rates of entropy generation per unit volume due to: (i) fluid friction, (ii) axial heat conduction in both media and (iii) interstitial heat transfer with a finite temperature difference between the phases.The total cycle-average entropy generation rate, Ṡg , was used as the objective function in an optimization procedure that evaluates the influence of parameters such as the mass flow rate, operating frequency, regenerator cross sectional area, housing aspect ratio, utilization factor and particle diameter according to the FA and VG PEC of Webb and Kim (2005).It is important to notice that, in cases where magnetic and thermal hystereses are absent and the MCE is reversible, there is no magnetic contribution to the cycle-average entropy generation.It is also clear that, by reducing the rates of entropy generation in the regenerator, the operating costs of the magnetic refrigerator as a whole are also reduced.However, a complete cost evaluation is outside the scope of this paper.Cost analyses were performed recently by other authors (Tura andRowe 2014, Bjørk et al. 2016), though without the appropriate restrictions to the optimization process (e.g., fixed cooling capacity and temperature span) used in the VG and FA PEC explored here.

-MATHEMATICAL MODELING
The AMR model comprises the momentum equation to calculate the fluid flow through the porous matrix and the energy equations for the fluid and solid phases.The MCE of the solid matrix was implemented via the discrete approach (Nielsen et al. 2011).Demagnetizing losses were implemented according to the procedure described in Trevizoli et al. (2012).The heat transfer fluid was treated as water and the solid material was assumed as commercial-grade (99.5% wt.) gadolinium (Gd), whose physical properties, namely adiabatic temperature change, specific heat capacity and magnetization were obtained from interpolating functions generated from experimental data.The reader is referred to Trevizoli (2015) and Trevizoli et al. (2016b) for a detailed presentation of the Gd physical property database.The mathematical model is described next, taking into account the following simplifying assumptions: one-dimensional, incompressible fluid flow, low-porosity porous medium (ε < 0.6) and absence of body forces (Trevizoli et al. 2014b, Trevizoli andBarbosa 2015).The mathematical model and the closure relationships for packed beds of spheres used in the present analysis were compared with experimental data for passive and active magnetic regenerators (Trevizoli et al. 2016a, 2016b, Trevizoli and Barbosa 2016).

-MOMENTUM EQUATION
After the simplifying assumptions listed above, the 1-D Brinkman-Frochheimer equation for momentum transfer in porous media is given by (Kaviany 1995, Nield andBejan 2006), where the term on the left is the macroscopic inertial force and those on the right are the pore pressure gradient, microscopic shear stress (Darcy term) and microscopic inertial force (Ergun inertial term), respectively.The permeability of the porous medium and the Ergun constant are given by K = ε 3 d p 150(1 − ε) 2 and c E = 1.75/(150ε 3 ) 0.5 , respectively (Ergun 1952, Kaviany 1995).

-Fluid phase
The energy equation for the one-dimensional, incompressible fluid flow in the porous medium is given by, An Acad Bras Cienc (2017) 89 (1 Suppl.) where the terms on the left are due to inertial (thermal capacity) effects and longitudinal advection, and those on the right are the transversal heat transfer term calculated using a convective heat transfer coefficient, axial conduction and viscous dissipation terms, respectively.k ef f f = εk f is the effective thermal conductivity of the fluid phase.D || /α f = 0.75P e is the longitudinal thermal dispersion in a bed of spheres, for P e >> 1 (Kaviany 1995, Koch and Brady 1985).P e = Re dp P r/2 is the Péclet number and Re dp is the Reynolds number based on the particle diameter (Kaviany 1995, Koch andBrady 1985).The interstitial convective heat transfer coefficient for a packed bed of spheres, , was calculated using the Nusselt number correlation of Pallares and Grau (2010).

-Solid phase
After the simplifying assumptions, the energy equation for the solid phase is given by: where the term on the left accounts for thermal inertia in the solid, and those on the right are due to interstitial heat convection and axial heat conduction, respectively.k ef f s is the solid effective thermal conductivity calculated using the Hadley correlation (Kaviany 1995, Hadley 1986).

-MAGNETOCALORIC EFFECT AND DEMAGNETIZING LOSSES IMPLEMENTATION
In an AMR model, the MCE can be quantified in terms of the adiabatic temperature change, ∆T ad (∆H, T ).For a magnetic material at an equilibrium state specified by an initial temperature T 0 and magnetic field H 0 , the total entropy is s 0 (H 0 , T 0 ).If a variation in the applied field, ∆H = H 1 −H 0 , is performed adiabatically, then the total entropy of the system remains constant, s 1 = s 0 , and a new equilibrium state is reached at s 1 (H 1 , T 1 ).The difference between the initial and final temperatures is the adiabatic temperature change, i.e., ∆T ad (∆H, T ) = T 1 − T 0 .If ∆H > 0, then ∆T ad (∆H, T ) is positive and, reversibly, when ∆H < 0, ∆T ad (∆H, T ) is negative.In materials with a continuous (second-order) magnetic phase transition, such as Gd, there are no hysteresis effects (magnetic or thermal) and the MCE can be treated as fully reversible.An extended discussion on the calculation of ∆T ad (∆H, T ) for magnetization and demagnetization and its impact on the regenerator performance is available in Trevizoli et al. (2012), Bahl and Nielsen (2009) and Nielsen et al. (2010).
In the present AMR performance analysis, an instantaneous (step) magnetic flux density change between 0 and 1.5 T was considered.Although this is the simplest waveform possible, it is difficult to be implemented in practice.An evaluation of the effect of the magnetic field waveform on the AMR performance was carried out in Trevizoli et al. (2014b).As mentioned above, the MCE was implemented via the so-called discrete approach, which is the simplest and most straightforward calculation method (Nielsen et al. 2011).It consists of a direct temperature increment (or reduction) during the magnetization (or demagnetization) step.Thus, for a one-dimensional regenerator domain: where the change in magnetic field, ∆H, was considered instantaneous and uniform along the regenerator length.
The demagnetizing losses were implemented considering a constant and uniform temperature.More information regarding the demagnetizing losses modeling and its effect on the AMR performance can be found in Nielsen et al. (2011) and Trevizoli et al. (2012, 2014b).Thus, for a homogeneous temperature solid, the effective magnetic field, H ef f , is the difference between the applied magnetic field, H app , and the demagnetizing field, H dem : where H dem = N D M .N D is the demagnetization factor and M is the magnetization.The effective magnetic field was calculated using a simplified method that consists of assuming a linear temperature profile (between T C and T H ) along the regenerator bed and calculating an average temperature for each finite volume in the regenerator.Eq.( 5) is used to calculate the local values of the effective magnetic field at each volume, which remain unchanged along the entire simulation.
The demagnetization factor of the regenerator matrix geometry was estimated using a correlation available in Coey (2010) for a bed of spherical particles: where N D,csg is the demagnetizing factor of the casing, which is a function of its geometric shape.For cylindrical geometries, the relationship due to Sato and Ishii (1989) can be used.Considering a cylinder of radius r, whose height is arbitrarily chosen as 2n r, the casing demagnetizing factor is given by N D,csg = n/(2n + 1).N D,sph = 1/3 is the demagnetization factor of a single sphere.

ENTROPY GENERATION
In the regenerator, entropy is generated due to heat transfer and viscous dissipation.The local rate of entropy generation per unit volume is given by (Steijaert 1999, Trevizoli andBarbosa 2015): where the first term on the right is the entropy generation rate per unit volume due to interphase heat transfer with a finite temperature difference.The second and third terms are the entropy generation rates due to axial conduction in the fluid and solid matrix, and the fourth term is the entropy generation rate per unit volume due to viscous friction.The cycle-average entropy generation rate in the regenerator, Ṡg , is defined as: In the present analysis, Ṡg , in W/K, is used as the objective function to be minimized in the regenerator optimization.In cases where the magnetocaloric effect is reversible, such as in materials with a continuous magnetic transition (e.g., gadolinium), the local rate of entropy generation due to the MCE is directly included in the interphase heat transfer.In this case, the ∆T ad resulting from the external magnetic flux density variation creates a finite temperature difference between the solid and fluid phases.Otherwise, if irreversibilities are present, such as magnetic and thermal hysteresis, additional terms have to be included in Eq.( 7) to compute these irreversibilities.However, these will not be necessary in the present paper, as Gd is the magnetocaloric material being considered.

-PERFORMANCE EVALUATION CRITERIA
The simulations for the minimum entropy generation calculation were carried out based on the PEC of Webb and Kim (2005).The first criterion is of Variable Geometry (VG), for which the regenerator housing cross sectional area, or the housing diameter, D h,Reg , and the length, L Reg , are allowed to vary, keeping a constant housing volume.The second criterion is of Fixed Face Area (FA), where the regenerator housing cross sectional area is kept constant and the regenerator length can vary, thus changing also the housing volume.The baseline (reference) housing geometry has the following geometric characteristics: D h,Reg = 25 mm and ζ = 2, where ζ = L Reg /D h,Reg is the aspect ratio of the regenerator housing.Thus, the baseline regenerator housing volume is 24.544 cm 3 , and the ranges of the variables explored in the analysis are presented in Table I.The operating conditions (i.e., frequency and mass flow rates) and ranges of geometric parameters evaluated in this paper were chosen so as to be consistent with figures currently encountered in magnetic regenerator systems (Nielsen et al. 2011, Yu et al. 2010).As can be seen from Table I, all of the possible D h,Reg and ζ combinations in the VG criterion resulted in a fixed housing volume of 24.544 cm 3 .The particle diameter was varied from 0.3 to 2 mm, for each possible combination.In the FA criterion, D h,Reg was fixed at 25 mm, while ζ was changed from 1 to 8. The particle diameter was also varied from 0.3 to 2 mm for each possible combination.
For each combination of D h,Reg , ζ and d p in each PEC, simulations were carried out considering two different scenarios.In the first scenario, the frequency was kept constant at 1 Hz and the mass flow rate was varied from 40 to 100 kg/h (steps of 10 kg/h) for PEC VG and from 50 to 200 kg/h (steps of 25 kg/h) for the FA PEC.In the second scenario, the frequency was varied from 1 to 4 Hz (steps of 0.5 Hz) and the mass flow rate was kept constant at 60 kg/h for the VG PEC and at 200 kg/h for the FA PEC.It is important to notice that the utilization factor, φ, can be constant or variable depending on the evaluation criterion (VG or FA), as presented in Trevizoli and Barbosa (2015).The utilization factor is defined as the ratio of the thermal masses of the fluid and solid phases: where c s and c p,f are taken as 350 J/(kg.K) and 3920 J/(kg.K), respectively, which correspond to reference values at 290 K (transition temperature of Gd).
In the VG cases, if the frequency is fixed the utilization factor is only a function of the mass flow rate because the housing volume and porosity do not change.Similarly, if the mass flow rate is fixed, the utilization factor is only a function of the frequency.On the other hand, for the FA PEC, the utilization factor decreases with ζ because the mass of the regenerative matrix increases with L Reg .Thus, for the FA PEC, φ varies with the mass flow rate (for a fixed frequency), or with the frequency (for a fixed mass flow rate) and with the housing aspect ratio.
In the present AMR optimization, the proposed performance constraint is of a fixed cycle-average cooling capacity for a fixed temperature difference between the hot and cold reservoirs.This constraint corresponds to a fixed cycle-average heat transfer rate at the cold heat exchanger, irrespective of the matrix size/mass and operating conditions.
The instantaneous cooling capacity, QC , is calculated by: and the cycle-average cooling capacity, QC , is determined by: where T CE is the time-dependent temperature exiting the regenerator at the cold end during the hot blow.
Considering balanced flow conditions (Schmidt and Willmott 1981), Qc can be approximated by the following relationship: where ṁ is the blow-average mass flow rate, T C − T CE is the blow-average temperature difference at the cold reservoir, ∆T C .
In the present analysis, a cycle-average cooling capacity, QC , of 20 W per regenerator was considered, for a temperature span, ∆T HEX , of 15 K, i.e., the hot reservoir temperature, T H , was set at 300 K and the cold reservoir temperature, T C , was set at 285 K.It should be noted that this choice of values is for illustration purposes only and has been chosen because it represents an operating condition that can be reached with several combinations of parameters.Other sets of performance constraints can be selected with qualitatively similar results.
For a specific operating condition (mass flow rate and operating frequency) the ranges of the geometric variables (e.g., housing and particle diameters and aspect ratio) were chosen according to each PEC (VG or FA), as shown in Table I.The temperature span and the hot and cold reservoir temperatures were kept constant for all operating conditions as part of the proposed constraint.For each condition in the range associated with a given PEC, the momentum and energy equations were solved numerically and the cycle-average entropy generation rate was calculated.
Dimensionless forms of the momentum and energy equations for the fluid and solid phases were solved using the finite volume method (Patankar 1980).The dimensionless form of the equations and the dimensionless parameters are presented in detail in Trevizoli (2015).The energy equations were implemented using a fully implicit scheme (in both the spatial and temporal terms) so that a coupled solution of these equations is performed at a given time.Since the momentum equation is not position-dependent, a fully explicit discretization scheme was adopted.Nevertheless, the Ergun inertial term contains a strong nonlinearity, which requires an iterative solution.The Weighted Upstream Differencing Scheme (WUDS) was used in the fluid energy equation and the Central Difference Scheme (CDS) was applied in the solid.The solver was based on a line-by-line method (Thomas Algorithm).The model was implemented in Delphi 7, with a user friendly interface.
In the AMR modeling, the solution of the energy equations can be strongly affected by variations in the solid and fluid thermophysical properties, especially the solid specific heat capacity.Thus, the dependence of the fluid and solid thermophysical properties on the temperature has been incorporated in the energy equations solver.For simplicity, however, this dependence has been omitted from the solution of the momentum equation.This influenced very little the heat transfer results, but significantly reduced the computational cost.
To accelerate the numerical convergence, the fluid and solid temperatures are initialized as linear profiles between T C at z = 0 and T H at z = L Reg .The initial condition for the momentum equation is u = 0 at t = 0.The boundary conditions for the fluid phase depend on the direction of the flow.If u > 0 (cold blow), then at z = 0 (cold boundary) T f is set equal to T C and, at z = L Reg , an outflow boundary condition is applied, i.e., ∂T f ∂z (t, z = L Reg ) = 0. On the other hand, if u < 0 (hot blow), at z = L Reg (hot boundary) T f is set equal to T H and, at z = 0, the outflow boundary condition is ∂T f ∂z (t, z = 0) = 0.The zero-derivative condition aims to represent a smooth continuation of the incompressible flow through the boundary.Its application is justified by the following: (i) the gradient changes in the flow direction at the boundary are expected to be very small; (ii) the conditions downstream of the outlet boundary have little influence on the flow upstream (thermal reservoirs); (iii) there is only one possible mass flow rate for the imposed instantaneous pressure gradient on the 1-D momentum equation.To solve the solid energy equation, the Central Differencing Scheme (CDS) was used as the interpolating function for the diffusion terms (Versteeg and Malalasekera 2007).The boundary conditions for the solid phase are: ∂Ts ∂z (t, z = 0) = ∂Ts ∂z (t, z = L Reg ) = 0.A numerical mesh consisting of 360 dimensionless time steps and 200 dimensionless control volumes has been used in all simulations.Based on a mesh study performed with several frequencies, mass flow rates, housing diameters, aspect ratios and particle diameters, the proposed mesh size was proven to be satisfactory (Trevizoli et al. 2014b, Trevizoli andBarbosa 2015).
For each scenario of the VG PEC, discrete data points were selected in the housing diameter and particle diameter ranges shown in Table I.The incremental steps were 2.5 mm for the range 12.5 ≤ D h,Reg ≤ 30 mm, 5 mm for 30 ≤ D h,Reg ≤ 50 mm and 0.1 mm for d p , which resulted in 216 independent cases for each value of mass flow rate or frequency.Similarly, for the FA PEC, the aspect ratio range was divided in incremental steps of 0.5, plus the 0.1-mm steps for d p , resulting in 270 cases for the variable frequency and variable mass flow rate scenarios for each value of mass flow rate or frequency.In total, 6804 different simulations were performed for the two PEC.As performed in Trevizoli and Barbosa (2015), the searches for the points of minimum entropy generation rate were refined further using 4th-order (or less) polynomial interpolations (R 2 > 0.9999) to guarantee the stability of the numerical solutions and provide a better resolution (finer than 0.5 mm for D h,Reg , 0.05 for ζ and 0.01 mm for d p ) for the minimum S g value at a reasonable computational cost.

-ENTROPY GENERATION VS. CYCLE-AVERAGE COOLING CAPACITY CONTOUR MAPS
As the numerical convergence of the governing equations is obtained for a specific set of simulated parameters associated with a PEC, the AMR configurations that yield the target value of Qc can be identified.As presented in Trevizoli and Barbosa (2015) for passive regenerators, Fig. 2 presents the contour maps for the VG PEC for ṁ = 60 kg/h and frequencies of 1 and 4 Hz.Fig. 3 shows the contour maps for the FA PEC for ṁ = 200 kg/h and frequencies of 1 and 4 Hz.Lines of constant Ṡg (red lines) are plotted together with lines of constant Qc (blue lines), as a function of the particle diameter and housing diameter for the VG PEC (Fig. 2) and particle diameter and housing aspect ratio for the FA PEC (Fig. 3).
From Figs. 2 and 3, it can be seen that the target value of Qc = 20 W per regenerator (blue solid line), i.e., the constraint, can be achieved with different combinations of d p and D h,Reg or ζ, for different frequencies and flow rates for each PEC.For the VG PEC with ṁ = 60 kg/h and f = 1 Hz, Fig. 2(a), Qc = 20 W per regenerator can be achieved with values of d p between 0.7 and 1.45 mm for housing diameters between 12.5 and 50 mm, which correspond to values of ζ between 16 and 0.25, respectively.In this case, the minimum value of Ṡg that guarantees Qc = 20 W per regenerator was around d p = 1.2 mm and D h,Reg = 20 mm.It is interesting to note that lower rates of entropy generation can be achieved with smaller particle diameters.For instance, a 0.4-mm diameter particle results in an entropy generation rate of approximately 0.005 W/kg.K.However, for the range of housing diameters shown in Fig. 2(a), the cooling capacities associated with d p = 0.4 mm are much greater than the (illustrative) 20-W constraint enforced to demonstrate the present analysis.This highlights the importance of correctly specifying the constraints of the analysis to avoid comparisons between cases that do not actually represent the same design condition.For an operating frequency of 4 Hz, Fig. 2 One may argue that the optimal particle sizes obtained above are somewhat large, given the current practice in AMR design and operation (Kitanovski et al. 2015).At these particular conditions, the optimal particle sizes were large for two reasons: (i) a small cooling capacity constraint was applied and (ii) external losses (such as void volumes and heat leaks) were not considered.Had a higher cooling capacity constraint been selected or had thermal losses been taken into account, a higher N T U matrix would be necessary to satisfy the operation constraints, which would certainly result in a smaller optimum particle size for a fixed mass flow rate.

-Fixed frequency
This section describes the results for the VG PEC for a fixed frequency and a variable mass flow rate.Performing an analysis similar to that presented in Trevizoli and Barbosa (2015) for passive regenerators, it is possible to find the optimal combinations of d p , D h,Reg and ζ that guarantee a Qc = 20 W per regenerator with ∆T HEX = 15 K.The results are shown in Fig. 4. A summary of the conditions associated with the points of minimum Ṡg is presented in Table II, which also presents the values of the utilization, φ, and the total demagnetization factor, N D,total -calculated from Eq.( 6) -which depends only on D h,Reg and ζ since the porosity of the regenerator is constant.Figs.4(a) and (b) are equivalent to the x-z and y-z projections of the conditions for which Qc = 20 W on the Ṡg surface shown in Fig. 2(a) for ṁ = 60 kg/h.The individual contributions to the total entropy at the minimum, Ṡg,min , due to interstitial heat transfer with a finite temperature difference, Ṡg,HT , axial conduction in the solid, Ṡg,SAC , axial conduction in the fluid, Ṡg,F AC , and viscous dissipation, Ṡg,V D , are presented in Fig. 5.
A comment is in order regarding the relative magnitudes of the individual contributions shown in Fig. 5. Conventional optimization logic dictates that optimal configurations should be encountered when the loss terms are more or less equalized.In our simulations, this was indeed the case near the overall minimum for each condition.For instance, for the condition shown in Fig. 2(a), i.e., VG PEC with ṁ = 60 kg/h and f = 1 Hz, it was observed that the individual contributions to Ṡg were similar in magnitude near the 0.4-mm particle diameter, where the entropy generation rate is close to 0.005 W/kg.K.However, this condition does not correspond to the 20-W cooling capacity design constraint used to illustrate the procedure.For this reason, the contribution due to interstitial heat transfer with a finite temperature difference dominates at the (local minimum) conditions of Fig. 5.An evaluation of the results in Figs. 4 and 5 and Table II reveals that for the smallest mass flow rate, ṁ = 40 kg/h, to achieve the performance constraint of Qc = 20 W per regenerator, a larger ∆T C is necessary (Eq.12), which is associated with a lower average temperature T CE (the reservoir temperature is fixed).To achieve a lower T CE , a higher effectiveness-N T U matrix is necessary, which is guaranteed by operating at the smallest utilization, with a small particle size (compared with 50-70 kg/h) and a longer matrix with a small D h,Reg .This results in higher superficial velocities.This optimum scenario is confirmed by the PAULO V. TREVIZOLI AND JADER R. BARBOSA JR  smallest absolute value of Ṡg,HT , meaning a highly effective heat transfer process.However, the use of a small particle size increased slightly the penalties associated with viscous dissipation.
As the mass flow rate is increased up to 70 kg/h, the optimal region shifts to bigger particle diameters, larger housing diameters and shorter aspect ratios (see Table II).Although the N T U is low, the small absolute values of utilization still result in a high effectiveness and a large ∆T C , which when multiplied by the low flow rate guarantees the 20-W constraint.This trend persists until a minimum N T U is found at ṁ = 70 kg/h.For mass flow rates higher than 70 kg/h, lower values of ∆T C are sufficient to guarantee the 20-W constraint.However, the increase in utilization needs to be compensated with a larger N T U in order to achieve the 15-K temperature span.Hence, the optimal conditions involve small particle diameters, large housing diameters and small values of aspect ratio.Since at high mass flow rates and small particle diameters the viscous losses become more important, a large housing diameter is required to decrease the superficial velocity, which results in shorter matrices.Nevertheless, smaller particle diameters are necessary to enhance the interstitial heat transfer, increasing N T U and decreasing Ṡg,HT (in %) (see insert in Fig. 5).Still, as the regenerator effectiveness decreases with the increasing utilization, the absolute value of Ṡg,HT increases  with ṁ.With shorter matrices, the axial conduction contributions become more important (in %) since the spheres size becomes smaller and the superficial velocity decreases.The contribution of Ṡg,SAC is small, around 1%. Ṡg,V D also increases with the mass flow rate and the reduction of the particle size.
Regarding the demagnetization factor, for a constant temperature difference between the thermal reservoirs, a bigger N D,total means a smaller effective magnetic field, H ef f -calculated using Eq.( 5) -and a smaller ∆T ad .To compensate the reduction of the MCE, a higher effectiveness regenerator is required.However, for the cases evaluated here, the variations in N D,total were little and, for the relatively small temperature span of 15 K used as a constraint, these effects were of secondary importance to the MCE when different optimal dimensions were compared.
The results presented in this section showed that the optimal configurations for f = 1 Hz and a variable mass flow rate are associated with small regenerator housing diameters (between 18 and 25 mm).Thus, the magnetic circuit size can be estimated considering a cylindrical Halbach permanent magnet array, for which the magnetic flux density, B, inside the bore of an infinite cylinder is given by Bjørk et al. (2008): where B rem is the remanence of the magnet grade and D e and D i are the external and internal diameters of the cylinder, respectively.Considering B = 1.5 T, B rem = 1.4 T for a Nd 2 Fe 14 B magnet grade, D i = D h,Reg and L mag = L Reg = ζD h,Reg , the optimal regenerator dimensions for each flow rate lead to the ideal Halbach cylinder dimensions presented in Table III.
The results in Table III show that the smaller mass flow rates require smaller diameters and longer magnets, while the higher flow rates are related to bigger diameters and shorter magnets.However, the total volume of permanent magnet necessary to create a flux density of 1.5 T in the bore of the Halbach cylinder is about the same for all cases.Therefore, in terms of the minimum entropy generated (less irreversibility), using the smallest mass flow rate possible is recommended for it satisfies the performance constraint with a more efficient AMR design.It should be recognized, however, that small flow rates are very difficult to realize in practice, as the system becomes very sensitive to heat leaks and carryover (dead volumes) losses.
As argued before, if external losses are included, the optimal values of d p are expected to be shifted to smaller figures and some of the operating conditions of low mass flow rate will no longer match the specific cooling capacity constraint.Thus, higher flow rates would be required in a more realistic scenario.

-Fixed mass flow rate
The results for the VG PEC for a fixed mass flow of 60 kg/h rate and frequencies between 1 and 4 Hz are presented as a function of changes in D h,Reg , ζ and d p for Qc = 20 W per regenerator at ∆T HEX = 15 K, as illustrated in Fig. 6.A summary of the conditions that led to the minimum Ṡg results is presented in Table IV and the individual contributions to the total entropy at the minimum, Ṡg,min , are shown in Fig. 7.In this scenario, with a fixed mass flow rate, ∆T C is also constant in order to achieve the performance constraint of Qc = 20 W per regenerator.As a result, the regenerator effectiveness is approximately constant and the output parameters vary according to the frequency/utilization.N D,total is also approximately constant.The results in Figs. 6 and 7 and Table IV show that for f = 1.0 Hz the optimal region involves comparatively smaller particle diameters, meaning a larger N T U for a higher utilization.Hence, larger housing diameters and shorter matrices are needed to decrease the superficial velocity and reduce the viscous losses, which leads to a more significant axial conduction contribution to Ṡg,min .
By increasing the frequency and reducing the utilization, the heat transfer is enhanced, which decreases the absolute values of Ṡg,HT , as seen in Fig. 7. On the other hand, the individual contributions of axial heat conduction and viscous dissipation at the minimum are reduced.As seen in Table IV, the longer matrices cause reductions in Ṡg,F AC and Ṡg,SAC , and bigger particle diameters are responsible for the reduction in Ṡg,V D .The optimum values encountered for D h,Reg , ζ and d p did not present significant variations with frequency, and, as a result, the N T U variation was also small.As seen in Fig. 7, the contributions of Ṡg,SAC are almost negligible, while Ṡg,HT dominates the total Ṡg,min .In conclusion, for the specified mass flow rate and for the cooling capacity and temperature span constraints presented here, operation at a higher frequency increases the level of irreversibility in the AMR, as Ṡg,min increases with frequency.The optimized points at ṁ = 60 kg/h and f between 1 and 4 Hz are associated with small housing diameters, around 19 mm, and aspect ratios of approximately 4.5.Regarding the magnetic circuit size, as in the previous scenario (fixed mass flow rate), all of the optimized conditions had virtually the same total volume of permanent magnet to generate a 1.5-T variation in magnetic flux density.The FA PEC is further explored in Fig. 8, which shows Ṡg as a function of d p and ζ.Table V summarizes the parameters associated with the optimal conditions for each mass flow rate.Fig. 9 shows the individual contributions to the total entropy at the minimum, Ṡg,min .
For ṁ = 50 kg/h, the minimum Ṡg configuration is out of the range of the particle diameter and aspect ratio.The values in Table V refer to the minimum value found (at the end of the interval).Thus, the ṁ = 50 kg/h condition is not considered in the following analysis.
For ṁ = 75 kg/h, the minimum Ṡg is associated with the smallest particle diameter, that corresponds to a large heat transfer area.Such a high effectiveness-N T U matrix is required to achieve the target cooling  capacity by increasing ∆T C , which results in a comparatively small value of Ṡg,HT (see Fig. 9).However, a shorter matrix is necessary to cut the viscous dissipation contribution, but it ends up increasing the penalties associated with the axial heat conduction.In this case, Ṡg,HT , Ṡg,V D and Ṡg,F AC have around the same weight (in % and absolute values) in the total Ṡg .
By increasing the mass flow rate, the minimum Ṡg is shifted to bigger particle diameters and longer matrices.As the mass flow rate is increased, ∆T C becomes smaller, and less effective matrices are neces-  sary, which raises Ṡg,HT .Therefore, the use of bigger particle sizes compensates the increasing superficial velocities and the effect of longer matrices on the viscous dissipation contribution, which increases in terms of its absolute value but not in relative terms (see insert in Fig. 9).Longer matrices help to reduce the impact of axial conduction.The increase of N D,total is also small for the optimized dimensions.
For a fixed frequency, Ṡg,HT and Ṡg,V D are the main contributions (in %) to the total entropy generation rate.As the mass flow rate increases, N T U decreases and ζ increases, making the heat transfer less effective (higher Ṡg,HT ) and, even though the particle diameter that minimizes the entropy generation increases, the viscous dissipation becomes more important.Also, the Ṡg,F AC and Ṡg,SAC contributions were found to be significant at lower mass flow rates, i.e., lower superficial velocities and shorter matrices.
Regarding the optimum regenerator and magnet dimensions, an external diameter D e = 73 mm was encountered for the ideal Halbach cylinder.Obviously, the amount of magnet material required to generate the 1.5 T flux density change increases linearly with the optimized regenerator length.In this case, in terms of the total minimum rate of entropy generation, again, using the lowest mass flow rate possible means a more efficient AMR design, which coincides with the smallest ideal magnetic circuit size (short Halbach cylinders).

-Fixed mass flow rate
This section presents the results for the FA PEC for a fixed mass flow rate of 200 kg/h and frequencies between 1 and 4 Hz.Again, the results are illustrated for a target Qc = 20 W per regenerator at ∆T HEX = 15 K. Fig. 10 shows the behavior of Ṡg as a function of d p and ζ.Table VI summarizes the optimal parameters for each frequency.Fig. 11 presents the individual contributions to the total entropy at the minimum, Ṡg,min .
As in the VG PEC, for a constant mass flow rate, the ∆T C that satisfies the performance constraint of Qc = 20 W per regenerator is also fixed.Hence, in the FA PEC for a constant mass flow rate, the  effectiveness changes only a little to compensate the small variations of N D,total and, thus, minor differences in the utilization factor and N T U are observed when compared with the FA PEC for a constant frequency.
For f = 4 Hz, the minimum Ṡg configuration is out of the range of the particle diameter and aspect ratio.The values in Table VI refer to the minimum value found (at the end of the interval).Thus, this case was not considered in the analysis below.
Since the thermal effectiveness is approximately constant, this AMR optimization scenario is quite similar to the passive regenerator optimization presented in Trevizoli and Barbosa (2015).For f = 1 Hz, the minimum Ṡg corresponds to a big particle diameter and a large aspect ratio, while for f = 3.5 Hz, the minimum Ṡg was characterized by a small particle size and a short matrix length.As the frequency is increased, the utilization factor decreases due to the shorter blow period.Nevertheless, the matrix length also decreases, reducing the thermal mass of the solid phase and the total interstitial area.Consequently, the utilization factor and the N T U did not present significant variations with frequency, resulting in a nearly constant contribution (in % and absolute value) of Ṡg,HT to Ṡg,min , as can be seen in Fig. 11.At low frequencies, Ṡg,V D was large even for bigger particle diameters due to the large value of the matrix length.On the other hand, higher frequencies can sustain shorter matrix lengths and smaller particle diameters, but this increases (in % and absolute value) the penalty associated with axial heat conduction in the fluid, Ṡg,F AC .Ṡg,SAC presented a minor contribution to Ṡg,min .
In the present scenario, for the specified mass flow rate and the cooling capacity and temperature spans constraints considered, the total minimum entropy generation rate was observed at the lowest operating frequency, but Ṡg,min slightly increases with frequency.On the other hand, higher frequencies would require a smaller ideal magnetic circuit size due to the shorter aspect ratios.Therefore, in this case, since the irreversibilities are not strongly affected by the increasing frequency, it may be a good idea to consider the magnet cost in the selection of the operating frequency.

-CONCLUSIONS
This paper proposed a procedure to quantify the performance of active magnetic regenerators (AMRs) based on the Entropy Generation Minimization (EGM) method (Bejan 1995) combined with performance evaluation criteria (PEC) of variable geometry (VG) and fixed face area (FA) (Webb and Kim 2005).A onedimensional model was developed to solve the fluid flow and coupled heat transfer in regenerative matrix composed of spherical Gd particles.A step change in the applied magnetic field between 0 and 1.5 T was considered.Entropy generation contributions due to axial heat conduction, fluid friction and interstitial heat transfer were taken into account in the mathematical model.
To perform a meaningful evaluation that enables a comparison of the system performance at different operating conditions, appropriate constraints had to be applied.Here, the optimal regenerator configurations were determined for a fixed average cooling capacity of 20 W per regenerator bed and a fixed temperature span between the thermal reservoirs of 15 K.However, virtually any choice of fixed cooling capacity and temperature span could be made.The constraints used to illustrate the present results gave enough room to span different mass flow rates and frequencies.Had a larger cooling capacity per regenerator been chosen, perhaps no results could be generated for low mass flow rates.
An important point regarding the search for the point of minimum entropy generation rate is that, in the present paper, the preferred approach was to evaluate the mass flow rate and the frequency separately.However, in designing a real device, mass flow rates and frequencies can be varied concertedly to find the minimum entropy production rate for the desired set of constraints.
The main practical conclusions arising from this work are as follows: 1.For the fixed frequency/variable flow rate VG PEC, the optimal AMR configuration for low flow rates involves a long regenerator with a small housing diameter (large aspect ratio) and larger particles, which act together to increase the interstitial heat transfer coefficient.At high flow rates, the optimal regenerator configurations are achieved with a small particle diameter and small aspect ratios to compensate the increasing viscous losses.Irrespective of the operating condition, interstitial heat transfer and axial heat conduction in the fluid presented the largest contributions to the entropy generation in the regenerator; 2. For the fixed flow rate/variable frequency VG PEC, the optimal AMR configurations for small frequencies are associated with smaller particle diameters and shorter regenerator lengths.This reduces the fluid superficial velocity and controls the viscous dissipation.At higher frequencies, the constraints of fixed average cooling capacity and temperature span can be satisfied by using longer matrices (larger values of aspect ratio) and bigger particles.However, the dimensions of the optimal regenerative matrix are much less affected by the operating frequency than by the flow rate; 3. For the fixed frequency/variable flow rate FA PEC, the optimal AMR configuration for low flow rates is associated with a short matrix and a small particle size.As the flow rate increases, the fixed cooling capacity and temperature span constraints are satisfied by larger particle sizes.As in the passive regenerator case (Trevizoli and Barbosa 2015), longer matrices are required to compensate for the decrease in insterstitial area.At low flow rates the contributions of the terms due to interstitial heat transfer, axial heat conduction in the fluid and viscous dissipation to the total entropy generation rate are all of the same order of magnitude.As the mass flow rate increases (longer matrices), the axial heat conduction becomes less important; 4. For the fixed flow rate/variable frequency FA PEC, the optimal AMR configurations for the lower frequencies have large regenerator aspect ratios and large particle diameters.As the frequency increases, smaller particles and shorter matrices are capable of satisfying the constraints of fixed average cooling capacity and temperature span.Interstitial heat transfer and viscous dissipation are the major effect contributing to the entropy generation in the matrix.

Figure 1 -
Figure 1 -Schematic T -S diagram of the thermo-magnetic regenerative Brayton cycle.
(b), the minimum Ṡg is shifted to larger values of d p , yielding other combinations of D h,Reg and ζ.For the FA PEC with ṁ = 200 kg/h and f = 1 Hz, Fig. 3(a), the minimum Ṡg for Qc = 20 W per regenerator was identified in the vicinity of d p = 1.1 mm, with an aspect ratio of around 4. As expected, for other frequencies the minimum Ṡg range changes and Qc = 20 W per regenerator is achieved with other combinations of d p and ζ, as presented in Fig. 3(b).

Figure 2 -Figure 3 -
Figure 2 -Lines of constant Qc (blue lines) and Ṡg (red lines) as a function of dp and D h,Reg for the VG PEC with ṁ = 60 kg/h and frequencies of (a) 1 Hz; (b) 4 Hz.

Figure 4 -
Figure 4 -Minimum entropy analysis for the VG PEC for Qc = 20 W per regenerator and ∆THEX = 15 K, for a frequency of 1 Hz and variable flow rates in the range of 40-100 kg/h: (a) Ṡg as a function of dp; (b) Ṡg as a function of D h,Reg ; (c) Ṡg as a function of ζ.

Figure 5 -
Figure 5 -Individual contributions to the total entropy at the minimum, Ṡg,min, for the VG PEC for Qc = 20 W per regenerator and ∆THEX = 15 K, for a fixed frequency of 1 Hz and flow rates between 40 and 100 kg/h.

Figure 6 -
Figure 6 -Minimum entropy analysis for the VG PEC for Qc = 20 W per regenerator at ∆THEX = 15 K, fixed mass flow rate of 60 kg/h and variable frequencies in the range of 1-4 Hz: (a) Ṡg as a function of dp; (b) Ṡg as a function of D h,Reg ; (c) Ṡg as a function of ζ.
5.3 -FIXED FACE AREA (FA) EVALUATION CRITERIA5.3.1 -Fixed frequencyThis section presents the results for the FA PEC for a fixed frequency of 1 Hz and mass flow rates between 50 and 200 kg/h as a function of ζ and d p .As mentioned previously, D h,Reg is constant and kept fixed at 25 mm.Again, the results are illustrated for a target Qc = 20 W per regenerator at ∆T HEX = 15 K.

ContributionsFigure 7 -
Figure 7 -Individual contributions to the total entropy at the minimum, Ṡg,min, for the VG PEC for Qc = 20 W per regenerator at THEX = 15 K, at a fixed mass flow rate of 60 kg/h and frequencies between 1 and 4 Hz.

Figure 8 -
Figure 8 -Minimum entropy analysis for the FA PEC for Qc = 20 W per regenerator at THEX = 15 K, fixed frequency of 1 Hz and variable mass flow rates in the range of 50-200 kg/h: (a) Ṡg as a function of dp; (b) Ṡg as a function of ζ.

Figure 9 -
Figure 9 -Individual contributions to the total entropy at the minimum, Ṡg,min, for the FA PEC for Qc = 20 W per regenerator at ∆THEX = 15 K, fixed frequency of 1 Hz and variable mass flow rates in the range of 50-200 kg/h.

Figure 10 -FrequencyFigure 11 -
Figure 10 -Minimum entropy analysis for the FA PEC for Qc = 20 W per regenerator at ∆THEX = 15 K, fixed mass flow rate of 200 kg/h and variable frequencies in the range of 1-4 Hz: (a) Ṡg as a function of dp; (b) Ṡg as a function of ζ.

TABLE I Ranges of the geometric variables for each PEC.
PEC D h,Reg [mm] ζ Housing volume [cm 3 ] d p [mm]