© The Author(s), Published by EDP Sciences, 2023

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 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–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, thermo-physical 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 heat-conductive Navier–Stokes equations, namely the conservation laws of mass, momentum, energy and the ideal state equation:

$$\frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}+\nabla \cdot \left(\rho \mathbf{V}\right)=0,$$(1)

$$\rho \left(\frac{\mathrm{\partial }\mathbf{V}}{\mathrm{\partial }t}+\mathbf{V}\cdot \nabla \mathbf{V}\right)=-\nabla P+\nabla \cdot \mathbf{\tau }+{\mathbf{f}}_{\mathbf{V}},$$(2)

$$\rho c_p\left(\frac{\mathrm{\partial }T}{\mathrm{\partial }t}+\mathbf{V}\cdot \nabla T\right)=\nabla \cdot \kappa \nabla T+\frac{\mathrm{d}P}{\mathrm{d}t}+\mathrm{\Phi },$$(3)

$$P=\rho \frac{R}{\mathcal{M}}T,$$(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 $\tau =\mu \left(\nabla \mathbf{V}{+}^T\right(\nabla \mathbf{V})-\left(\frac{2}{3}\nabla \cdot \mathbf{V}\right)\mathbb{I})$, 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⁻¹.kg⁻¹) and $\mathcal{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 λ₀ 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=\frac{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 $\mathcal{O}\left(M^2\right)$. 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:

$$\frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}+\nabla \cdot \left(\rho \mathbf{V}\right)=0,$$(5)

$$\rho \left(\frac{\mathrm{\partial }\mathbf{V}}{\mathrm{\partial }t}+\mathbf{V}\cdot \nabla \mathbf{V}\right)=-\nabla P_{\mathrm{dyn}}+\nabla \cdot \tau ,$$(6)

$$\rho c_p\left(\frac{\mathrm{\partial }T}{\mathrm{\partial }t}+\mathbf{V}\cdot \nabla T\right)=\nabla \cdot \kappa \nabla T+\frac{\mathrm{d}{P}_{\mathrm{th}}}{\mathrm{d}t},$$(7)

$$P_{\mathrm{th}}=\rho \frac{R}{\mathcal{M}}T.$$(8)

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:

$$P_{\mathrm{th}}=\frac{R}{\mathcal{M}}\frac{{\int }_{\mathrm{\Omega }}\mathrm{}\rho \mathrm{d}\mathrm{\Omega }}{{\int }_{\mathrm{\Omega }}\mathrm{}\frac{1}{T}\mathrm{d}\mathrm{\Omega }},$$(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 thermo-physical properties of the stack plates:

$${\rho }_s{c}_{\mathrm{ps}}\frac{\mathrm{\partial }T}{\mathrm{\partial }t}=\nabla \cdot {\kappa }_s\nabla T.$$(10)

Continuity of heat flux is applied at fluid-solid 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 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:

$${\nabla }^2\cdot \mathrm{\Psi }=\frac{\nabla \cdot \left({\rho }^{n+1}{\mathbf{V}}^{\mathrm{*}}-\frac{\mathrm{\partial }{\rho }^{n+1}}{\mathrm{\partial }t}\right)}{\Delta t},$$(11)where $\mathrm{\Psi }=P_{\mathrm{dyn}}^{n+1}-P_{\mathrm{dyn}}^n$ 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), $\frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}$ its time derivative and Δt the numerical time step. In order to satisfy the mass conservation Equation (5), the velocity field is updated as:

$${\mathbf{V}}^{n+1}={\mathbf{V}}^{\mathrm{*}}-\frac{\Delta t}{{\rho }^{n+1}}\nabla \mathrm{\Psi }.$$(12)

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.

2. 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).

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.

2.3 Boundary conditions

2.3.1 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.

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:

$$Q_m\left(i\right)={\rho }_0{U}_i{S}_i\sin (2\pi f_{0}t-{\varphi }_i),$$(13)where Q_m(i) is the mass flow rate at each inlet i with an amplitude ρ₀U_iS_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₀ 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 back-to-back inlets which have the same parameters ρ₀, 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 ρ₀, 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 ${\mathrm{\Phi }}_u-{\mathrm{\Phi }}_p=\frac{\pi }{2}-{\varphi }_{\mathrm{OFG}}\pm m\cdot \pi$ (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 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 ρ_s = 2000 kg.m⁻³. and the thermal conductivity is set to κ_s = 0.2 W.m⁻¹.K⁻¹ (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⁻¹.kg⁻¹). 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⁻¹.kg⁻¹, κ = 8.56·10⁻² W.m⁻¹.K⁻¹ and μ = 2.11·10⁻⁵ kg.m⁻¹.s⁻¹ respectively.

The frequency of both the acoustic driver and the oscillating fluid generator is fixed at f₀ = 50 Hz. By considering a reference temperature of T₀ = 298 K and mass density of ρ₀ = 0.597 kg.m⁻³ (corresponding pressure P₀ = 10⁵ Pa), the thicknesses of the dynamic and thermal boundary layers along the stack plates are estimated to be δ_ν = 4.75·10⁻⁴ m and δ_κ = 8.07·10⁻⁴ 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 κ₁ = 0.25 W.m⁻¹.K⁻¹ (case 3a : moderate conductivity specific to ceramics [15]) or κ₂ = 14 W.m⁻¹.K⁻¹ (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⁻⁵ 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⁻⁴ m).

The time step is fixed at Δt = 2·10⁻⁵ 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₀ = 10⁵ Pa and T₀ = 298 K and the mass density is equal to ρ₀ = 0.597 kg.m⁻³.

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 half-distance 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:

$$\mathrm{Dr}=\frac{P_{\mathrm{th}}^{\max }-P_{\mathrm{th}}^{\min }}{P_{\mathrm{th}}^{\max }+P_{\mathrm{th}}^{\min }}.$$

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 $u_s=\frac{U_{\mathrm{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 Δ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 quasi-convergence, 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, 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⁻¹ (U_OFG = 0.25 m.s⁻¹ 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 thermo-physical properties

Results shown in this paragraph were obtained for the configuration case 1. The considered thermo-physical properties are the heat capacity ρ_sc_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 ρ_sc_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⁻¹.kg⁻¹. 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 cps, 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 (non-zero) 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⁻¹ 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/ΔTref 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 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 = 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 $u_s^{\mathrm{*}}=u_s/u_{\mathrm{opt}}$ 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 $u_s^{\mathrm{*}}$, up to $u_s^{\mathrm{*}}$ about 0.85, then reaches the same plateau. There are however some quantitative discrepancies for $u_s^{\mathrm{*}}$ < 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 $u_s^{\mathrm{*}}$ decreases. However the normalization allows to obtain analogous behavior for all these physically very different configurations.

Figure 5 Normalized mean temperature difference ΔT/ΔTmax between stack ends as a function of the normalized stack velocity scale us/uopt.

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⁻¹, 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 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 Φ 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/ΔTmax 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 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.

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 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 Δ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 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.

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 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.

Acknowledgments

This research was financially supported by Agence Nationale de la Recherche (Projet ANR-17-CE06-0 0 07-01).

References

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 cps, 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/ΔTref 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/ΔTmax between stack ends as a function of the normalized stack velocity scale us/uopt.
In the text

	Figure 6 Normalized mean temperature difference ΔT/ΔTmax 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