A numerical model of thermoacoustic heat pumping inside a compact cavity

– This paper presents a numerical study of thermoacoustic heat pumping along a stack of solid plates placed inside a compact cavity submitted to an oscillating ﬂ ow. Velocity and pressure ﬁ elds are controlled by two acoustic sources: a main “ pressure ” source monitoring the ﬂ uid compression and expansion phases, and a secondary “ velocity ” source generating the oscillating ﬂ uid motion. Numerical simulations are performed with an “ in-house ” code solving Navier – Stokes equations under a Low Mach number approximation in a two-dimensional geometry. In the linear regime, thermoacoustic heat pumping is correctly described with this model for different sets of parameters such as thermo-physical properties of the stack plates, amplitude of pressure oscillation or of the velocity source, phase shift between both sources. Numerical results on the normalized temperature difference established between the ends of stack plates are in excellent agreement with analytical estimates and experimental results published in the literature. Several con ﬁ gurations corresponding to different thermal conditions applied on the outside wall and an inside separation plate are then considered. If the separation plate is adiabatic, temperature varies linearly along the stack, recovering classical linear theory ’ s results. If the separation plate is thermally conductive, the model, providing detailed description of local heat and mass transfer, shows that the temperature ﬁ eld becomes fully two-dimensional and thermoa-coustic heat pumping is less ef ﬁ cient. The model is well adapted to explore the in ﬂ uence of local heat transfer constraints on the heat pump ef ﬁ ciency and thus well suited for detailed analyses of more complex mechanisms such as buoyancy effects.


Introduction
Thermoacoustic effect is a phenomenon of heat transport driven by acoustic oscillations. It occurs in the gas oscillating in the vicinity of a conducting solid, within viscous and thermal boundary layers. Thermoacoustic devices are used to generate either heat pumping for refrigerating purposes or mechanical work for engine operating purposes. The energy conversion takes place in a thermoacoustic core consisting of a stack of solid plates or a porous medium made of mesh grids or another geometry, placed between two heat exchangers. The thermoacoustic core is usually placed inside a closed resonator in which the gas is working at the fundamental resonant frequency associated with the geometry of the resonator [1][2][3].
Thermoacoustic engines and refrigerators are interesting for industrial applications because they are an ecological alternative of classical thermal machines, since they have very few moving parts, use inert gas and are simple to build.
They are however less competitive because of a much lower efficiency, which is limited by acoustics. For more than 30 years, different configurations have been studied for various types of resonators and thermoacoustic cores and coupling methods between engines and refrigerators. A wide variety of thermoacoustic devices were constructed in order to improve the overall performance for each application: cryogenics [4], air conditioning [5], domestic refrigeration, solar powered thermoacoustic engines [6], electric energy generators [7], waste heat recovery systems [8].
In the last 20 years efforts have focused on reducing the dimensions of the systems, especially for applications in cooling microelectronic components [9]. The present study deals with thermoacoustic refrigerating systems. When reducing the size of a device with a classical configuration, the resonance working frequency is significantly increased, and it is difficult to construct thermoacoustic cores well adapted to high frequencies, so that the resulting systems have limited performances [10]. A different design based on a co-axial travelling wave refrigerator [11] allowed to reduce the resonator's length and to obtain a more compact device with better performances than previously developed standing-wave refrigerators (COP equal to 25% of the Carnot COP, representing improvement by a factor of two). Smith et al. [12] developed a more compact travelling-wave system in which the resonator is replaced by a hybrid acoustic-mechanical system.
Another solution was to eliminate completely the resonator and create the acoustic field needed in the thermoacoustic core with two acoustic sources inside a non resonant compact cavity [13][14][15] (i.e. the cavity size is much smaller than the wavelength associated to the frequency of the acoustic sources). The main driver acts as a pressure source and the second source is a loud-speaker acting as a velocity source. The thermoacoustic core is a regenerator (mesh of stainless steel screens) or a stack of solid plates placed between heat exchangers. The overall geometry is coaxial with a peripheral channel allowing a feedback flow. Therefore the fluid can oscillate in this closed cavity in a similar way as in a travelling wave system. The acoustic sources are working at the same frequency and their amplitudes and phases are tuned independently. The device can be directly optimized by controlling the magnitude and phase shift of acoustic sources.
Following this idea a compact thermoacoustic refrigerator prototype for automobile air conditioning applications was designed and built recently [16]. The prototype was designed using the numerical code DeltaEC [17], based on the linear theory of thermoacoustics (Rott's linear theory [18,19], later extended by Swift [2], Arnott et al. [20]). The numerical implementation is based upon a lumped element approximation, each element being a part of the system: stack, heat exchanger, duct, etc. For each element, the geometry of the transverse section and the characteristics of the solid are accounted for through mean values over the element section, in the quasi-one dimensional set of equations, written in the frequency domain. The code provides an estimate of one-dimensional distributions of pressure, velocity and temperature, for a given system working in the periodic regime. Until now, this approach has been widely used to design thermoacoustic devices, despite observed discrepancies on the thermal quantities [16,21]. These could be associated with nonlinear phenomena such as streaming [22], end effects or natural convection [23,24], known to be responsible for reducing the efficiency of thermoacoustic systems. Those phenomena are partially taken into account or simply neglected in the context of the linear thermoacoustic theory. Previous CFD studies of thermoacoustic standing wave refrigerator have been based upon a fully viscous, conductive, and compressible solution [25] for high frequencies (20 KHz), simulating devices of extremely small size (on the order of millimeters). Although promising, this approach would lead to extremely large simulations for the present application (cavity size on the order of ten centimeters, and low frequency of 50 Hz). Simulations based on the linearized Navier-Stokes equations (LNSE) are commonly used to implement multidimensional acoustic simulations, in which a background flow interacts with and modifies an acoustic field [26]. This approach could be interesting but is restricted to the linear regime.
In the present work a numerical simulation of the flow and heat transfer in a compact simplified thermoacoustic heat pump is proposed. It consists in solving Navier-Stokes equations under a low Mach number approximation (using the numerical solver SUNFLUIDH [27]). The low Mach number approximation has already been used to describe a thermoacoustic stack refrigerator in the years 1998-1999 [28,29] and later by [30]. The simulation of the oscillation flow and the heat pumping effect was performed in a region restricted between two stack plates with boundary conditions that meant to represent the acoustic resonator. More recently, the same approach was used to simulate idealized thermoacoustic engines, by coupling a numerical model of nonlinear flow and heat exchange inside the stack with linear acoustics in the resonator and a variable load [31,32]. In the present study the low Mach number approach can be applied to describe finely the entire device since the frequency of the acoustic sources is such that the associated wavelength is much greater than the cavity size. As a consequence, pressure oscillations can be considered as spatially uniform in a first approximation. This model is capable of simulating non-linear hydrodynamics interactions with lower numerical cost than CFD simulations.
The main objective of this work is to show the relevance and interest of such an approach for the simulation of heat pumping in compact cavities. For the sake of simplicity a 2D geometry of the cavity is considered. The simulations provide a fine description of the flow fields and heat transfer inside the cavity allowing to account for their two-dimensional features and to analyse their effect on thermoacoustic heat pumping. Moreover, we point out that thermoacoustic effects can be correctly represented by coupling two mechanisms: a compression-relaxation cycle and an oscillating fluid motion. The former is related to a homogeneous thermodynamic pressure component whose time variation is governed by a first source acting as an "acoustic driver". The latter is sustained by a second source, the "oscillating fluid generator" that imposes a sinusoidal time-variation of the flow rate. This point of view is quite different from the linear theory of thermoacoustics where the main role of acoustic propagation and local pressure gradients is emphasized. Part of this work was presented during the CFA 2022 [33]. The paper is organized as follows: After a short presentation of the low Mach number approach coupled with specific boundary conditions modeling both acoustic sources, and of numerical methods, a reference simulation setup is described (geometry, thermophysical properties, . . .). In the result section, thermoacoustic heat pumping is first analyzed for different key parameters such as the stack thermal conductivity, the flow velocity magnitude, the drive ratio, the relative phase between pressure oscillations and particle velocity oscillations. Numerical results are compared with analytical estimates and experimental results presented by Poignand et al. [13,15]. Several configurations corresponding to different thermal conditions applied on the outside wall and the separation plate are then considered and their impact on thermoacoustic heat pumping is discussed. For a perfect gas, the flow dynamics is commonly described by the compressible, viscous, and heat-conductive Navier-Stokes equations, namely the conservation laws of mass, momentum, energy and the ideal state equation: where q stands for the density, V for the velocity field, f V for the volume forces (i.e. gravity force) and s for the stress tensor, which is defined as s ¼ lðrVþ T ðrVÞÀ 2 3 r Á V À Á IÞ, with l the dynamic viscosity. In the energy equation, T represents the temperature, c p the specific heat at constant pressure, j the thermal conductivity and U = s : rV the viscous dissipation rate. In the state equation, R is the universal gas constant (R = 8.314 J.K À1 .kg À1 ) and M the molar mass of the gas. The viscous dissipation rate is often negligible compared to the other terms for low speed flows and this was checked in the present study. It will therefore no longer be considered here.
Solving the fully compressible Navier-Stokes equations can be too time consuming since it requires a numerical time step very small in order to verify an acoustic Courant-Friedrichs-Lewy condition, based on the speed of sound as velocity scale.
As already mentioned in introduction, the dimensions of the prototype under study are much smaller than the wavelength associated with the oscillation frequency of the two sources, i.e. the Helmholtz number He defined as the ratio of a reference length scale L of the cavity to the wavelength k 0 is small. Acoustic wave propagation can therefore be neglected and the whole prototype considered as acoustically compact. Navier-Stokes equations can then be simplified using asymptotic expansions in powers of Mach number (M ¼ V 0 c 0 ) yielding the so-called low Mach number model that filters acoustic waves, following the approach proposed in [34,35]. As a result, the pressure P is split into two terms : the first order thermodynamic pressure component P th and a much smaller dynamic pressure component P dyn of order OðM 2 Þ. The first term P th is homogeneous in space, depends only on time and is related to the global time-variation of the thermodynamic state of the system during an expansion/relaxation cycle. The second term, which depends on both space and time, represents the small variations of pressure that directly acts on the flow dynamics. The Navier-Stokes equations are rewritten as: The low Mach number model allows to finely describe flow dynamics and heat transfer for compact thermoacoustic configurations, while being less time-consuming than a fully compressible model. In this study, the buoyancy effect is not considered (f V = 0). This point will be discussed in Section 4. Since the pressure is split into two components an added equation is required to close the system. The relation yielding P dyn is described in the following Section 2.2.
The relation for P th is obtained by integrating over the whole domain the state Equation (8), yielding: taking into account that the entire mass of fluid is known at any time in the domain X.
In the solid parts of the stack, the following heat equation is solved, denoting with the subscript s the thermophysical properties of the stack plates: Continuity of heat flux is applied at fluid-solid interfaces.
The main steps of the algorithm are outlined in the next section.

Numerical methods
The governing equations are solved following a finite difference approach of second order in time and space. Convective/advective fluxes and conductive/viscous terms are discretized with a second order scheme on a staggered grid. The time scheme is based on a Crank-Nicholson method where the diffusive terms are implicit in order to increase the numerical stability of the method with respect to the time step, which is based on the CFL criterion only. The dynamic pressure component is solved by means of an incremental projection method commonly used for the simulation of incompressible flows [36] and adapted for the low-Mach number approach. For that, the divergence-free constraint is replaced by the mass conservation constraint by using Equation (5). The projection method also allows to update the velocity field estimated from the momentum Equation (6) in such a way that mass conservation is satisfied. More information can be found in [37], for example. To sum up, this leads to a Poisson's equation of the form: where W ¼ P nþ1 dyn À P n dyn is the time increment of pressure, the index n denotes the time iteration, V * is the velocity field estimated from the momentum Equation (6), q the mass density obtained from the state Equation (8), oq ot its time derivative and Dt the numerical time step. In order to satisfy the mass conservation Equation (5), the velocity field is updated as: The main steps of the algorithm used for the time advancement of physical quantities can be summarized as: 1. Solve Equation (10) to update the temperature field. (9) and (8) to obtain the thermodynamic pressure component P th and the density q respectively/in that order. The entire mass of fluid in the domain is known at each time step from initial and boundary conditions (see Sect. 2.3.2 below). 3. Solve Equation (6) to estimate the velocity field V * (which does not satisfy mass conservation). 4. Use the projection method to compute the dynamic pressure component P dyn and correct the velocity field in order to satisfy mass conservation.

Wall boundary conditions
Wall boundary conditions considered in this numerical study are classical. The usual impermeability and no-slip conditions (V = 0 on the solid walls) are used for the velocity. Since velocity values are fixed, the projection method naturally imposes a zero normal derivative for the dynamic pressure component P dyn . Several cases are considered about heat transfer. Either walls are adiabatic, which implies to impose a zero normal derivative for temperature, or wall temperature is imposed and set to a constant value T w . For thermally conductive solids, continuity of temperature and heat flux at the solid-fluid interface is ensured in the model, since the enthalpy equation naturally reduces to the heat equation within the solid.

Modeling the acoustic sources
The experimental technology used to monitor the governing pressure variation and oscillating fluid motion consists in various acoustic sources, like pistons and loudspeaker membranes. These complex elements are difficult to model in a CFD code. This would require specific methods managing moving objects in a fluid, like ALE or immersed body methods, as well as a finer mesh refinement in the zones of interest. As a consequence, numerical simulations would be particularly complex and time consuming. To overcome this difficulty, we have opted to model the acoustic sources by using oscillating flowrate inlets (see Fig. 1). For each acoustic source, the boundary conditions are defined as: where Q m (i) is the mass flow rate at each inlet i with an amplitude q 0 U i S i (reference values of fluid density, velocity and surface of the inlet i) and a time variation defined from the oscillation frequency of the device f 0 and the phase / i . This formulation ensures mass conservation over each thermoacoustic cycle and allows us to know the mass of fluid at any time. For others quantities (temperature, density, dynamic pressure component) a zero normal derivative is imposed. For sake of clarity the inlets are reported in red color in the sketch of the computational domain (see Fig. 1). The oscillating fluid generator (OFG) is constructed as two back-to-back inlets which have the same parameters q 0 , U OFG , S OFG , / OFG . The OFG model therefore produces a zero mass balance at any time and plays a role in the fluid motion only. The acoustic driver (AD), which governs the the compression/relaxation cycle through the variation of pressure, is modeled as a single inlet whose parameters are q 0 , U AD , S AD , / AD . By considering the phase of the acoustic driver as a reference (/ AD = 0) and the quadratic phase delay between the time variation of pressure and the acoustic driver velocity, the usual phase shift between velocity and pressure can be expressed as where m is an integer). The drive ratio (i.e. the ratio between the amplitude of pressure oscillation and the mean pressure) and the velocity scale in the stack are directly related to the flow rates of the acoustic driver and the oscillating fluid generator.  In this section, we describe the 2D configuration used in our numerical simulations (in the context of this study, we consider that any phenomenon associated with the third direction can be neglected). It includes the main components of compact thermoacoustic refrigerators designed in laboratory experiments [14,16]: an acoustic driver (AD) controling the compression and relaxation cycle, an oscillating fluid generator (OFG), a peripheral channel (PC) allowing circulating fluid and a stack of 16 solid plates along which is created the temperature gradient (see Fig. 1). The geometry is cartesian and the computational domain is restricted to the upper part of the device by assuming symmetry of the flow behavior with respect to the median-plane. The main geometrical features of the device are summarized in Table 1.
The physical properties of the stack plates are chosen to be representative of materials commonly used in thermoacoustic devices. The mass density is fixed at q s = 2000 kg.m À3 . and the thermal conductivity is set to j s = 0.2 W.m À1 .K À1 (except in Section 4.1.1 where different values of j s are considered). The specific heat capacity is however set to an artificially small value (c ps = 10 J.K À1 .kg À1 ). The objective is to shorten the transition time before reaching a converged time-average state of flow and heat transfer, and thus reduce the CPU time. The value of c ps has a negligible impact on the final temperature gradient in the stack as will also be shown in Section 4.1.1.
The fluid is a gas mixture composed of 30% argon and 70% helium as in the TACOT experiments [16]. This choice is made in order to verify the ability of the numerical model to supply a high temperature gradient of the same order of magnitude as those obtained with experimental prototypes [16]. The perfect gas state law is considered with constant physical properties. The specific heat capacity, thermal conductivity and dynamic viscosity are set to c p = 1404.97 J.K À1 .kg À1 , j = 8.56Á10 À2 W.m À1 .K À1 and l = 2.11Á10 À5 kg.m À1 .s À1 respectively.
The frequency of both the acoustic driver and the oscillating fluid generator is fixed at f 0 = 50 Hz. By considering a reference temperature of T 0 = 298 K and mass density of q 0 = 0.597 kg.m À3 (corresponding pressure P 0 = 10 5 Pa), the thicknesses of the dynamic and thermal boundary layers along the stack plates are estimated to be d m = 4.75Á10 À4 m and d j = 8.07Á10 À4 m respectively.
The value of the resulting temperature difference along the stack plate and therefore of the heat pump efficiency is potentially dependent upon constraints on local heat transfer within the cavity. Several configurations corresponding to different thermal boundary conditions applied on the outside wall and the separation plate are thus considered: Case 1: All solid walls are adiabatic (except for the stack plates). This configuration corresponds to an ideal situation with no heat transfer through external walls or the separation plate between the stack region and the peripheral channel (locations are shown on the sketch of the device, Fig. 1). Case 2: All external walls are at fixed temperature T w and the separation plate between the stack region and the peripheral channel is adiabatic. It is also a simplified configuration corresponding to isothermal outer walls, but it allows for heat loss to the ambient surroundings. In real situations there is also convective heat transfer with the surroundings and a fine determination of the overall heat loss would require complex numerical simulations which is beyond the scope of this study. The wall temperature T w is arbitrarily set to 285 K. Case 3: All external walls are at fixed temperature T w = 285 K and the separation plate between the stack region and the peripheral channel is conductive. For the sake of simplicity, the values of heat capacity and mass density are the same as for the stack. The thermal conductivity is set to either j 1 = 0.25 W.m À1 .K À1 (case 3a : moderate conductivity specific to ceramics [15]) or j 2 = 14 W.m À1 .K À1 (case 3b: large conductivity specific to stainless steel [16]) accounting for the variability of materials used in thermoacoustic devices.

Numerical parameters
The computational domain is discretized according to a cartesian grid 384 Â 320 in the longitudinal (x) and the normal (y) directions respectively (with respect to the stack plates orientation). In the normal direction, the cell size distribution is regular over the stack area. There are 16 cells across the plate thickness as well as across the fluid interspace between plates, the resulting cell size being Dy = 8Á10 À5 m. This ensures a correct spatial resolution of the boundary layers.
In the longitudinal direction, the cell size varies following a specific distribution law so that the minimum values are positioned nearby the stack ends (Dx min = 2.5Á10 À4 m).
The time step is fixed at Dt = 2Á10 À5 s. This choice ensures that the CFL criterion is satisfied (maximum value is about 0.5) and each period of the thermoacoustic cycle is exactly described with 1000 time steps. Initially the fluid is at rest. Initial conditions for pressure and temperature are set to P 0 = 10 5 Pa and T 0 = 298 K and the mass density is equal to q 0 = 0.597 kg.m À3 .

Numerical results
The numerical results presented below are obtained for physical configurations adapted from experiments [14,16] so that the operating mode mainly corresponds to the linear regime. The working gas and frequency were chosen to be the same as in the TACOT prototype [16] in order to obtain a thermoacoustic heat pumping effect large enough for industrial applications. However, the mean thermodynamic pressure was set to atmospheric pressure in order to keep boundary layers large enough to be discretized without excessive numerical cost. Here the ratio between the halfdistance between plates and the thermal boundary layer thickness is equal to 0.8 which means that the stack is neither quasi-adiabatic nor quasi-isothermal, but in between.
Numerical simulations are then compared with analytical results obtained from linear theory [13] for various sets of parameters as well as with experimental measurements of Poignand et al. [15], with either a quasi-adiabatic stack (Exp. 1) or a quasi-isothermal regenerator (Exp. 2). Equations for the linear theory are solved with parameter values presented in Section 3.1, in order to describe the ideal case 1. The main input parameters characterizing the different simulations are the drive ratio, the stack velocity scale and the phase shift between the acoustic driver and the oscillating fluid generator. The drive ratio is the commonly used parameter to characterize the acoustic operating regime. It is defined with respect to the thermodynamic pressure component as: In the present study, the drive ratio is small enough (Dr < 2.5%) so that Rott's linear theory should give the correct value for the temperature difference between stack ends (DT). The stack velocity scale u s is commonly defined as the mean velocity amplitude over the stack channels. It can be simply estimated as u s ¼ U OFG R b (the density variation in space being neglected as a first approximation), where U OFG is the velocity amplitude of the OFG, and R b is the stack porosity (ratio between the fluid part and entire stack volumes).
In the following sections, results obtained on the three configurations presented in Section 3.1 are detailed. These results mainly focus on the mean value of DT computed from time series recorded at locations reported in Figure 1. This quantity is normalized in order to compare numerical results with experimental measurements whose operating regimes are not the same as those used for computations. Temperature fields are also presented to point out how thermoacoustic heat pumping can be influenced by the previously described thermal configurations. In order to smooth the time variations related to each thermoacoustic cycle, all physical variables are averaged over 100 oscillation cycles. Time series show the evolution of the averaged quantities. Simulations are run until convergence is obtained for the mean value of DT (or quasi-convergence, see discussion in Sect. 4.2).

Parametric analysis of thermoacoustic heat pumping
In this section, the sensitivity of the thermoacoustic heat pumping effect, given by DT, with respect to some input parameters (driver ratio, stack-velocity scale, . . .) is shown. The main objective is to validate the model in the linear regime by comparison with experimental measurements and analytical results obtained from the linear theory. Results corresponding to the configuration case 1 only were already described in CFA 2022 [33]. In the following discussion the reference operating conditions correspond to a drive ratio equal to Dr = 0.0236, a velocity scale in the stack on the order of u s = 0.5 m.s À1 (U OFG = 0.25 m.s À1 and R b = 0.5), and the mass flow rates associated with the AD and OFG sources are adjusted in such a way that the phase shift between the velocity and the pressure in the stack is equal to U u À U p = 5p/6, as explained in the end of Section 2.3.2. These reference values were chosen in order to obtain a significant thermoacoustic effect for all configuration cases. All other input parameters are described in Section 3. The sensitivity to each parameter is analysed separately.

Influence of stack thermo-physical properties
Results shown in this paragraph were obtained for the configuration case 1. The considered thermo-physical properties are the heat capacity q s c ps and the thermal conductivity j, since the heat flux is essentially determined by these parameters. The time evolution of the temperature difference DT is plotted in Figure 2 for several values of the stack plates specific heat capacity c ps (the actual parameter should be q s c ps but it is obvious that considering only the variation of c ps is sufficient, q s being fixed). There is a transient phase before reaching an asymptotic regime giving DT ' 51.5 K, which, as expected, is longer as the value of the heat capacity is larger (see Fig. 2). The low sensitivity of DT in the asymptotic regime to the heat capacity of the stack allows us to consider a small value in order to reduce drastically the computation time and justifies our choice to set c ps = 10 J.K À1 .kg À1 . Similar observations were made for the other configurations (cases 2 and 3). Therefore the same value of c ps was kept for all simulations. Figure 3 shows the value of DT between stack ends for six values of the thermal conductivity j s of the solid stack plates (case 1). The solid line shows the corresponding prediction of the linear theory of thermoacoustics. The agreement is excellent. As expected, when j s is very large, DT becomes negligible because of enhanced heat conduction in the solid: thermoacoustic heat pumping becomes ineffective. For j s = 0 (adiabatic stack) the thermoacoustic effect is canceled, as predicted by Swift [1]. Note that for very small (non-zero) values of j s the numerical prediction departs for the linear theory, suggesting that the evolution as j s ? 0 is sharp but continuous and that there is an optimum j s value for which DT is maximum. Again similar observations were made for the other configurations (cases 2 and 3).

Influence of the Drive ratio
In this paragraph, the used relative velocities and phase shift are fixed at their reference values u s = 0.5 m.s À1 and / u À / p = 5p/6. Figure 4 shows the variation of the normalized temperature difference along the stack as a function of the drive ratio for all considered configurations (cases 1 to 3). Numerical points are plotted as well as the linear theory prediction and experimental results from [15]. In order to compare our results with those of Poignand et al. [15], which are performed in very different physical conditions, the reference scale DT ref used to normalize DT is chosen at Dr = 0.01. Both sets of experimental points (Exp. 1) and (Exp. 2) are shown on the figure, and the evolutions are consistent with the numerical simulation results which fall in between as expected, since in the numerical configuration the stack is neither fully adiabatic nor isothermal. Very good agreement with the linear theory is observed (see Fig. 4). We notice the numerical results are very similar independently from the considered configuration.

Influence of the stack velocity-scale
Another important parameter of the thermoacoustic heat pumping effect is the characteristic velocity scale within the stack. The variation of this parameter is obtained by changing the input velocity U OFG of the oscillating fluid generator. The other relevant parameters, the drive ratio and the phase shift, are fixed at their respective reference value (Dr = 0.0236 and / u À / p = 5p/6). In order to compare numerical, theoretical and experimental results, the stack velocity scale corresponding to each data set is normalized with respect to its own optimum value u opt providing the maximum value of the DT in each case. Figure 5 shows the evolution of DT between stack ends as a function of the normalized stack velocity scale u Ã s ¼ u s =u opt for all considered configurations. Again a similar behavior is observed in numerical or experimental results and in the linear theory. The value of DT resulting from the different approaches increases with u Ã s , up to u Ã s about 0.85, then reaches the same plateau. There are however some quantitative discrepancies for u Ã s < 0.6. While the experiment 1 fits well with the linear theory, the experiment 2 yields lower values of DT. The numerical results present larger values which are similar for the three considered cases. The scattering of experiment 2 and numerical results around the linear theory curve is of the same order and is more noticeable as u Ã s decreases. However the normalization allows to obtain analogous behavior for all these physically very different configurations.   In this paragraph, the used relative velocities and drive ratio are fixed at their reference values (u s = 0.5 m.s À1 , Dr = 0.0236). The influence of the phase shift between velocity and pressure is studied and results are plotted in Figure 6 for the three numerical configurations. Here again a similar evolution is observed in numerical simulations, experimental results and the linear theory. In the [Àp, 0] interval, numerical results are bounded by linear theory and experimental results, while in the [0, p] interval all curves are in close agreement. Numerical simulations, for every cases, confirm the experiment results in that the phase shift has a major influence on the thermoacoustic heat pumping effect. Indeed it is possible to control the orientation of heat pumping (and invert DT between the stack ends) with the phase shift. The observed asymmetry of the curve (for example on extreme values) is probably correlated to the device geometry with respect to the stack middle cross-section. Note that the model deduced from the linear theory does not predict this asymmetry. It is observed that for all considered configurations, the numerical value of the phase shift U leading to the maximum value of DT varies between 5p/6 and p which is slightly larger than the value found in the linear model (3p/4) and close to the experimental results reported in [15].

Influence of heat transfer constraints (choice of configuration)
In addition to the very good agreement between the numerical results, experimental measurements and linear theory, we found that the relative variation of the thermoacoustic heat pumping effect with respect to parameters defining the operating regime (drive ratio, stack velocity scale or phase shift) does not depend significantly on the considered configuration (i.e. on the heat transfer constraints). In other words, the normalized temperature difference across the stack plates follows the same scaling laws.
We now analyse how the heat loss through external walls and the thermal conductivity of the separation plate act on the magnitude of the thermoacoustic heat pumping effect. We use the reference operating conditions described in Section 4.1 for the three considered configurations. Figure 7 describes the time evolution of DT for these different cases.
In both configurations for which the separation plate is adiabatic (cases 1 and 2), the established values of DT are close to each other. Moreover, the mean temperature fields for these two cases present very similar spatial distributions over the computational domain (see Figs. 8a and 8b). The temperature distribution is unidirectional and quasi linear along the plates in the stack region with same gradient value while it is rather uniform in other parts of the device. From this observation we can deduce that the established time-average heat flux transported by the flow along the peripheral channel is relatively small and the adiabatic separation plate is sufficient to ensure an efficient thermal insulation of the stack region, thereby strongly reducing   the impact of the heat loss on the thermoacoustic heat pumping efficiency. Only the global temperature level is larger in the adiabatic external walls configuration (since no heat loss is permitted) than that obtained for fixed temperature outer walls, with heat losses in the latter configuration occurring mainly during the transient stage.
Note that the adiabatic external wall configuration (case 1) induces a small but continuous variation of DT over a long time range. As a consequence, perfectly converged state of the system (in the sense of time average quantities) cannot be obtained especially for operating regime close to optimal conditions. This can be explained by the fact that the adiabatic conditions coupled with the mechanism of expansion-compression of the gas cannot maintain the same thermodynamic state of the fluid from one thermoacoustic cycle to another. In other words, ensuring time-average equilibrium between the enthalpy variation of the fluid over a thermoacoustic cycle and the heat exchange with the stack over any thermoacoustic cycle is very difficult. This small shortcoming however does not affect our results if the simulation takes place over a limited time (for example, 400 s in configuration case 1, see Fig. 7).
What about the absolute magnitude of the temperature difference established between stack ends? The established value of DT obtained numerically in case 2 is DT ' 56 K. Rott's linear theory developed in [13] yields DT lin = 49 K with the numerical input data set (case 2), in accordance with our numerical result.
We now examine the influence of the thermal conductivity of the separation plate, by comparing results obtained in cases 2, 3a and 3b. As expected, the larger the thermal conductivity of the plate, the smaller the thermoacoustic heat pumping efficiency, as shown in Figure 7. DT is reduced by a factor of about 2-2.5 between the adiabatic case (case 2) and both thermally conductive cases (3a and 3b). Moreover, the time required to reach a converged state is drastically shorter in cases 3a and 3b. When the separation plate is thermally conductive, a significant heat flux contribution appears in the normal direction which enables heat loss through the peripheral channel and the external walls and reduces the magnitude of the mean temperature gradient along the stack plates. As a consequence, the temperature distribution over the computational domain is drastically modified. When the separation plate is adiabatic, it is organized following a linear distribution along the plates inside the stack plates and a quasi-uniform distribution elsewhere (as seen in Fig. 8b). When the separation plate is thermally conductive, the temperature field presents a 2D pattern resulting from heat flux contributions along each direction (see Fig. 8c).
If we consider devices with larger span, the vertical distribution of the temperature field in the stack region could trigger buoyancy effects and in turn strongly modify the thermoacoustic heat pumping effect along the stack. As a consequence, thermal insulation of the stack region appears to be a key point in the device optimization, minimizing the vertical mean temperature gradient.
These observations constitute a significant added contribution of numerical simulations compared with linear theory and experiments.

Conclusion
A numerical model based on a Low Mach number approximation was designed to simulate the thermoacoustic heat pumping effect inside a compact cavity consisting of an acoustic driver, an oscillating fluid generator and a stack of solid plates. In this fluid mechanics approach, the physical coupling between the oscillating fluid and homogeneous compression-relaxation cycles results in thermoacoustic heat pumping along the stack plates with no account of acoustic propagation. The main characteristics of thermoacoustics are recovered: the normalized temperature difference as a function of a wide range of operating conditions is in very good agreement with the linear theory predictions and experimental measurements, at least in the linear regime.
The model was used to explore a wide range of input parameters in the linear regime for several numerical configurations, and results were compared with the linear theory as well as experimental measurements. The main results are: Firstly, the normalized temperature difference between stack plate ends follows a scaling law with respect to the main parameters used to determine the operating conditions (drive ratio, stack velocity scale or phase shift) independently from the configuration.
Secondly, the influence of heat loss through external walls and of insulation characteristics of the stack region (through the separation plate) is studied. It is shown that the latter has a prevailing effect on the resulting thermoacoustic heat pumping efficiency, since the temperature distribution inside the stack is no longer unidirectional when the thermal conduction properties of the separation plate become significant. This leads naturally to investigate the relevance of gravity effects when the temperature field inside the stack region becomes multidimensional, especially for compact cavities with larger spans in which buoyancy could arise. The present model is very well adapted for such a study which is the object of our work in progress. In addition the model is also well adapted to explore other nonlinear effects, which can be associated for example with jet flow instabilities around stack ends.

Conflict of interest
Author declared no conflict of interests.