Development and Validation of a 3-D Analytical Fluidic Model in Cartesian Coordinates for a Magnetic Refrigeration Application

This paper focuses on the development and validation of a three-dimensional (3-D) analytical model based on the formal resolution of the incompressible Navier-Stokes equations in Cartesian coordinates. This analytical fluidic model calculates the evolution of the fluid flow velocity for various pressure gradient shapes (i.e., square, trapezoidal, ramp, triangle, sinusoidal signals, etc.), with or without mean value, using Fourier series. The model can consider square and rectangular channels of any dimensions in terms of height, width, and depth. It allows for the study of different types of fluid flow and geometries. The introduction highlights the significance of analytical modeling, particularly in the context of Magnetic Refrigeration (MR) technologies. The model assumptions and governing equations are presented. The results obtained from the analytical model are validated and show good convergence with those obtained using the COMSOL Multiphysics® software. The comparison demonstrates a high level of agreement between the two models. The analytical model performs well in terms of computation time and accuracy. It will be coupled with an analytical thermic model and used to optimize Active Magnetic Regenerative Refrigeration (AMRR) cycles. The ultimate objective is to determine an optimal pressure gradient evolution in order to maximize heat exchange in a regenerator channel.


Introduction
The field of analytical modeling holds significant importance in the development and optimization of new technologies.In this context, MR emerges as a promising alternative to conventional vapor compression systems, particularly at room temperature.MR offers a high-performance potential characterized by superior thermal and thermodynamic efficiencies [1].Analytical models provide a valuable tool for understanding and predicting the behaviour of complex systems, it becomes crucial to explore their application in the context of MR.This technology uses MagnetoCaloric Materials (MCM) and water as coolant to achieve heat exchanges and a temperature gradient between two heat sources.Unlike traditional vapor compression systems that use hydro(chloro)fluorocarbons as refrigerants, MR systems are environmentally friendly.To improve the performances of AMRR, one of the main problems is to find MCMs with a higher magnetocaloric effect.Numerous studies have focused on this issue [2].Nowadays, the magnetocaloric effect, which describes the change of the MCM temperature in response to the external magnetic field evolution, has been widely described and investigated in many papers [3].Innovative physics modelling and optimization are necessary to improve operating frequencies, heat exchanges, and achieve higher magnetic field generation.
There are many reviews on AMRR concerning modeling, applications, perspectives and studies on MCM [4]- [6].A comprehensive review of models has been published by Silva et al. [7] (2021), which includes a classification and highlighting of optimization studies for different types of systems.Eustache et al. [8] (2021) published a detailed review focusing on (semi-)analytical models used to study AMRR cycles.This review emphasized the lack of purely analytical and 3-D models, as the majority of the presented models were one-dimension (1-D) and mainly focused on thermo-fluidic aspects, neglecting the magnetic factor.Only few researches concern the magnetic behaviour of the MR system [9]- [11].Last year, Plait et al. [12] (2022) developed a complete two-dimensional (2-D) semi-analytical multiphysics model (magneto-thermo-fluidic).Computation time, an important parameter, is rarely discussed in studies, although [8]- [13] have touched upon this topic.(Semi-)analytical models are generally faster than the numerical ones, which often require supercomputers to avoid excessively long computation time.An interesting analytical fluidic model was developed by Fan and Chao [14] (1965).They solved the Navier-Stokes equations for several forms of pressure gradient.This type of analytical model is less frequent now in the literature.In fact, this approach has largely been abandoned in favor of numerical models, which are more convenient to develop and utilize.
The analytical model presented in this paper focuses on the velocity modeling of an oscillating fluid flow within a channel of a parallel plate regenerator.This model represents the first step towards analytical and systemic modeling of AMRR, which combines magnetic, thermal and fluidic physics.It determines the magnitude and velocity profile in a channel as a function of the imposed pressure gradient.The user defined the pressure gradient, determined using Fourier series, serves as an input parameter and initial condition for the model.The evolution of velocity over time is directly influenced by this pressure gradient, and the shape of the velocity profile is determined by the oscillation frequency of the fluid.This analytical model can serve as a tool to understand and simulate the oscillating flow of a fluid in a steady-state behaviour.The second section of the paper introduces the assumptions and geometry of the model, establishing a framework for understanding the resolution methodology and the operational principles of the model.Then, the third section highlights the obtained solutions and explains the steps taken to derive the results.Finally, in the fourth section, the model's validity is assessed by comparing its outcomes with those obtained using a numerical software, in order to evaluate its performance and accuracy.

Considered geometry and assumptions
Figure 1 represents a schematic view of a parallel plate magnetocaloric regenerator positioned between two heat sources.A mechanical system, such as a reciprocating hydraulic cylinder, is employed to induce oscillating fluid motion between both heat sources.The model objective is to determine the  velocity profile of the fluid moving between fixed parallel plates.Aprea et al. [15] (2017) have demonstrated that parallel plate geometry is highly effective in enhancing heat exchange between the fluid and the solid.Figure 2 illustrates the geometry considered for the study.
It has been shown by Lionte et al. [16] (2015) that the fluidic effects at the inlet of a channel are negligible compared to the channel length.In the studied case, the fluid is assumed to be incompressible, and all channels have the same dimensions and are long enough to neglect the inlet effects.The fluid flow profile is already quasi established at the channel inlet, and the pressure is homogeneous in all sections of the regenerator.Under these assumptions, the third dimension can be neglected in the equations.Consequently, the study can be limited to a single channel, since the velocity profile and fluid behaviour are similar in each channel.The following assumptions are applied to the system of Navier-Stokes' equations, where all bold characters represent vectors with components in three dimensions: • The fluid has a laminar flow, so the velocity vector has only one component along the z-axis (i.e., the velocity vector  = {0; 0; u z }, the value of the velocity still depends on the spatial and temporal dimensions, so u z = f(x, y, t)).• The only active force on the fluid is a pressure gradient created by a pressure difference between channel inlet and outlet.The channel filled with water, is very small and gravity is therefore negligible.• The pressure gradient generates a steady state flow velocity along a channel of constant crosssection, so a 2-D analytical resolution is sufficient to produce a 3-D model.• Since the fluid is considered incompressible, the velocity divergence is equal to zero (i.e., ∇.  = 0).These assumptions simplify the solution of the Navier-Stokes' system of Partial Differential Equations (PDE).The Boundary Conditions (BCs) are: u z (x i , y, t) = 0, u z (x f , y, t) = 0, u z (x, y i , t) = 0 and u z (x, y f , t) = 0. To ensure that the fluid flow remains laminar and don't become turbulent it is important to control the value of the Reynolds number (i.e., noted Re): with L the characteristic length (m) depending on the geometry, u the fluid flow velocity (m s ⁄ ) and ν = μ ρ 0 ⁄ the kinematic viscosity (m 2 s ⁄ ) where μ and ρ 0 are respectively the dynamic viscosity (Pa • s) and the mass density (kg m 3 ⁄ ).There is a critical Reynolds numbers beyond which the fluid flow is considered turbulent.Beyond this value, the analytical model is no more reliable [17].

Governing PDE
The Navier-Stokes' equations are a set of equations that the conservation of mass, momentum, and energy for a Newtonian fluid.By solving these equations, it's possible to calculate the flow velocity profile of a fluid.( 2) presents the Navier-Stokes' system of PDE: with p the pressure (Pa), and p the source term of the model which is the pressure gradient.Applying the assumptions, the system of equations to be solved reduces to: where x, y and z are the spatial coordinates as shown in Figure 2 and t represents time.To simplify the study, the time variable t is replaced by the electrical angle ϑ = ωt, where ω = 2πf is the electrical angular frequency (rad s ⁄ ) with f the electrical frequency (Hz).This adjustment allows us to simulate results over a period from 0 to 2π, independent of the fluid's oscillating frequency.(3) thus becomes: where δ is the skin depth defined by: This affects the performance of the numerical simulation and the physical phenomena.This point will be discussed further in the model validation section.

Pressure gradient
The second element of the governing PDE is the pressure gradient.It can be defined as the product between the pressure amplitude difference and a time pressure signal, representing the evolution of the pressure gradient: The time pressure signal p z (ϑ) is a dimensionless function of ϑ, which can be defined by a Fourier series as (7.a).Since the cross-section of the regenerator remains constant (see Figure 1 and Figure 2), the pressure varies linearly along the channel between the inlet and the outlet.Therefore, the pressure evolution in the channel depth p z (z) can be defined as an affine function given by (7.b).
For p z (z), p out and p in represent respectively the pressure at the channel outlet and inlet.For p z (ϑ), p 0 and p n represent respectively the mean value and the harmonic term with n the time harmonics.These various terms depend on the type of pressure gradient studied.Two Fourier series has been considered to generate different pressure gradient signals, with and without mean value (i.e., Figure 3.a and 3.b, respectively).Pulsating flows can be studied by adding a mean value to the Fourier series, otherwise oscillating flow are studied.Moreover, the signal form, such as square, triangle or sinusoidal, can be studied by adjusting the values of ϑ s and ϑ rc or by changing the number of harmonics n in the summation.Finally, by replacing the pressure gradient term (4) becomes:

Analytical solution
This section highlights the main equations used in the model resolution.The velocity can be expressed as a Fourier series, similar to the form of the pressure gradient acting as the source term.8) can be solved in two steps due to the influence of the second member.A first step for the mean value of the Fourier series and a second step considering the harmonics sum of the time.Thus, by identification with the source term, the velocity solution takes the following form: where u z0 and u zu are respectively the mean value and the harmonic term.To find the final solution, the general solution and the particular solution for the mean value and the harmonic term must be determined and added.However, some unknown coefficients will still remain.The application of the Dirichlet BCs allows the determination of these remaining unknown coefficients.Furthermore, a verification process is carried out at each step of the resolution.For example, the final solution obtained for u z is substituted back into 8 ) to confirm that equality is still respected.

Solution for the mean value
The resolution of (10) allows to study the behaviour of the fluid with a non-zero average velocity over one period.This kind of fluid flow are called pulsed flow.
The general solution (11.b) is determined after using the separation of variables method for a second member equal to zero.The particular solution (11.c) is the unique solution taken by u z0 to respect the equality of (10).
where k are the spatial harmonics in the x-axis, λ k = k • π τ x ⁄ is the spatial periodicity in the x-axis, and I 0, is a coefficient obtained after integration by Fourier series expansion with BCs (i.e., see (18) in Appendix).The geometric coefficients τ x = x f − x i and τ y = y f − y i represent the height and width dimensions, respectively (see Figure 2).(11.a) respects the BCs, u z0 is cancelled by replacing x by x i or x f and y by y i or y f .

Solution for time harmonics
The resolution of (11) allows to study periodic fluid flow.
The same method is applied to obtain the solution with time harmonics.The solution for ( 12) is the following: u zn = u zn,k + u zn,0 (13.a) where λ n,k = √ λ n 2 + λ k 2 , with the time periodicity term λ n = j • n • 2 δ 2 ⁄ and I , is an harmonic dependant term determined by using Fourier series expansion (i.e., see (19) in Appendix).

Final notation
The final solution ( 14) is obtained by replacing all the previous term by their value in (9): The above equation and all its corresponding terms, have been written under a Matlab® code on the R2019a version.Axis values and the harmonics terms are defined as vectors.The results obtained were with those obtained with a numerical model by extracting the results in Matlab® in matrix form.The results are subsequently plotted and juxtaposed compared in the following section.

Numerical model
The analytical results were compared against numerical results obtained using the finite-element method on COMSOL Multiphysics® software.This software can be used to solve multiphysics problems in 2-D or 3-D.In order to validate our analytical model, two types of channel geometry have been modelled, for different frequencies (i.e., 0.01 Hz, 0.1 Hz, 1 Hz, 5 Hz and 10 Hz) and different pressure gradient shapes (i.e., trapezoidal and sinusoidal).The fluidic problem was defined in COMSOL Multiphysics® as a transient single-phase laminar flow.For each simulation, the computation time and the obtained results were compared in order to show the accuracy and performance of the analytical model.The studied fluid motion represents an oscillating fluid flow without a mean value.The results are presented in absolute values.For ϑ varying from 0 to π, the fluid flows from the inlet to the outlet and in the opposite direction for ϑ varying from π to 2π.The fluid flow direction is imposed by the pressure gradient.
For the comparison, the numerical model must be created in 3-D.Two geometries have been studied to verify the performances of the model for several channel sizes.The first channel geometry studied is similar to that one usually found in the conventional parallel plate structures of magnetocaloric regenerators.This geometry is a thin rectangular channel for a length of 20 mm.The flow cross section is (1x13) mm².The second geometry is also a (13x13) mm² square channel of the same length as the initial geometry.
The same BCs were applied to both models.However, the numerical model requires initial conditions for velocity and pressure.Inlet and outlet conditions are also required.For the inlet, a pressure waveform with the half amplitude of the one used in the analytical model was applied.Similarly, the same waveform was used for the outlet, but with a half-period offset.

Geometry and frequency influences
In COMSOL Multiphysics®, several mesh levels are available (i.e., from extra-coarse to extra-fine).For the same level of mesh accuracy, the number of meshes is greater for thin dimensions.However, using a very high mesh accuracy leads to a significant increase in computation time.Having a mesh accuracy that exceeds what is considered "finer" can lead to problems such as "out of memory" or excessively long computation times, especially on less powerful computers.This is clearly visible in the results comparison presented in Table 1, where the model exhibits longer computation times for the thin width channel geometry compared to the square channel.The necessity for a more precise mesh in thin geometries is due to the skin effect of the fluid [18], similar to the skin effect observed in electrical conductors [19].This effect becomes more important at higher operating frequencies, as indicated in (5) due to the smaller depth.The finite-element method requires a minimum number of nodes in a dimension equal or smaller than the skin depth to ensure a minimal difference between the model and reality.The analytical resolution does not encounter these issues.
For low frequencies and small cross-section channels, the steady-state is reached more quickly, so fewer periods need to be calculated to compare the two models.However, in the case where a larger   volume of fluid is transported, a higher number of simulation periods will be required to achieve a steady-state for an oscillating flow.Furthermore, the oscillation frequency of the fluid flow will also influence the fluid flow profile.It has been observed that above a certain value of the dimensionless Womersley number [20], denoted as α, an annular effect appears (as shown Figure 4).This phenomenon has been described by Richardson and Tyler [21] (1929) when they studied the velocity of an alternating flow.The Womersley's number is defined by: where D H is the characteristic length (i.e., the hydraulic diameter of a circular pipe).The interest of MR for this annular effect concerns a potentially better heat exchange.This issue will be investigated in a future study with a thermal model coupled with the presented analytical fluidic model.

Model validation
To compare the models, two equations have been used to calculate the difference between the obtained results.The first is the relative error calculated with The relative error equation has been applied to compare the velocity amplitude.For laminar flow regime, the velocity amplitude is always more important at the channel center, considering that the speed of the fluid is zero at the wall.Figure 5 compares the velocity evolution at the center of the channel for two different geometries and for the same frequency.The results show good agreement with the maximum velocity amplitude with a rather small relative error.(17) represents the MSRE used to compare the difference between both models for each point on a surface over a whole period as shown in Figure 6.
where u a is a vector of velocity analytical value and u n a vector of velocity numerical value.Both models return the results in matrix form.To compare the values, the matrices are transcribed into vector and s u is the size of the vectors u a and u n .To determine this error, some aberrant values have been removed from the calculation.For some points calculated close to the wall and when the direction of the fluid flow changes, the velocity value is sometimes very close to zero.The models give sometimes very large difference for these points, probably due to a lack of precision in the mesh or to approximate values.On average, about 4 % of the calculated points were not considered in the final MSRE calculation due to this extreme difference.This difference is more representative of the accuracy of the model because the points are comparated on a whole period and all over the channel surface as shown in for a given time value.A slider is used to move the plot over time and see how the shape of the fluid flow profile moves over time and then returns.

Results synthesis
All numerical simulations were performed with the so-called "fine" mesh of COMSOL Multiphysics® excepted for the simulations at 0.01 Hz and 0.1 Hz for a square channel where a "finer" mesh has been used.
In terms of computation time, the analytical model is solved faster than the numerical one.The computation time of the analytical model is approximatively 5.8 seconds.For each simulation, the discretization of the analytical model is 50 points for x and y axes.For z-axis, the discretization is 20 points (i.e, with 1 point for the discretization in z-axis, the model becomes a 2-D model).The number of temporal and spatial harmonics is 100.Except in the case of a sinusoidal pressure signal where the number of temporal harmonics is zero (i.e., this is why the computation time is slightly faster).The time required for postprocessing the results has not been considered to compare the models.However, for the numerical model, post-processing also requires a waiting time of a few minutes.For the analytical model, the results can be observed instantly.Another performance presented by the analytical model is the accuracy of the results.Comparisons between the two models show good agreement for any simulation with a MSRE always below 4 %.Simulations with an even finer mesh size would have been desirable, but due to time and technical limitations of the computers, the simulation was limited to "fine" and "finer" mesh sizes.

Conclusions and prospects
As presented in the introduction, extensive bibliographic research have not revealed any 3-D analytical fluidic models in cartesian coordinates studying AMRR cycles.Most models are only 1-D models or numerical, which often leads to significant inaccuracies in the results and longer computation times.
The resolution of the incompressible Navier-Stokes' PDE applied to a Cartesian geometry has been presented and detailed.The velocity evolution, directly in 3-D steady-state, is determined for any type of pressure gradient.The transient phase up to steady state is not interesting for studying AMRR cycles and is not taken into account by the model.This is an advantage in terms of computation time.The computation time and results of the two models were compared for one period.The analytical model has an extremely short computation time compared to the numerical method and also a good accuracy compared with the results of the COMSOL Multiphysics® software in "fine" and "finer" precision.The required time and the number of simulated periods until the steady-state are important factors to consider when evaluating the performance of a model intended for the study of AMRR cycles.
It can be observed that the channel dimensions and frequency have a significant impact on the evolution of the fluid flow profile.Further investigations need to be carried out to find optimal operating parameters to enable better heat exchanges for AMRR.A coupling of thermic and magnetic analytical models would allow a systemic study of AMRR cycles with a considerably reduced computation time, as demonstrated in this paper.

Figure 1 .
Figure 1.Schematic of a parallel plate AMRR device.Figure 2. Channel geometry considered by the analytical model.

Figure 2 .
Figure 1.Schematic of a parallel plate AMRR device.Figure 2. Channel geometry considered by the analytical model.

Figure 4 .
Figure 4. Annular effect observed for a square channel, for f = 2Hz and α = 17.71 with (a) x-y-view and (b) side view.

Figure 5 .
Figure 5.Comparison of analytical and numerical results over one period for one point of the space with f = 0.1 Hz for (a) a square channel geometry (ε r = 1.01%) and (b) a thin width channel (ε r = 3.55%).

Figure 6 .
Figure 6.Comparison of analytical and numerical results over the whole surface in the channel depth and one point in time.

Table 1 -
Comparison synthesis between both models Number of simulated periods until to reach steady state for the numerical model.b Computation time has been calculated only for one period.c MSRE is the Mean Square Relative Error, detailed further in the paper.