Numerical study of the effects of natural convection in a thermoacoustic configuration

. This study focuses on natural convection flows within a cylindrical guide containing a porous medium. This configuration is applicable to standing-wave thermoacoustic engines, usually composed of an acoustic resonator where a (short) stack (or porous medium) is inserted, with a heat exchanger placed at one of its ends. The resulting horizontal temperature gradient, when high enough, triggers the onset of an acoustic wave. Natural convection effects are usually neglected in thermoacoustics so that axisymmetry is often assumed. Here a 3D numerical study of natural convection flow is performed using a finite volume code for solving mixed Navier-Stokes and Darcy-Brinkman equations under Boussinesq approximation. The influence of the porous medium’s physical characteristics (permeability, thermal conductivity, anisotropy) on the flow and temperature fields is investigated. It is shown that such flows are fully three-dimensional and therefore can modify significantly starting as well as steady operating conditions of the thermoacoustic engine.


Introduction
In thermoacoustic devices, dynamic and thermal interactions between acoustic oscillations of a working gas and solid walls generate either heat pumping for refrigerating purposes or mechanical work for engine operating purposes [1]. These devices represent an ecological alternative for industrial applications but they suffer from poor performance. This is generally caused by the existence of several nonlinear phenomena interfering with the process of energy transfer and penalizing the global performance of the device.
The energy conversion takes place in the vicinity of solid-fluid interfaces and in order to increase the exchange surface, a porous medium is used for this conversion. The porous medium is a significant component of thermoacoustic devices and its characteristics greatly influence the device performance. The pores are usually of various shapes and sizes. In order to maintain the temperature distribution along the stack stable, thermal conductivity of the solid matrix should be lower than that of the working gas. Arnott et al. [2] developed a general thermoacoustic linear theory for the elements of thermoacoustic devices * e-mail: weisman@limsi.fr that have a microstructure and compared their performances. The porous media considered were solid stacks of plates or solids with many capillary tubes of circular, square or triangular geometry. They concluded that heat and work flows are approximately 10 percent greater for the parallel plate stack geometry than for other pore geometries. This model is still used in DeltaEC software [3], which is the reference code for computing a complete thermoacoustic system.
In thermoacoustic systems several nonlinear phenomena can take place, reducing their performance, such as: jet streaming, acoustic streaming, end effects and natural convection. These effects are still under investigation. Among them, acoustic streaming has already been the subject of many studies during the last decades [4][5][6]. Natural convection is a phenomenon that has mostly been neglected until now in thermoacoustic applications, although it does exist because of the presence along the resonator of heat exchangers at different temperatures: two heat exchangers are generally placed at each end of the porous medium, in order to maintain a temperature gradient along the porous medium in thermoacoustic engines, or to benefit from heat pumping effects. Therefore, if the device is horizontal, there will be horizontal temperature gradients along the porous medium saturated with gas, as well as along the two cavities filled with gas on both sides of the porous medium. Natural convection flows can then develop as in differentially heated enclosures.
Natural convection in fluids and porous cavities has been treated extensively in the literature (see for example [7,8]). A more recent study has been conducted by Weisman et al. [9] in a geometry of a thermoacoustic engine. The goal was to predict the temperature distribution in a 2D differentially heated horizontal cavity, filled with a stack of plates which constitutes an anisotropic porous medium. This study showed that a darcean porous model gives correct orders of velocity scales and allows for a qualitative description of the solutions. Some recent experimental studies of natural convection in thermoacoustic configurations show the influence of natural convection on Rayleigh streaming [10,11].
In the present work, natural convection flows are investigated in a simplified geometry of a standing-wave thermoacoustic engine, in order to assess the flow and temperature distributions and their order of magnitude. A 3D numerical study of natural convection flow is performed using a finite volume code for solving Navier-Stokes equations under Boussinesq approximation in the fluid domain. In the porous medium, the equations are modified using a Darcean term, according to Darcy-Brinkman formulation [7,8]. The Darcean term becomes dominant as permeability decreases. The influence of the porous medium's physical characteristics (permeability, thermal conductivity, anisotropy) on the flow and temperature fields is investigated. The results are then analysed in order to explore thermo-physical characteristics adequate for thermoacoustic devices [12].

Geometry and boundary conditions
A simplified geometry for a standing-wave thermoacoustic engine is introduced: a horizontal cylindrical tube of length L and inner radius R, closed at both ends, is filled with air. A cylindrical porous medium of length L p and radius R is inserted in the tube. In a standing-wave thermoacoustic engine, the fundamental mode is of acoustic wavelength λ = 2L. The porous medium length L p is usually small compared to the wavelength, such that L p L. The computational domain is shown in Figure 1.
Attached to the left end of the porous medium is a hot heat exchanger (positioned at a distance L L from the left end of the tube), which can be a coiled resistance modeled here as solid rings of width e h , maintained at given (hot) temperature T h . Usually the "hot" part of the tube (left of coiled resistance) is insulated and therefore in the present model, all walls in this part are considered adiabatic. There is no cold heat exchanger in the present model, although there usually is one in experiments. Here the "cold" temperature is the ambient air temperature outside the cylinder, noted T c . In the present model this cold temperature is imposed on all walls at the right of the porous medium. Finally, a linear temperature profile is imposed on the tube wall around the porous medium. Thermal boundary conditions are shown in Figure 1 (bottom). At all solid-fluid interfaces, we assume continuity of temperature and heat flux as well as no slip conditions. The porous-fluid interface is not treated as such since the equations are solved in both porous and fluid domains (see next section), a phase variable allowing to select the appropriate equation. Continuity of velocity, temperature, heat flux and stresses are implicitely ensured with the use of a staggered mesh.
The porous medium is usually positioned inside the tube at about 1/8 of the total length from the left end. However, in this natural convection study, with the imposed boundary conditions, we have shown that the gas in the right half of the tube is almost stagnant at cold temperature. Therefore it is not necessary to extend the computational domain too much. In the following simulations the porous medium will be positioned at about 1/4 of the total tube length.
The porous medium, consisting of its solid matrix and the fluid pores, is modeled as an equivalent homogeneous porous medium, with porosity , permeability tensor K, thermal conductivity tensor k p . Also, it is common in thermoacoustics to use porous media composed of a solid (ceramic) matrix with horizontal cylindrical pores. In that case the porous medium is strongly anisotropic and this possibility will also be considered.

Governing equations
In the fluid domain the Navier-Stokes equations are classically written under Boussinesq assumptions. The reference gas density ρ, thermal conductivity k, heat capacity c p , dynamic viscosity µ and volume expansion coefficient at constant pressure β = 1/T for an ideal gas are constant and equal to their value for cold temperature T c .
In the porous medium, the momentum equations are modified using a Darcy term. In the resulting Darcy-Brinkman equation, two simplifications are made related to the fact that in the application, porosity is rather large (around 0.8): the 1/ factor is omitted in the unsteady term (which does not affect the steady solution) and the effective viscosity is close to the fluid viscosity µ eff = µ/ µ [7]. The advection and Forchheimer terms are also omitted in the formulation since the pore Reynolds number Re p is assumed to be very small (Re p 1). This assumption will be verified on results. The energy equation in the porous medium is deduced from the fluid and solid energy equations averaged over an elementary volume [7]. The value of (ρc p ) p in the porous medium, given by (ρc p ) p = (ρc p ) + (1 − )(ρc p ) s , does not affect steady state. A phase variable X is introduced to identify the nature of the medium (X = 1 for gas, X = 0 for the porous medium). Therefore the equations can be written in a unified form as: In the above equations V is the velocity vector (equal to the seepage velocity in the porous medium), p is the nonhydrostatic pressure, T is the temperature, t is time, e x is the unit vector in the ascending vertical direction x, g is the acceleration of gravity and I the identity tensor.

Numerical solution
Computations are performed with the code SUNFLUIDH, which has been developped at LIMSI for simulations of unsteady incompressible or low Mach number flows.
The previous governing equations are solved using a finite volume approach with a 2nd order discretization in space and time. An implicit treatment of the viscous and thermal conductive terms ensures better numerical stability in regard to the time step which only depends on the CFL criterion. The equation set is solved by means of an ADI method [13]. An incremental projection method is used in order to ensure the divergence-free velocity field [14,15]. This implies to solve a Poisson's equation at each time step which is carried out with a multigrid-SOR algorithm [16]. The solution corresponds to the time-variation of pressure from which the correction term is estimated in order to get the divergence-free velocity field.
The mesh is body-fitted irregular, with 64 points regularly spaced along the radial direction, 48 points regularly spaced along the azimutal direction and 256 points placed along the axis, with refined spacing near the heat exchanger and in the porous medium. The minimum and maximum mesh sizes in the axial direction are ∆z min = 3.9 × 10 −4 m and ∆z max = 4.2 × 10 −3 m. In the radial direction, ∆r = 3.9 × 10 −4 m. This spatial resolution is adequate for the given problem configurations.
The time step is adjusted throughout the simulation in order to satisfy a CFL<0.3 criterion.

Reference configuration
The geometrical and thermal characteristics chosen for all cases are shown in Table 1. The thermophysical characteristics for air (fixed) and for the porous medium (for the reference case) are shown in Table 2. These parameters correspond (mostly) to some existing experiments in thermoacoustics [11]. The porous medium is considered to be isotropic and the permeability and thermal conductivity tensors are reduced to scalar values (K = KI and k p = k p I). The permeability value 10 −6 m 2 constitutes a high limit to practical values found in the literature [12]. If the porous medium is made of horizontal capillaries of square cross-section, the permeability is K = a 2 /8 corresponding to a hydraulic radius a equal to 3 mm. The equivalent thermal conductivity k p is calculated from that of the solid matrix constituting the porous medium (k s ), the gas thermal conductivity and the porosity: The nondimensional parameter controlling natural convection is the Rayleigh number, which is classically calculated as Ra = ρgβ∆T H 3 µα , with H a reference height, ∆T the imposed temperature difference and α = k ρcp the fluid thermal diffusivity [7]. The associated vertical velocity scale is U ref = α H Ra 1/2 . For the reference configuration, with H = R and ∆T = 207 K, we find that Ra = 2 × 10 5 and U ref = 0.48 m/s. The solution should thus be stationary. However as we will see in the results section, several convection flows can be distinguished within the domain, and this choice of Rayleigh number may not be representative of the overall flow. Note that for this value of ∆T the ratio ∆T /(2T ref ) = 0.3 is rather large and falls outside of the validity range for the Boussinesq approximation (for air). Non Boussinesq simulations would require the use of a low Mach number model which is harder to implement. Since the flow is steady, expected differences should be small [17]. Therefore Boussinesq approximation is used throughout this paper.
In a porous cavity, the flow slows down and therefore a modified Rayleigh number is introduced, Ra p = ρgβ∆T HK µαp , with α p = kp ρcp the porous medium equivalent thermal diffusivity [7]. The associated vertical velocity scale is U refp = αp H Ra p . For the reference configuration, using the same values for ∆T and H = R, this yields Ra p = 6.71, U refp = 0.35 m/s and we expect a stationary solution of the natural convection problem. Again these reference values are just indicators since the problem here is a mixed fluid and porous domain with prescribed temperature differences unevenly distributed.

Fig. 2. Reference configuration, vertical middle plane section: (a) Velocity vectors with contours of axial velocity (m/s). (b)
Contours of temperature (in K), white dotted lines show the extent of the porous medium.

Numerical results for the reference configuration
The flow and temperature fields are first described in the reference configuration presented in the previous section. They are fully three-dimensional. Figure 2 shows on the vertical middle plane section the velocity vectors with contours of axial velocity (Fig. 2a) and contours of temperature (Fig. 2b). White dotted lines show the extent of the porous medium. The flow and temperature fields on this plane are characteristic of natural convection in a differentially heated fluid cavity on the left of the heat exchanger (with adiabatic conditions on all boundaries except for the heat exchanger section). On the right of the heat exchanger, the flow and temperature field resemble that of natural convection in a differentially heated porous-fluid cavity (with prescribed temperature on all boundaries). However the heat exchanger is not a wall and therefore there is fluid passing from one cavity to the other. Also the heat exchanger is not ideal and the condition T = T h is not imposed homogeneously on the cross-section. Figure 3 shows on several vertical cross-sections the contours of vertical velocity (left) and of temperature (right). The corresponding cross sections are shown on the sketch at the bottom of Figure 3. The flow and temperature fields at z = z 1 (Fig. 3ab1), on the left of the heat exchanger, are characteristic of natural convection flow in a differentially heated cylindrical cavity, with adiabatic boundaries. On a section cutting through the heat exchanger, at z = z 2 (Fig. 3ab2) the flow is constrained between the solid rings that constitute the heat-exchanger and hot temperature T h is maintained on the rings. On a section cutting through the porous medium, at z = z 3 (Fig. 3ab3), the flow and temperature field are essentially characteristic of natural convection in a differentially heated porous cylindrical cavity, with fixed temperature on the tube wall. Finally, on the right part of the porous medium, at z = z 4 (Fig. 3ab4) (8K/ ) = 3 mm, this yields the pore Reynolds number Re p 1.7. This value is actually border line for justifying the use of Darcy-Brinkman model without the introduction of the quadradic Forchheimer drag and the advective terms. The transition region between models corresponds to Re p between 1 and 10 [7].
For this reference configuration, the boundary layer thickness can be estimated for the porous domain using δ ν δ κ HRa −1/2 K = 9.6 × 10 −3 m [7]. Therefore since the mesh is refined near the porous medium's left and right sections there should be about 20 mesh points within the boundary layer, which is plenty. In the fluid domain on the right of the porous medium, the temperature difference is much lower (it can be estimated to be about 25 K, based . There should be about 5 points within the boundary layer, which is sufficient. The reference configuration being most restrictive (for all other cases considered in the next section the boundary layer thickness is larger) numerical results can be considered as converged.

Parametric study
In this section, several parameters characterizing the porous medium are varied, keeping a fixed porosity: permeability, anisotropy, equivalent thermal conductivity.
In order to analyse the influence of each parameter, we define several flow and temperature output parameters: The maximum vertical velocity amplitude in the tube is noted U max and (x max , y max , z max ) U is the corresponding location. The maximum axial velocity amplitude in the tube is noted W max and (x max , y max , z max ) W is the corresponding location. The total kinetic energy within the tube is noted E c . Also, the ratio of vertical versus horizontal temperature gradient in the core of the porous medium is noted γ = | ∂T ∂x |/| ∂T ∂z |. It is a measure of the strength of natural convection: if conduction is dominant γ → 0 and if convection is dominant γ → ∞.
The value of the permeability K depends on the geometry of the porous medium. Usually in thermoacoustic configurations, the porosity is fairly large (between 0.75 and 0.85), and the permeability is kept small to maximize surface exchange (between 10 −6 m 2 and 10 −9 m 2 ). In the following paragraph, the permeability is varied between these values, considering the porous medium to be homogeneous and isotropic with constant porosity = 0.84. The hydraulic radius a corresponding to each value of K can be estimated from the relationship K = a 2 /8 (for capillaries of square cross-section). The hydraulic radius varies between 3 mm (K = 10 −6 m 2 ) to 0.1 mm (K = 10 −9 m 2 ). Table 3 presents the results for 4 different values of the porous medium permeability, when the porous medium is isotropic. As K decreases, the total kinetic energy and the γ ratio also decrease, which is to be expected because the porous medium acts more and more like a solid barrier and natural convection flow looses strength. It is interesting to note that the maximum vertical velocity decreases and that its location shifts from just right of the heat exchanger (on the central horizontal plane, facing the middle of the outer heat exchanger fluid space) to just outside the porous medium as permeability decreases (slightly higher than the tube center). The maximum axial velocity also decreases with K and its location (in the right fluid cavity, on the central vertical plane, facing the middle of the outer heat exchanger fluid space) shifts to just outside the porous medium as permeability decreases. This is consistent with natural convection flow developing outside the porous medium when K is very small. The orders of magnitude of the maximum velocities vary with the value of K, but not proportionally as in a simpler configuration with natural convection through a porous medium (see the reference velocity scales U ref and U refp presented in the previous section).   Figure 4 shows on the vertical middle plane section, for K = 10 −7 m 2 , the velocity vectors with contours of axial velocity (Fig. 4a) and contours of temperature (Fig. 4b).
White dotted lines show the extent of the porous medium. Although velocities are much smaller than in the reference case, the effect of natural convection on the temperature distribution is clear and a much lesser fraction of the porous medium is submitted to a horizontal temperature gradient. Therefore the thermoacoustic energy conversion should be lowered significantly. In this configuration the hydraulic radius is on the order of 1 mm, the average velocity amplitude estimated in the porous medium is about 5.5 mm/s and the pore Reynolds number is 0.3, which validates the use of Darcy-Brinkman equation. Figure 5 shows on the vertical middle plane section, for K = 10 −8 m 2 the velocity vectors with contours of axial velocity (Fig. 5a) and contours of temperature (Fig. 5b). White dotted lines show the extent of the porous medium. It is clear on this figure that the porous medium behaves almost like a conducting solid. On the right end of the porous medium, the temperature is only slightly higher than T c , so that convection is greatly diminished and the maximum velocities U max and W max are much smaller than in the reference case (on the order of 0.05 m/s, which corresponds to the fluid velocity scale U ref for H = R and ∆T = 2 K). On the left of the heat exchanger, the temperature is almost homogeneous (equal to T hot ) and there is nearly no flow. Similar observations can be made for K = 10 −9 m 2 . The pore Reynolds number are 10 −2 for K = 10 −8 m 2 and 3 × 10 −4 for K = 10 −9 m 2 , which validates the use of Darcy-Brinkman equation.
The equivalent porous medium thermal conductivity is varied next, for homogeneous isotropic porous media, keeping all other parameters corresponding to the reference case. Four values were tested, from k p = 0.028 to k p = 1.4 W.m −1 .K −1 , characteristic of porous media used in thermoacoustics [12]. It was shown that the flow and temperature field variations are negligible within this range of values.
Anisotropy is considered next: in the case of pores that are horizontal capillary tubes, the porous medium is almost impermeable in the transverse direction, so that K xx 0. The permeability in the axial direction K zz is varied between 10 −6 m 2 and 10 −9 m 2 . Throughout this study the equivalent porous medium thermal conductivity is assumed to be homogeneous and isotropic. Figure 6 shows on the vertical middle plane section the velocity vectors with contours of axial velocity (Fig. 6a) and contours of temperature (Fig. 6b), for the anisotropic case, K zz = 10 −6 m 2 (all other parameters being maintained the same as in the reference configuration). The flow structure is strongly modified inside the porous medium, since it is constrained horizontally and flows through the fluid space between the heat exchanger rings (compare with Fig. 2a). The temperature field is also modified (compare with Fig. 2b). The solution inside the porous medium is analogous to natural convection flow in a shallow horizontal cavity, with on the top half hot fluid flowing to  the right and on the bottom half cold fluid flowing to the left. Then on the right of the porous medium, a more standard natural convection fluid cavity flow is recovered. For this configuration, the maximum vertical velocity is 0.12 m/s measured at (x = R/2, y = 0, z = 0.151 m) (just on the right of the porous medium) and the maximum axial velocity is 0.17 m/s measured at (x = 9R/10, y = 0, z = 0.165 m) (also on the right of the porous medium). The total kinetic energy is E c = 2.84 × 10 −4 J (smaller than in the reference configuration) but the ratio γ is equal to γ = 4.56 (larger than in the reference configuration).

Consequences for thermoacoustic applications
The previous simulations of natural convection flows were performed on configurations (geometry, materials, applied temperature differences) applicable to thermoacoustic engines, even though there was no acoustics.
Based on these simulations and observations, steady natural convection flows should not be neglected in the design of standing-wave thermoacoustic engines. The onset of an acoustic wave is based on the existence of a temperature gradient regularly distributed along the porous medium. Since natural convection flows induce very distorted temperature fields, they can delay or even prevent the onset of acoustic waves. Natural convection flows could interfere with acoustics and also with other nonlinear phenomena occuring in such thermoacoustic devices, such as streaming flow. It would be interesting to run some simulations of acoustics and natural convection together.
Finally such flows may have applications for other devices involving thermoacoustics, such as thermoacoustic heat pumps or traveling-wave thermoacoustic engines, as long as there is a horizontal temperature difference between the constituting elements.