Issue 
Acta Acust.
Volume 7, 2023
Topical Issue  CFA 2022



Article Number  15  
Number of page(s)  11  
DOI  https://doi.org/10.1051/aacus/2023008  
Published online  12 May 2023 
Scientific Article
A numerical model of thermoacoustic heat pumping inside a compact cavity
^{1}
Université ParisSaclay, CNRS, LISN, 91400 Orsay, France
^{2}
Sorbonne Université, CNRS, Institut Jean le Rond d’Alembert, UMR 7190, 4 Place Jussieu, 75005 Paris, France
^{*} Corresponding author: dianageorgiana.balteancarles@sorbonneuniversite.fr
Received:
24
October
2022
Accepted:
24
March
2023
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 flow. Velocity and pressure fields are controlled by two acoustic sources: a main “pressure” source monitoring the fluid compression and expansion phases, and a secondary “velocity” source generating the oscillating fluid motion. Numerical simulations are performed with an “inhouse” code solving Navier–Stokes equations under a Low Mach number approximation in a twodimensional geometry. In the linear regime, thermoacoustic heat pumping is correctly described with this model for different sets of parameters such as thermophysical 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 configurations 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 field becomes fully twodimensional and thermoacoustic heat pumping is less efficient. The model is well adapted to explore the influence of local heat transfer constraints on the heat pump efficiency and thus well suited for detailed analyses of more complex mechanisms such as buoyancy effects.
Key words: Thermoacoustics / Low Mach number / Compact cavity / Numerical simulation / Thermoacoustic heat pump
© The Author(s), Published by EDP Sciences, 2023
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 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–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 coaxial travelling wave refrigerator [11] allowed to reduce the resonator’s length and to obtain a more compact device with better performances than previously developed standingwave refrigerators (COP equal to 25% of the Carnot COP, representing improvement by a factor of two). Smith et al. [12] developed a more compact travellingwave system in which the resonator is replaced by a hybrid acousticmechanical 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–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 loudspeaker 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 quasione dimensional set of equations, written in the frequency domain. The code provides an estimate of onedimensional 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 nonlinear 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 twodimensional 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 compressionrelaxation 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 timevariation 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.
2 Numerical models
2.1 Governing equations in the low Mach number approach
For a perfect gas, the flow dynamics is commonly described by the compressible, viscous, and heatconductive Navier–Stokes equations, namely the conservation laws of mass, momentum, energy and the ideal state equation:
(4)where ρ stands for the density, V for the velocity field, f_{V} for the volume forces (i.e. gravity force) and τ for the stress tensor, which is defined as , with μ the dynamic viscosity. In the energy equation, T represents the temperature, c_{p} the specific heat at constant pressure, κ the thermal conductivity and Φ = τ:∇V the viscous dissipation rate. In the state equation, R is the universal gas constant (R = 8.314 J.K^{−1}.kg^{−1}) and 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 λ_{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 () yielding the socalled 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 . The first term P_{th} is homogeneous in space, depends only on time and is related to the global timevariation 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 timeconsuming 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:
(9)taking into account that the entire mass of fluid is known at any time in the domain Ω.
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 fluidsolid interfaces. The main steps of the algorithm are outlined in the next section.
2.2 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 lowMach number approach. For that, the divergencefree 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:
(11)where is the time increment of pressure, the index n denotes the time iteration, V^{*} is the velocity field estimated from the momentum Equation (6), ρ the mass density obtained from the state Equation (8), its time derivative and Δt 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:
Solve Equation (10) to update the temperature field.
Use the state equation in both forms (9) and (8) to obtain the thermodynamic pressure component P_{th} and the density ρ 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).
Solve Equation (6) to estimate the velocity field V^{*} (which does not satisfy mass conservation).
Use the projection method to compute the dynamic pressure component P_{dyn} and correct the velocity field in order to satisfy mass conservation.
2.3 Boundary conditions
2.3.1 Wall boundary conditions
Wall boundary conditions considered in this numerical study are classical. The usual impermeability and noslip 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.
2.3.2 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:
(13)where Q_{m}(i) is the mass flow rate at each inlet i with an amplitude ρ_{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.
Figure 1 Sketch of the 2D TACOT configuration. The computational domain is restricted to the upper part above the symmetry plane. Circle tags show probe locations used for recording time series of temperature. 
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 backtoback inlets which have the same parameters ρ_{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 ρ_{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.
3 Numerical configurations
3.1 Description of the thermoacoustic device
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 medianplane. The main geometrical features of the device are summarized in Table 1.
Geometrical features of the thermocoustic device defined in the top halfdomain.
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 ρ_{s} = 2000 kg.m^{−3}. and the thermal conductivity is set to κ_{s} = 0.2 W.m^{−1}.K^{−1} (except in Section 4.1.1 where different values of κ_{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 timeaverage 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}, κ = 8.56·10^{−2} W.m^{−1}.K^{−1} and μ = 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 ρ_{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 δ_{ν} = 4.75·10^{−4} m and δ_{κ} = 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 κ_{1} = 0.25 W.m^{−1}.K^{−1} (case 3a : moderate conductivity specific to ceramics [15]) or κ_{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.
3.2 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 Δy = 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 (Δx_{min} = 2.5·10^{−4} m).
The time step is fixed at Δt = 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 ρ_{0} = 0.597 kg.m^{−3}.
4 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 quasiadiabatic nor quasiisothermal, 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 quasiadiabatic stack (Exp. 1) or a quasiisothermal 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 (ΔT). The stack velocity scale u_{s} is commonly defined as the mean velocity amplitude over the stack channels. It can be simply estimated as (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 ΔT 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 ΔT (or quasiconvergence, see discussion in Sect. 4.2).
4.1 Parametric analysis of thermoacoustic heat pumping
In this section, the sensitivity of the thermoacoustic heat pumping effect, given by ΔT, with respect to some input parameters (driver ratio, stackvelocity 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} − Φ_{p} = 5π/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.
4.1.1 Influence of stack thermophysical properties
Results shown in this paragraph were obtained for the configuration case 1. The considered thermophysical properties are the heat capacity ρ_{s}c_{ps} and the thermal conductivity κ, since the heat flux is essentially determined by these parameters. The time evolution of the temperature difference ΔT is plotted in Figure 2 for several values of the stack plates specific heat capacity c_{ps} (the actual parameter should be ρ_{s}c_{ps} but it is obvious that considering only the variation of c_{ps} is sufficient, ρ_{s} being fixed). There is a transient phase before reaching an asymptotic regime giving ΔT ≃ 51.5 K, which, as expected, is longer as the value of the heat capacity is larger (see Fig. 2). The low sensitivity of ΔT 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 2 Mean temperature difference ΔT between stack ends computed by numerical simulation for several values of c_{ps}, configuration case 1. 
Figure 3 shows the value of ΔT between stack ends for six values of the thermal conductivity κ_{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 κ_{s} is very large, ΔT becomes negligible because of enhanced heat conduction in the solid: thermoacoustic heat pumping becomes ineffective. For κ_{s} = 0 (adiabatic stack) the thermoacoustic effect is canceled, as predicted by Swift [1]. Note that for very small (nonzero) values of κ_{s} the numerical prediction departs for the linear theory, suggesting that the evolution as κ_{s} → 0 is sharp but continuous and that there is an optimum κ_{s} value for which ΔT is maximum. Again similar observations were made for the other configurations (cases 2 and 3).
Figure 3 Mean temperature difference ΔT between stack ends for 6 values of κ_{s}. Numerical results case 1 (blue dots), linear theory (solid line). 
4.1.2 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} = 5π/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 ΔT_{ref} used to normalize ΔT 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.
Figure 4 Normalized mean temperature difference ΔT/ΔT_{ref} between stack ends as a function of the Drive ratio Dr. Numerical results (cases 1, 2, 3), linear theory (black solid line), Exp. 1 (green crosses) and Exp. 2 (red triangles) from [15]. 
4.1.3 Influence of the stack velocityscale
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} = 5π/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 ΔT in each case.
Figure 5 shows the evolution of ΔT between stack ends as a function of the normalized stack velocity scale for all considered configurations. Again a similar behavior is observed in numerical or experimental results and in the linear theory. The value of ΔT resulting from the different approaches increases with , up to about 0.85, then reaches the same plateau. There are however some quantitative discrepancies for < 0.6. While the experiment 1 fits well with the linear theory, the experiment 2 yields lower values of ΔT. 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 decreases. However the normalization allows to obtain analogous behavior for all these physically very different configurations.
Figure 5 Normalized mean temperature difference ΔT/ΔT_{max} between stack ends as a function of the normalized stack velocity scale u_{s}/u_{opt}. 
4.1.4 Influence of the phase shift between velocity and pressure.
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 [−π, 0] interval, numerical results are bounded by linear theory and experimental results, while in the [0, π] 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 ΔT 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 crosssection. 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 Φ leading to the maximum value of ΔT varies between 5π/6 and π which is slightly larger than the value found in the linear model (3π/4) and close to the experimental results reported in [15].
Figure 6 Normalized mean temperature difference ΔT/ΔT_{max} as a function of Φ_{u} − Φ_{p}. 
4.2 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 ΔT for these different cases.
Figure 7 Time evolution of the temperature difference between stack ends for the different configurations. 
In both configurations for which the separation plate is adiabatic (cases 1 and 2), the established values of ΔT 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 timeaverage 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.
Figure 8 Mean temperature field over the computational domain for cases 1, 2 and 3a. 
Note that the adiabatic external wall configuration (case 1) induces a small but continuous variation of ΔT 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 timeaverage 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 ΔT obtained numerically in case 2 is ΔT ≃ 56 K. Rott’s linear theory developed in [13] yields ΔT_{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. ΔT 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 quasiuniform 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.
5 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 compressionrelaxation 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.
Acknowledgments
This research was financially supported by Agence Nationale de la Recherche (Projet ANR17CE060 0 0701).
References
 G.W. Swift: Thermoacoustics: A unifying perspective for some engines and refrigerators. Journal of the Acoustical Society of America, Melville, New York, 2002. [Google Scholar]
 G.W. Swift: Thermoacoustic engines. Journal of the Acoustical Society of America 84 (1988) 1145–1180. [Google Scholar]
 T. Biwa: Introduction to thermoacoustic devices. World Scientific Publishing Company, 2021. [CrossRef] [Google Scholar]
 W. Dai, E. Luo, J. Hu, H. Ling: A heat driven thermoacoustic cooler capable of reaching liquid nitrogen temperature. Applied Physics Letters 86 (2005) 224103. [CrossRef] [Google Scholar]
 L. Zoontjens, C. Howard, A. Zander, B. Cazzolato: Feasibility study of an automotive thermo acoustics refrigerator, in Proceedings of Acoustics, 9–11 November 2005, Busselton, Western Australia. 2005. [Google Scholar]
 J.A. Adeff, T.J. Hofler: Design and construction of a solarpowered, thermoacoustically driven, thermoacoustic refrigerator. Journal of the Acoustical Society of America 107, 6 (2005) L37–L42. [Google Scholar]
 S. Backhaus, E. Tward, M. Petach: Traveling wave thermoacoustic electric generator. Applied Physics Letters 85 (2004) 1085. [CrossRef] [Google Scholar]
 S. Spoelstra, M.E.H. Tijani: Thermoacoustic heat pumps for energy savings, in seminar “Boundary crossing acoustics” of the Acoustical Society of the Netherlands, 23 November 2005, ECN. 2005. [Google Scholar]
 O. Symko, E. AbdelRahman, Y. Kwon, M. Emmi, R. Behunin: Design and development of highfrequency thermoacoustic engines for thermal management in microelectronics. Microelectronics Journal 35 (2004) 185–191. [Google Scholar]
 B. Lihoreau, P. Lotton, G. Penelet, M. Bruneau: Thermoacoustic, small cavity excitation to achieve optimal performance. Acta Acustica United with Acustica 97 (2011) 926–932. [CrossRef] [Google Scholar]
 M.E.H. Tijani, S. Spoelstra: Study of a coaxial thermoacousticstirling cooler. Cryogenics 48 (2008) 77–82. [CrossRef] [Google Scholar]
 R. Smith, M. Poese, S. Garett, R. Wakeland: Thermoacoustic device. US Patent no. 0192324 A1, 2003. [Google Scholar]
 G. Poignand, B. Lihoreau, P. Lotton, E. Gaviot, M. Bruneau, V. Gusev: Optimal acoustic fields in compact thermoacoustic refrigerators. Applied Acoustics 68 (2007) 642–659. [CrossRef] [Google Scholar]
 G. Poignand, P. Lotton, G. Penelet, M. Bruneau: Small cavity excitation to achieve optimal performance. Acta Acustica United with Acustica 97 (2011) 926–932. [CrossRef] [Google Scholar]
 G. Poignand, A. Podkovskiy, G. Penelet, P. Lotton, M. Bruneau: Analysis of a coaxial, compact thermoacoustic heat pump. Acta Acustica United with Acustica 99 (2013) 898–904. [CrossRef] [Google Scholar]
 I.A. Ramadan, H. Bailliet, G. Poignand, D. Gardner: Design, manufacturing and testing of a compact thermoacoustic refrigerator. Applied Thermal Engineering 189 (2021) 116705. [CrossRef] [Google Scholar]
 B. Ward, J. Clark, G.W. Swift: Design Environment for Lowamplitude Thermoacoustic Energy Conversion (DeltaEC Version 6.2). Users Guide, Los Alamos National Laboratory, 2008. http://www.lanl.gov/thermoacoustics/UsersGuide.pdf. [Google Scholar]
 N. Rott: Damped and thermally driven acoustic oscillations in wide and narrow tubes. Zeitschrift für Angewandte Mathematik und Physik 20 (1969) 230–243. [Google Scholar]
 N. Rott: Thermoacoustics. Advances in Applied Mechanics 20 (1980) 135–243. [CrossRef] [Google Scholar]
 W.P. Arnott, H.E. Bass, R. Raspet: General formulation of thermoacoustics for stacks having arbitrarily shaped pore cross section. Journal of the Acoustical Society of America 90 (1991) 3228–3237. [Google Scholar]
 E. Matthew, E. Poese, S.L. Garrett: Performance measurements on a thermoacoustic refrigerator driven at high amplitudes. Journal of the Acoustical Society of America 107, 5 (2000) 2480–2486. [Google Scholar]
 C. Scalo, S.K. Lele, L. Hesselink: Linear and nonlinear modelling of a theoretical travellingwave thermoacoustic heat engine. Journal of Fluid Mechanics 766 (2015) 368–404. [Google Scholar]
 C. Shen, Y. He, Y. Li, H. Ke, D. Zhang, Y. Liu: Performance of solar powered thermoacoustic engine at different tilted angles. Applied Thermal Engineering 29, 13 (2000) 2745–2756. [Google Scholar]
 N. Pan, S. Wang, C. Shen: Visualization investigation of the flow and heat transfer in thermoacoustic engine driven by loudspeaker. International Journal of Heat and Mass Transfer 55 (2012) 7737–7746. [CrossRef] [Google Scholar]
 D. Marx, P. BlancBenon: Computation of the temperature distortion in the stack of a standingwave thermoacoustic refrigerator. Journal of the Acoustical Society of America 118 (2005) 2993–2999. [Google Scholar]
 A. Kierkegaard, S. Boij, G. Efraimsson: A frequency domain linearized NavierStokes equations approach to acoustic propagation in flow ducts with sharp edges. Journal of the Acoustical Society of America 127 (2010) 710–719. [Google Scholar]
 Y. Fraigneau: SUNFLUIDH : A software for computional fluid dynamics, User guide. LIMSI, 2013. https://sunfluidh.lisn.upsaclay.fr/, 2013. [Google Scholar]
 A. Worlikar, O. Knio: Numerical study of oscillatory flow and heat transfer in a loaded thermoacoustic stack. Numerical Heat Transfer, Part A 35 (1999) 49–65. [Google Scholar]
 A.S. Worlikar, O.M. Knio, R. Klein: Numerical simulation of a thermoacoustic refrigerator. II: Stratified flow around the stack. Journal of Computational Fluids 144 (1998) 299–324. [Google Scholar]
 P. Duthil, C. Weisman, E. Bretagne, M.X. François: Experimental and numerical investigation of heat transfer and flow within a thermoacoustic cell. International Journal of Transport Phenomena 6 (2004) 265–272. [Google Scholar]
 O. Hireche, C. Weisman, D. BalteanCarlès, P. Le Quéré, L. Bauwens: Low Mach number analysis of idealized thermoacoustic engines with numerical solution. Journal of the Acoustical Society of America 128, 6 (2010) 3438–3448. [Google Scholar]
 L. Ma, C. Weisman, D. BalteanCarlès, I. Delbende, L. Bauwens: Effect of a resistive load on the starting performance of a standing wave thermoacoustic engine: A numerical study. Journal of the Acoustical Society of America 138, 2 (2015) 847–857. [Google Scholar]
 Y. Fraigneau, N. de Pinho Dias, C. Weisman, D. BalteanCarlès: Numerical simulation of thermoacoustic heat pumping inside a compact cavity, in 16ème Congrès Français d’Acoustique, 11–15 April 2022, Marseille. 2022. [Google Scholar]
 S. Paolucci: On the filtering of sound from the NavierStokes equations. Report No. SAND828257, Sandia National Laboratories, 1982. [Google Scholar]
 R. Klein: Multiple spatial scales in engineering and atmospheric low Mach number flows. ESAIM: Mathematical Modelling and Numerical Analysis 39 (2005) 537–559. [CrossRef] [EDP Sciences] [Google Scholar]
 K. Goda: A multistep technique with implicit difference schemes for calculating two or threedimensional cavity flows. Journal of Computational Physics 30 (1979) 76–95. [Google Scholar]
 R. Knikker: A comparative study of highorder variableproperty segregated algorithms for unsteady low Mach number flows. International Journal for Numerical Methods in Fluids 66 (2011) 403–427. [CrossRef] [Google Scholar]
Cite this article as: Fraigneau Y. Weisman CB. & BalteanCarlès DB. 2023. A numerical model of thermoacoustic heat pumping inside a compact cavity. Acta Acustica, 7, 15.
All Tables
All Figures
Figure 1 Sketch of the 2D TACOT configuration. The computational domain is restricted to the upper part above the symmetry plane. Circle tags show probe locations used for recording time series of temperature. 

In the text 
Figure 2 Mean temperature difference ΔT between stack ends computed by numerical simulation for several values of c_{ps}, configuration case 1. 

In the text 
Figure 3 Mean temperature difference ΔT between stack ends for 6 values of κ_{s}. Numerical results case 1 (blue dots), linear theory (solid line). 

In the text 
Figure 4 Normalized mean temperature difference ΔT/ΔT_{ref} between stack ends as a function of the Drive ratio Dr. Numerical results (cases 1, 2, 3), linear theory (black solid line), Exp. 1 (green crosses) and Exp. 2 (red triangles) from [15]. 

In the text 
Figure 5 Normalized mean temperature difference ΔT/ΔT_{max} between stack ends as a function of the normalized stack velocity scale u_{s}/u_{opt}. 

In the text 
Figure 6 Normalized mean temperature difference ΔT/ΔT_{max} as a function of Φ_{u} − Φ_{p}. 

In the text 
Figure 7 Time evolution of the temperature difference between stack ends for the different configurations. 

In the text 
Figure 8 Mean temperature field over the computational domain for cases 1, 2 and 3a. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.