### 3.2.4. Problem 3: third source term

We now consider the third source term, yielding the following equation

$$ \frac{\pi}{2} c_0 \hat{L} \frac{\partial^3 \bar{u}_3}{\partial \hat{y}^3} = - \frac{\partial^2 (\overline{u'v'})}{\partial \hat{x}^2}, \quad (3.31) $$

with boundary conditions (3.11a,b)–(3.11d) for $\bar{u}_3$. Integrating (3.31) gives the simplified solution

$$ \bar{u}_3(\hat{x}, \hat{y}) = \frac{\pi^2 U_{ac}^2}{480 c_0 \hat{L}^2} \sin(2\pi \hat{x}) \left\{ \begin{aligned} & 4\hat{R}^3 + 150 - \frac{585}{\hat{R}} + \frac{720}{\hat{R}^2} \\ & -3\hat{y}^2 \left( \frac{240}{\hat{R}^4} - \frac{315}{\hat{R}^3} + \frac{150}{\hat{R}^2} + 8\hat{R} \right) + 60e^{2(\hat{y}-\hat{R})} + 20 \frac{\hat{y}^4}{\hat{R}} \\ & + 20e^{\hat{y}-\hat{R}} \left( \left( 12 - \frac{18}{\hat{R}} \right) \cos(\hat{R}-\hat{y}) - \frac{18}{\hat{R}} \sin(\hat{R}-\hat{y}) + 12 \frac{\hat{y}}{\hat{R}} \sin(\hat{R}-\hat{y}) \right) \end{aligned} \right\}. \quad (3.32) $$

Along the axis we obtain, neglecting terms of order $O(e^{-\hat{R}})$ and smaller

$$ \frac{\bar{u}_3(\hat{x}, 0)}{\bar{u}_R} = -\frac{2\pi^2}{45} \left(\frac{\hat{R}}{\hat{L}}\right)^2 \left[ \hat{R} + \frac{75}{2} - \frac{585}{4\hat{R}} + \frac{180}{\hat{R}^2} \right]. \quad (3.33) $$

Since the term in brackets is always positive, velocity $\bar{u}_3$ is of the opposite sign to Rayleigh's solution. Thus the associated vortices are reversed.

Note that in this equation the dominant term $\hat{R}(\hat{R}/\hat{L})^2$ could increase significantly for large values of $\hat{R}$. However, our working hypotheses are that the linear regime of streaming is respected, that is $Re_{NL} = (M\hat{R})^2 \ll 1$, and also that $(\hat{R}/\hat{L})^2 \ll 1$. These relations imply that the product $\hat{R}(\hat{R}/\hat{L})^2$ and consequently the velocity component $\bar{u}_3$ will both remain bounded for large values of $\hat{R}$.

The second-order pressure gradient in the $x$ direction can be calculated from

$$ \frac{\partial \bar{p}_3}{\partial \hat{x}} = \rho_0 \frac{\pi}{2} c_0 \frac{\partial^2 \bar{u}_3}{\partial \hat{y}^2}, \quad (3.34) $$

yielding

$$ \frac{\partial \bar{p}_3}{\partial \hat{x}} = -\frac{\rho_0 \pi^3 U_{ac}^2}{2 \hat{L}^2} \sin(2\pi \hat{x}) \left\{ \begin{aligned} & \frac{3}{\hat{R}^4} - \frac{63}{16\hat{R}^3} + \frac{15}{8\hat{R}^2} + \frac{\hat{R}}{10} - \frac{\hat{y}^2}{2\hat{R}} - \frac{1}{20}e^{2(\hat{y}-\hat{R})} \\ & + \frac{1}{2}e^{\hat{y}-\hat{R}} \left( \frac{-1+2\hat{y}}{\hat{R}} \cos(\hat{R}-\hat{y}) + \left(\frac{1}{\hat{R}}-2\right) \sin(\hat{R}-\hat{y}) \right) \end{aligned} \right\}. \quad (3.35) $$

For this source term the axial pressure gradient (which is very small) does depend on $\hat{y}$.
### 3.2.5. Problem 4: fourth source term
Considering the fourth source term yields

$$ \frac{\pi}{2} c_0 \hat{L} \frac{\partial^3 \bar{u}_4}{\partial \hat{y}^3} = - \hat{L} \frac{\partial^2 (\overline{v'v'})}{\partial \hat{x} \partial \hat{y}}, \qquad (3.36) $$

with boundary conditions (3.11a,b)–(3.11d) for $\bar{u}_4$. Integrating equation (3.36) and using the boundary conditions, the following simplified solution is found:

$$ \bar{u}_4(\hat{x}, \hat{y}) = \frac{\pi^2 U_{ac}^2}{120c_0} \sin(2\pi\hat{x}) \left(\frac{\hat{R}}{\hat{L}}\right)^2 \left[1 - 6\left(\frac{\hat{y}}{\hat{R}}\right)^2 + 5\left(\frac{\hat{y}}{\hat{R}}\right)^4\right]. \qquad (3.37) $$

Along the axis we obtain the following expression for the normalized axial streaming velocity

$$ \frac{\bar{u}_4(\hat{x}, 0)}{\bar{u}_R} = \frac{2\pi^2}{45} \left(\frac{\hat{R}}{\hat{L}}\right)^2. \qquad (3.38) $$

Comparing equations (3.38) and (3.33) shows that $\bar{u}_4(\hat{x}, 0)/\bar{u}_R$ is equal to $-1/\hat{R}$ times the first term of $\bar{u}_3(\hat{x}, 0)/\bar{u}_R$. The contribution of the fourth source term is thus negligible for large enough guides.

### 3.2.6. Superposition of solutions

Gathering together the four components $\bar{u}_1$, $\bar{u}_2$, $\bar{u}_3$ and $\bar{u}_4$, and neglecting the lower-order terms, the total streaming velocity along the axis is given by

$$ \frac{\bar{u}(\hat{x}, 0)}{\bar{u}_R} = \frac{\sum_{i=1}^{4} \bar{u}_i(\hat{x}, 0)}{\bar{u}_R} = 1 - \frac{29}{6} \frac{1}{\hat{R}} - \frac{2\pi^2}{45} \frac{\hat{R}^3}{\hat{L}^2} \left[ 1 + \frac{73}{2\hat{R}} - \frac{585}{4\hat{R}^2} + \frac{180}{\hat{R}^3} \right]. \qquad (3.39) $$

This expression shows that, under certain circumstances, the streaming flow can be reversed on the axis $(\bar{u}(\hat{x}, 0)/\bar{u}_R < 0)$. This occurs when

$$ \hat{L} < \hat{L}_{\text{limit}} = \pi \sqrt{\frac{2}{45}} \left( \frac{\hat{R}^3 \left( 1 + \frac{73}{2\hat{R}} - \frac{585}{4\hat{R}^2} + \frac{180}{\hat{R}^3} \right)}{1 - \frac{29}{6} \frac{1}{\hat{R}}} \right)^{1/2}. \qquad (3.40) $$

Figure 4 shows the variation with $\hat{R}$ of the limit value $\hat{L}_{\text{limit}}/\hat{R}$. For geometries that fall under the curve, there exists a reversed flow on the axis. Minimum values are $\hat{R} = 13.52$, $\hat{L}/\hat{R} = L/R = 5.24$. These values are compatible with the hypothesis $(\hat{R}/\hat{L})^2 \ll 1$. For example, for $\hat{R} = 100$, reversed outer streaming velocity occurs if $\hat{L} < 789$, or $L/R < 7.89$.

In order to characterize the possible applications, let us consider a practical case of reverse flow along the guide axis. Let $\hat{L}$ and $\hat{R}$ be fixed. For resonance conditions, $\lambda = 2L$ and $\omega = \pi c_0/L$. Using $\delta_v^2 = 2\nu/\omega$ and $\hat{L} = L/\delta_v$, the following equations are obtained

$$ L = \frac{2\nu}{\pi c_0} \hat{L}^2, \quad R = L \frac{\hat{R}}{\hat{L}}. \qquad (3.41a,b) $$
<figure><img src="image_127.png" /><figcaption>FIGURE 4. Limit value $\hat{L}_{limit}(\hat{R})/\hat{R}$ under which the flow is reversed on the axis.
Minimum values: $\hat{R} = 13.52$, $\hat{L}_{limit}/\hat{R} = 5.24$.</figcaption></figure>

<figure><img src="image_128.png" /><figcaption>FIGURE 5. (Colour online) Axial streaming velocities $\bar{u}_1(\hat{x}, \hat{y})/\bar{u}_R$ (red solid line), $\bar{u}_2(\hat{x}, \hat{y})/\bar{u}_R$ (blue dashed line), $\bar{u}_3(\hat{x}, \hat{y})/\bar{u}_R$ (green dotdashed line) and $\bar{u}_4(\hat{x}, \hat{y})/\bar{u}_R$ (purple dotted line). $\hat{R} = 50$, $\hat{L} = 500$. (a) $\hat{y} = 0$; (b) $\hat{x} = -1/4$.</figcaption></figure>

Let $\hat{L} = 789$ and $\hat{R} = 100$, as previously. Taking standard values corresponding to air, $\nu = 1.496 \times 10^{-5} \text{ m}^2 \text{ s}^{-1}$ and $c_0 = 343.82 \text{ m s}^{-1}$, we obtain $L = 17.2 \text{ mm}$, $R = 2.15 \text{ mm}$ and the resonance frequency $f = 10 \text{ kHz}$. These dimensions are small and could concern microfluidics systems, where Rayleigh streaming is used, for example within the field of acoustic particle trapping and manipulation or efficient fluid mixing. It is interesting to remark that some experiments in microfluidics have shown Rayleigh streaming vortices rotating in the opposite direction to that expected. Our results could give an explanation to this phenomenon that is a subject of debate (Wiklund, Green & Ohlin 2012).

Figure 5 shows the four components $\bar{u}_1/\bar{u}_R$, $\bar{u}_2/\bar{u}_R$, $\bar{u}_3/\bar{u}_R$ and $\bar{u}_4/\bar{u}_R$ along the axis (a) and along $\hat{y}$ for $\hat{x} = -1/4$ (b). Here we have $\hat{R} = 50$ and $\hat{L} = 500$. We notice that the fourth contribution is negligible, and that the third contribution is of the same order of magnitude as the second one. Figure 6 shows the total axial velocity $\bar{u}$ profile for $\hat{R} = 100$ and two values of $\hat{L}$: $\hat{L} = 600$ and $\hat{L} = 8000$. The difference between the two profiles is important in the core of the guide, whereas in the very near wall region, both profiles are very close.

This feature was also observed (even though the acoustic regime and the geometric configuration were different) in experiments as well as numerical simulations, that revealed situations where the outer streaming flow was greatly modified (up to occurrence of reverse flow), while the inner streaming flow was nearly unmodified (Reyt et al. 2013; Reyt, Bailliet & Valière 2014).
<figure><img src="image_129.png" /><figcaption>FIGURE 6. Total normalized axial streaming velocity transverse profile $u(-1/4, \hat{y})/\bar{u}_R$ for $\hat{R}=100$ and $\hat{L}=600$ (solid line), $\hat{L}=8000$ (dashed line).</figcaption></figure>

At this stage we have shown that the usually neglected third source term of outer streaming is responsible for drastic change of behaviour of acoustic streaming for certain ($\hat{L}$, $\hat{R}$) conditions. The ratio $R/\delta_\nu$ is usually considered as the criterion to differentiate between different streaming regimes, from narrow guides with only inner cells to large guides with mostly outer cells. The present study implies that not only $R/\delta_\nu$ but also $R/L$ ratios should be considered in order to examine the organization of the streaming flow.

In the following sections we consider the impact of this finding on other characteristics of streaming flow: slip velocity, mass transport velocity, vertical component of streaming velocity.

### 3.2.7. Slip velocity

The viscous boundary layer vortex thickness (inner streaming) was estimated to be approximatively equal to $1.9\delta_\nu$ by Schlichting (1932), who did not consider streaming outside the boundary layer. As shown by experiments (Moreau *et al.* 2008), the limit for large guides between inner and outer streaming cells is $3\delta_\nu$ from the wall, as illustrated in figure 1. At $3\delta_\nu$ from the guide wall, the axial streaming velocity was indeed found to reach its maximum. The flow in the core of the guide is initiated by transverse diffusion of momentum from there. This position is the limit between inner and outer streaming cells and the value of the axial streaming velocity at non-dimensional $\hat{y}=\hat{y}^l=\hat{R}-3$ is thus generally considered as a slip velocity when inner streaming is not described. Rayleigh's solution for the slip velocity is given by $\bar{u}_R^{SV} = (3/8)U_{ac}^2/c_0 \sin(2\pi\hat{x})$.

The actual non-dimensional transverse location $\hat{y}^{SV}$ of the maximum of the axial streaming velocity can however deviate from this standard value, depending on values of $\hat{R}$ and $\hat{L}$. Figure 7 shows the difference (or offset) $\hat{y}^{SV} - \hat{y}^l$ between the actual and standard location of the slip velocity as a function of $\hat{R}$ (for $6 < \hat{R} < 200$), for two values of $\hat{L}$. For the largest value of $\hat{L}$, we note that the actual location shifts slightly away from the wall as $\hat{R}$ increases (it is placed at the standard value for $\hat{R} = 28$). This is in agreement with previous results showing the evolution of the axial streaming velocity transverse profile with the ratio $R/\delta_\nu$ (e.g. figure 4 of Bailliet *et al.* (2001) and figure 4 of Hamilton *et al.* (2003a)).
<figure><img src="image_130.png" /><figcaption>FIGURE 7. Difference (or offset) between the actual location $\hat{y}^{SV}$ and the standard location $\hat{R} - 3$ of the slip velocity, as a function of $\hat{R}$. Solid line: $\hat{L} = 8000$, dashed line: $\hat{L} = 500$.</figcaption></figure>

<figure><img src="image_131.png" /><figcaption>FIGURE 8. Normalized axial slip velocity at position $\hat{y} = \hat{y}^{SV}$, as a function of $\hat{R}$. Solid line: $\hat{L} = 8000$, dashed line: $\hat{L} = 500$.</figcaption></figure>

For the smallest value of $\hat{L}$, the actual location of the slip velocity shifts closer to the wall as $\hat{R}$ increases. This can be related to the fact that the $\bar{u}_3$ contribution becomes large (see figure 5).

Figure 8 shows the axial slip velocity $\bar{u}$, normalized by $\bar{u}_R^{SV}$, at the actual position $\hat{y} = \hat{y}^{SV}$ (corresponding to figure 7). It can be noted that the normalized value goes naturally to almost 1 as $\hat{R}$ becomes large for the largest value of $\hat{L}$. However this is not the case for the smallest value of $\hat{L}$ and the normalized value is reduced significantly for $\hat{R} > 50$.

Here again we explore the limit of validity of the classical approach by showing that both the position and the value of the slip velocity are dependent not only on $\hat{R}$ but also on $\hat{L}$. This is stressed also by results shown in figure 9 that gives a comparison of the slip velocity and the maximum normalized axial velocity on the guide axis, for $\hat{L} = 8000$ and $\hat{L} = 500$. We can see that the axial slip velocity departs from Rayleigh's solution more than that on the axis for small values of $\hat{R}$ (the solid curve is below
<figure><img src="image_132.png" /><figcaption>FIGURE 9. Normalized axial slip velocity versus $\hat{R}$ (solid line), and maximum normalized axial streaming velocity on the axis (dashed line). (a) $\hat{L} = 8000$, (b) $\hat{L} = 500$.</figcaption></figure>

the dashed curve). For large values of $\hat{R}$, the opposite is true: the slip velocity departs from Rayleigh's solution less than the axial streaming velocity on the guide axis. Comparing figure 9 (a,b) shows that this tendency is emphasized for smaller values of $\hat{L}$, due to increased importance of $\bar{u}_3$.

### 3.2.8. Average mass transport velocity

The streaming velocity considered up to now in this study, is often referred to as ‘Eulerian’ time-average velocity. The ‘average mass transport’ velocity, associated with streaming mass transport, is often considered, especially for thermoacoustic applications where heat transport due to streaming is the relevant quantity of interest. The axial component of the average mass transport velocity $\bar{u}^M$, at the first order of approximation, is given by

$$ \bar{u}^M(\hat{x}, \hat{y}) = \bar{u} + \frac{\overline{\rho' u'}}{\rho_0}. \qquad (3.42) $$

Along the axis, using (3.12) and neglecting terms of order $O(e^{-\hat{R}})$ and lower, it follows that

$$ \frac{\bar{u}^M(\hat{x}, 0)}{\bar{u}_R} = 1 - \frac{11}{2} \frac{1}{\hat{R}} - \frac{2\pi^2 \hat{R}^2}{45 \hat{L}^2} \left[ \hat{R} - 1 + \frac{30}{\hat{R}^2} - \frac{60}{\hat{R}^3} + \frac{45}{\hat{R}^4} \right]. \qquad (3.43) $$

Figure 10 shows the transverse profiles of axial streaming velocities $\bar{u}$ and $\bar{u}^M$, for three values of $\hat{R}$ (6, 50 and 90), and two values of $\hat{L}$ (8000 and 500). The difference between $\bar{u}$ and $\bar{u}^M$ is most significant in the viscous boundary layer, as usually considered. For small values of $\hat{R}$, differences are significant everywhere, including on the axis. These differences between $\bar{u}$ and $\bar{u}^M$ are obviously associated with the additional term $\overline{\rho' u'}$. As shown by expressions (3.4) and (3.6) this term is nearly zero since it is the product of two terms nearly in quadrature apart from viscous boundary layers effects. Indeed on the axis the acoustic velocity oscillates as $\cos(2\pi\hat{t})$ (see dominant terms in equation (3.4)) and the acoustic density oscillates as $\sin(2\pi\hat{t})$ (see dominant terms in equation (3.6)). This feature is associated with the standing wave
<figure><img src="image_133.png" /><figcaption>FIGURE 10. Axial velocities $\bar{u}^M(-1/4, \hat{y})/\bar{u}_R$ (dashed line), $\bar{u}(-1/4, \hat{y})/\bar{u}_R$ (solid line). (a) $\hat{R} = 6$, $\hat{L} = 8000$, (b) $\hat{R} = 50$, $\hat{L} = 8000$, (c) $\hat{R} = 90$, $\hat{L} = 8000$, (d) $\hat{R} = 6$, $\hat{L} = 500$, (e) $\hat{R} = 50$, $\hat{L} = 500$, (f) $\hat{R} = 90$, $\hat{L} = 500$.</figcaption></figure>

character of the oscillation in the guide under study that is closed at both ends and excited at its resonance frequency. In such resonators acoustic pressure (and thus density oscillations) are in quadrature with acoustic velocity due to boundary layers effects. Therefore for large guides (see figure 10b,c,e,f) the difference between $\bar{u}$ and $\bar{u}^M$ is only significant close to the wall. For narrow guides, the inner streaming takes up most of the guide, thus $\bar{u}$ and $\bar{u}^M$ are different through the whole guide section (see figure 10a,d).

Also in the case $\hat{R} = 90$ and $\hat{L} = 500$ (figure 10f), the change of sign of $\bar{u}$ at $\hat{y} \approx 25$ indicates the appearance of a new vortex with reverse flow near the axis.

### 3.2.9. Vertical streaming velocity

It can also be useful to consider the vertical component of streaming velocity, especially to detect the appearance of new cells. Let us define the average mass transport vertical velocity variable given by, in the first order of approximation

$$ \bar{v}^M(\hat{x}, \hat{y}) = \bar{v} + \frac{\bar{\rho}'\bar{v}'}{\bar{\rho}_0}. \qquad (3.44) $$

The continuity equation (first equation in (3.10)) becomes

$$ \frac{\partial \bar{u}^M}{\partial \hat{x}} + \hat{L} \frac{\partial \bar{v}^M}{\partial \hat{y}} = 0. \qquad (3.45) $$

Integrating equation (3.45) and using the boundary conditions (3.11c) yields the simplified expression

$$ \bar{v}^M(\hat{x}, \hat{y}) = \frac{\pi U_{ac}^2}{16c_0 \hat{L}} \cos(2\pi \hat{x}) \left\{ \frac{\hat{y}}{\hat{R}} \left( -33 + 6\hat{R} - 4\pi^2 \frac{\hat{R}^4}{\hat{L}^2} \right) + \frac{\hat{y}^3}{\hat{R}^3} \left( 11 - 6\hat{R} + \frac{8\pi^2 \hat{R}^4}{15 \hat{L}^2} \right) \right. \\ - \frac{4\pi^2 \hat{R}^4 \hat{y}^5}{15 \hat{L}^2 \hat{R}^5} + e^{\hat{y}-\hat{R}} \left[ \cos(\hat{R}-\hat{y}) \left( 16 + 4\frac{\hat{y}}{\hat{R}} \right) + \sin(\hat{R}-\hat{y}) \left( 16 - 4\frac{\hat{y}}{\hat{R}} \right) \right] \\ \left. + 2e^{2(\hat{y}-\hat{R})} \right\}. \qquad (3.46) $$
<figure><img src="image_134.png" /><figcaption>FIGURE 11. Vertical average mass transport velocities $\bar{v}^M(-1/2, \hat{y})/\bar{u}_R$ (dashed line), $\bar{v}(-1/2, \hat{y})/\bar{u}_R$ (solid line). (a) $\hat{R}=6$, $\hat{L}=8000$, (b) $\hat{R}=50$, $\hat{L}=8000$, (c) $\hat{R}=90$, $\hat{L}=8000$, (d) $\hat{R}=6$, $\hat{L}=500$, (e) $\hat{R}=50$, $\hat{L}=500$, (f) $\hat{R}=90$, $\hat{L}=500$.</figcaption></figure>

The Eulerian velocity $\bar{v}$ can be obtained using (3.44) and (3.13). Figure 11 shows the vertical velocities $\bar{v}$ and $\bar{v}^M$ for several values of $\hat{R}$ and two values of $\hat{L}$. For small values of $\hat{R}$ (figure 11a,d), $\bar{v}$ and $\bar{v}^M$ are very different. Also in the case $\hat{R}=90$ and $\hat{L}=500$ (figure 11f), the change of sign of $\bar{v}$ and $\bar{v}^M$ at $\hat{y}\approx 46$ indicates the appearance of a new vortex with reverse flow near the axis.

Finally, as a summary, figure 12 shows the streamlines of the streaming flow plotted using (3.42) and (3.46), for $\hat{R}=90$ and $\hat{L}=500$ and $8000$. For such large values of $\hat{R}$, the inner streaming is restrained in a very small region near the wall. For $\hat{L}=8000$, Rayleigh streaming is composed of two counter-rotating vortices as usual, while for $\hat{L}=500$ it is composed of four vortices, the flow along the axis being reversed.

## 4. The axisymmetric case

Streaming flow in cylindrical guides shows common features with the one between plates. It is composed of toroidal vortices that span $\lambda/4$ in the axial direction and inner and outer vortices coexist. Figure 1 corresponds to the streaming patterns observed in an axial section. However the value of the axial streaming velocity on the axis, and the position of the centre of the outer streaming vortex, are quite different from the case of parallel plates. We have performed an analysis similar to that presented above for the case of a cylindrical guide. One goal is to deduce the condition for $L/R$ such that additional counter-rotating streaming cells appear, associated with the third source term.

### 4.1. Acoustics

With the same change of variables as in the plane case ($\hat{r}=r/\delta_v$), where r is the radial coordinate), the linearized equations governing acoustics in the axisymmetric case are

$$ \frac{1}{2} c_0 \frac{\partial \rho'}{\partial \hat{t}} + \rho_0 \left( \frac{\partial u'}{\partial \hat{x}} + \hat{L} \frac{1}{\hat{r}} \frac{\partial}{\partial \hat{r}} (\hat{r} v') \right) = 0 \quad (4.1a) $$

$$ \rho_0 \frac{c_0}{2} \frac{\partial u'}{\partial \hat{t}} + c_0^2 \frac{\partial \rho'}{\partial \hat{x}} = \rho_0 \frac{\pi}{2} c_0 \frac{1}{\hat{r}} \frac{\partial}{\partial \hat{r}} \left( \hat{r} \frac{\partial u'}{\partial \hat{r}} \right) \quad (4.1b) $$

$$ \frac{\partial \rho'}{\partial \hat{r}} = 0, \quad (4.1c) $$
<figure><img src="image_135.png" /><figcaption>FIGURE 12. Streamlines of the mass transport velocity for $\hat{R} = 90$; $\hat{L} = 500$ (a), $\hat{L} = 8000$ (b). (c) Zoom on the inner region of top figure.</figcaption></figure>

where $u'$ and $v'$ are the axial and radial components of the acoustic velocity. The boundary conditions for $u'$ and $v'$ remain given by (3.3), replacing $\hat{y}$ by $\hat{r}$.

We start again with a given axial velocity solution of (4.1), approximated following Schuster & Matz (1940) for moderate and large channels ($\hat{R} \ge 6$)

$$u'(\hat{x}, \hat{r}, \hat{t}) = U_{ac} \cos(\pi\hat{x}) \left[ \left( -1 + \frac{1}{2} \frac{I_0((1-i)\hat{r})}{I_0((1-i)\hat{R})} + \frac{1}{2} \frac{I_0((1+i)\hat{r})}{I_0((1+i)\hat{R})} \right) \cos(\omega\hat{t}) \right. \\ \left. - \left( \frac{1}{2i} \frac{I_0((1-i)\hat{r})}{I_0((1-i)\hat{R})} - \frac{1}{2i} \frac{I_0((1+i)\hat{r})}{I_0((1+i)\hat{R})} \right) \sin(\omega\hat{t}) \right], \quad (4.2)$$

where $I_0$ is the modified Bessel function of the first kind.

Using the hypothesis that $\rho'$ depends on $\hat{x}$ and $\hat{t}$ only, together with the boundary conditions, equation (4.1a) is integrated with respect to $\hat{r}$ to obtain the radial acoustic
velocity component

$$v'(\hat{x}, \hat{r}, \hat{t}) = \frac{\pi}{4} U_{ac} \frac{\hat{r}}{\hat{L}} \sin(\pi \hat{x}) \left[ \frac{\hat{0} \tilde{F}_1 \left( ; 2, \frac{i \hat{r}^2}{2} \right) - \hat{0} \tilde{F}_1 \left( ; 2, \frac{i \hat{R}^2}{2} \right)}{I_0((1+i)\hat{R})} (\cos(\omega \hat{t}) - i \sin(\omega \hat{t})) \right. \\ \left. + \frac{\hat{0} \tilde{F}_1 \left( ; 2, \frac{-i \hat{r}^2}{2} \right) - \hat{0} \tilde{F}_1 \left( ; 2, \frac{-i \hat{R}^2}{2} \right)}{I_0((1-i)\hat{R})} (\cos(\omega \hat{t}) + i \sin(\omega \hat{t})) \right], \quad (4.3)$$

where $\hat{0} \tilde{F}_1(\cdot; b, z)$ is the regularized confluent hypergeometric function.

## 4.2. Streaming flow equations

Using the same change of variables as in the plane case, the equations to be solved are

$$\left\{
\begin{aligned}
\rho_0 \left( \frac{\partial \bar{u}}{\partial \hat{x}} + \hat{L} \frac{1}{\hat{r}} \frac{\partial}{\partial \hat{r}} (\hat{r} \bar{v}) \right) &= - \frac{\partial}{\partial \hat{x}} (\bar{\rho}' u') - \hat{L} \frac{1}{\hat{r}} \frac{\partial}{\partial \hat{r}} (\hat{r} \bar{\rho}' v') \\
\rho_0 \frac{\pi}{2} c_0 \frac{1}{\hat{r}} \frac{\partial}{\partial \hat{r}} \left( \hat{r} \frac{\partial \bar{u}}{\partial \hat{r}} \right) &= \frac{\partial \bar{p}}{\partial \hat{x}} + \rho_0 \frac{\partial}{\partial \hat{x}} (\bar{u}' u') + \rho_0 \hat{L} \frac{1}{\hat{r}} \frac{\partial}{\partial \hat{r}} (\hat{r} \bar{u}' v') \\
\rho_0 \frac{\pi}{2} c_0 \frac{1}{\hat{r}} \frac{\partial}{\partial \hat{r}} \left( \hat{r} \frac{\partial \bar{v}}{\partial \hat{r}} \right) &= \hat{L} \frac{\partial \bar{p}}{\partial \hat{r}} + \rho_0 \frac{\partial}{\partial \hat{x}} (\bar{u}' v') + \rho_0 \hat{L} \frac{1}{\hat{r}} \frac{\partial}{\partial \hat{r}} (\hat{r} \bar{v}' v')
\end{aligned}
\right.
\quad (4.4)$$

and the boundary conditions for $\bar{u}$ and $\bar{v}$ are now

$$\bar{u} = 0 \quad \text{at} \quad \hat{r} = \hat{R}; \quad \frac{\partial \bar{u}}{\partial \hat{r}} = 0 \quad \text{at} \quad \hat{r} = 0, \quad (4.5a,b)$$

$$\bar{v} = 0 \quad \text{at} \quad \hat{r} = 0 \quad \text{and} \quad \hat{r} = \hat{R}, \quad (4.5c)$$

$$\int_0^{\hat{R}} \hat{r} \bar{u} \, d\hat{r} = 0. \quad (4.5d)$$

The average products of fluctuations are calculated with the expressions of $u'$ and $v'$ from equations (4.2) and (4.3). Since the complete expressions are complicated, we only give here $\bar{u}'\bar{u}'$ and $\bar{u}'\bar{v}'$:

$$\bar{u}'\bar{u}'(\hat{x}, \hat{r}) = \frac{U_{ac}^2 \cos^2(\pi \hat{x})}{2} \frac{(I_0((1+i)\hat{r}) - I_0((1+i)\hat{R})) (I_0((1-i)\hat{r}) - I_0((1-i)\hat{R}))}{I_0((1-i)\hat{R}) I_0((1+i)\hat{R})} \quad (4.6)$$

$$\begin{aligned}
\bar{u}'\bar{v}'(\hat{x}, \hat{r}) = & \frac{1}{16}(1-i) \frac{\pi U_{ac}^2}{\hat{R}} \frac{1}{\hat{L}} \frac{\sin(2\pi\hat{x})}{I_0((1-i)\hat{R})I_0((1+i)\hat{R})} [ (I_0((1-i)\hat{r}) \\
& - I_0((1-i)\hat{R})) (\hat{R}I_1((1+i)\hat{r}) - \hat{r}I_1((1+i)\hat{R})) \\
& + (I_0((1+i)\hat{r}) - I_0((1+i)\hat{R})) (\hat{R}J_1((1+i)\hat{r}) - \hat{r}J_1((1+i)\hat{R})) ],
\end{aligned} \quad (4.7)$$

where $J_1$ is the Bessel function of the first kind.
<figure><img src="image_136.png" /><figcaption>FIGURE 13. (a) $T(\hat{r}, 20)$ (solid line), approximated function (dashed line). (b) Maximum relative error as a function of $\hat{R}$ (log-log plot).</figcaption></figure>

Taking the curl of the momentum equations results in

$$ -\frac{\pi}{2}c_0 \hat{L} \frac{\partial}{\partial \hat{r}} \left( \hat{r} \frac{\partial}{\partial \hat{r}} \left( \hat{r} \frac{\partial \bar{u}}{\partial \hat{r}} \right) \right) = -\hat{L} \frac{\partial^2(\overline{u'u'})}{\partial \hat{r} \partial \hat{x}} \\ -\hat{L}^2 \frac{\partial}{\partial \hat{r}} \left( \hat{r} \frac{\partial(\hat{r}u'v')}{\partial \hat{r}} \right) + \frac{\partial^2(\overline{u'v'})}{\partial \hat{x}^2} + \hat{L} \frac{1}{\hat{r}} \frac{\partial^2(\hat{r}v'v')}{\partial \hat{x} \partial \hat{r}}. \quad (4.8) $$

### 4.3. Results

In the same way as in the plane case, we consider separately the four problems associated with each source term in equation (4.8).

However, there is no symbolic form for the integral of terms of the form $T(\hat{r}, \hat{R}) = \text{Re}[(1 + i)I_0((1 + i)\hat{r})/I_0((1 + i)\hat{R}) \cdot I_1((1 - i)\hat{r})/I_0((1 - i)\hat{R})]$ that appear in the successive integrations. Nevertheless these terms can be approximated with good accuracy using an exponential function that can be integrated symbolically. Figure 13 (a) shows $T(\hat{r}, 20)$ together with the approximated function (a purely exponential function based on a two-point interpolation) that is used instead in the successive integrations. The curves are nearly superimposed and indistinguishable. The maximum relative error between the two functions is represented in a log-log plot in figure 13 (b) for values of $\hat{R}$ between 6 and 200, showing excellent accuracy and a variation in $1/\hat{R}^2$. The full expressions for the streaming velocity components are not developed below, since the formulae are very lengthy. Instead, the following fitted formula is obtained for the axial streaming velocity along the axis, normalized by Rayleigh's solution, which is now $\bar{u}_R = -(3/8)U_{ac}^2/c_0 \sin(2\pi\hat{x})$

$$ \frac{\bar{u}}{\bar{u}_R} = 1 - 5.676 \frac{1}{\hat{R}^{0.98}} - 0.296 \frac{\hat{R}^{2.98}}{\hat{L}^2}. \quad (4.9) $$

This expression shows a similar but slightly different dependence on $\hat{R}$ compared to that obtained in the plane case (3.39). It also shows that an extra cell of reversed flow may still appear near the axis in circumstances where the ratio $\hat{R}^{2.98}/\hat{L}^2$ in the third term is not negligible. More precisely, the condition for this is

$$ \hat{L} \le \hat{L}_{\text{limit}}(\hat{R}) = \left( \frac{0.296 \hat{R}^{2.98}}{1 - \frac{5.676}{\hat{R}^{0.98}}} \right)^{1/2}. \quad (4.10) $$
<figure><img src="image_137.png" /><figcaption>FIGURE 14. Limit value $\hat{L}_{limit}(\hat{R})/\hat{R}$ under which the flow is reversed on the axis.
Minimum values: $\hat{R} = 11.93$, $\hat{L}/\hat{R} = 2.6$.</figcaption></figure>

<figure><img src="image_138.png" /><figcaption>FIGURE 15. Axial velocities $\bar{u}(-1/4, \hat{y})/\bar{u}_R$, $\hat{L} = 500$ (solid line), $\hat{L} = 8000$ (dashed line).
(a) $\hat{R} = 6$, (b) $\hat{R} = 50$, (c) $\hat{R} = 100$.</figcaption></figure>

<figure><img src="image_139.png" /><figcaption>FIGURE 16. Radial velocities $\bar{v}(-1/2, \hat{y})/\bar{u}_R$, $\hat{L} = 500$ (solid line), $\hat{L} = 8000$ (dashed line).
(a) $\hat{R} = 6$, (b) $\hat{R} = 50$, (c) $\hat{R} = 100$.</figcaption></figure>

Figure 14 shows, for the axisymmetric case, the variation with $\hat{R}$ of $\hat{L}_{limit}/\hat{R}$ for the appearance of reverse flow on the axis. For example, for $\hat{R} = 100$, the extra cell will appear for $\hat{L} < 536$ in the axisymmetric case, that is for a shorter guide than in the plane case. Note that for small values on the vertical axis of figure 14 the limit of validity of the hypothesis $(R/L)^2 \ll 1$ is reached.

This has to be kept in mind when assessing geometries compatible with the existence of the extra cell.

Figures 15 and 16 show the radial profiles of axial and radial streaming velocities $\bar{u}$, $\bar{v}$, for three values of $\hat{R}$ (6, 50 and 100), and two values of $\hat{L}$ (500 and 8000). The behaviour of the streaming flow is similar to that in the plane case. The inner
<figure><img src="image_140.png" /><figcaption>FIGURE 17. Streamlines of the mass transport velocity in an axial plane section for $\hat{R} = 100$; $\hat{L} = 500$ (a), $\hat{L} = 8000$ (b).</figcaption></figure>

streaming is of lower amplitude than in the plane case. It is also seen that the centre of the outer vortex is closer to the wall than in the plane case, as expected (Rayleigh 1884). The influence of the third source term is visible in the case $\hat{R} = 50$, $\hat{L} = 500$. An extra vortex appears and the flow is reversed on the axis for $\hat{R} = 100$, $\hat{L} = 500$. Figure 17 shows the corresponding streamlines of the mass transport velocity in an axial plane section. The inner vortex is very narrow in this case and not visible on this figure.

## 5. Conclusion

Rayleigh's outer streaming flow in standing wave guides (infinite parallel plates or cylindrical tubes) was revisited, by solving the full linear streaming equations for moderate to large guides, using symbolic computations. The effect of each Reynolds stress source term was analysed separately and several features that were not previously reported were exhibited: analysis of the effect of the two classically dominant Reynolds stress source terms highlighted the commonly accepted fact that the streaming flow is generated only by viscous boundary layer effects. The outer streaming is created from momentum diffusion of the velocity created near the guide wall, and not from pressure gradient effects. The dependency with the guide's length and the radius of the axial streaming velocity normalized by Rayleigh's solution along the axis was quantified. It was shown that the effect of the third Reynolds stress source term which is usually neglected can become important for some geometrical parameters that were identified. In these situations, a vortex with reverse flow is created along the guide axis. Several characteristics of the streaming flow were analysed, such as slip velocity, mass transport velocity and transverse
component of streaming flow, in connection with the effect of the third source term. For all these characteristics it was shown that not only the ratio of wave guide width (or radius) to viscous boundary layer thickness, but also the ratio of guide width (or radius) to guide length, are driving criteria for the organization of streaming.

REFERENCES

AMARI, M., GUSEV, V. & JOLY, N. 2003 Temporal dynamics of the sound wind in acoustitron. *Acta Acust. united with Acustica* **89**, 1008–1024.

BAILLIET, H., GUSEV, V., RASPET, R. & HILLER, R. A. 2001 Acoustic streaming in closed thermoacoustic devices. *J. Acoust. Soc. Am.* **110**, 1808–1821.

BOLURIAAN, S. & MORRIS, P. J. 2003 Acoustic streaming: from Rayleigh to today. *Intl J. Aeroacoust.* **2** (3–4), 255–292.

DARU, V., BALTEAN-CARLÈS, D., WEISMAN, C., DEBESSE, P. & GANDIKOTA, G. 2013 Two-dimensional numerical simulations of nonlinear acoustic streaming in standing waves. *Wave Motion* **50**, 955–963.

DARU, V., REYT, I., BAILLIET, H., WEISMAN, C. & BALTEAN-CARLÈS, D. 2017a Acoustic and streaming velocity components in a resonant wave guide at high acoustic levels. *J. Acoust. Soc. Am.* **141** (1), 563–574.

DARU, V., WEISMAN, C., BALTEAN-CARLÈS, D., REYT, I. & BAILLIET, H. 2017b Inertial effects on acoustic Rayleigh streaming flow: transient and established regimes. *Wave Motion* **74**, 1–17.

FAVRE, A. 1965 Équations des gaz turbulents compressibles ii.-méthode des vitesses moyennes; méthode des vitesses macroscopiques pondérées par la masse volumique. *J. Méch.* **87**, 391–421.

HAMILTON, M. F., ILINSKII, Y. A. & ZABOLOTSKAYA, E. A. 2003a Acoustic streaming generated by standing waves in two-dimensional channels of arbitrary width. *J. Acoust. Soc. Am.* **113** (1), 153–160.

HAMILTON, M. F., ILINSKII, Y. A. & ZABOLOTSKAYA, E. A. 2003b Thermal effects on acoustic streaming in standing waves. *J. Acoust. Soc. Am.* **114** (6), 3092–3101.

LIGHTHILL, J. 1978 Acoustic streaming. *J. Sound Vib.* **61**, 391–418.

MENGUY, L. & GILBERT, J. 2000a Non-linear acoustic streaming accompanying a plane stationary wave in a guide. *Acta Acust. united with Acustica* **86**, 249–259.

MENGUY, L. & GILBERT, J. 2000b Weakly nonlinear gas oscillations in air-filled tubes; solutions and experiments. *Acta Acust. united with Acustica* **86** (5), 798–810.

MOREAU, S., BAILLIET, H. & VALIÈRE, J.-C. 2008 Measurements of inner and outer streaming vortices in a standing waveguide using laser doppler velocimetry. *J. Acoust. Soc. Am.* **123** (2), 640–647.

NYBORG, W. L. 1953 Acoustic streaming due to attenuated plane waves. *J. Acoust. Soc. Am.* **25**, 68–75.

NYBORG, W. L. 1965 Acoustic streaming. In *Physical Acoustics* (ed. W. P. Mason), vol. 2B, pp. 265–331. Academic Press.

OLSON, J. R. & SWIFT, G. W. 1997 Acoustic streaming in pulse tube refrigerators: tapered pulse tubes. *Cryogenics* **37**, 769–776.

PARIDAENS, R., KOUIDRI, S. & JEBALI JERBI, F. 2013 Investigation on the generation mechanisms of acoustic streaming in a thermoacoustic prime mover. *Cryogenics* **58**, 78–84.

PIERCE, A. D. 1989 *Acoustics: An Introduction to its Physical Principles and Applications*. ASA.

QI, Q. 1993 The effect of compressibility on acoustic streaming near a rigid boundary for a plane traveling wave. *J. Acoust. Soc. Am.* **94**, 1090–1098.

RAYLEIGH, LORD 1884 On the circulation of air observed in Kundt tubes, and on some allied acoustical problems. *Phil. Trans. R. Soc. Lond.* **175**, 1–21.

REYT, I., DARU, V., BAILLIET, H., MOREAU, S., VALIÈRE, J.-C., BALTEAN-CARLÈS, D. & WEISMAN, C. 2013 Fast acoustic streaming in standing waves: generation of an additional outer streaming cell. *J. Acoust. Soc. Am.* **134**, 1791–1801.
REYT, I., BAILLIET, H. & VALIÈRE, J.-C. 2014 Experimental investigation of acoustic streaming in a cylindrical wave guide up to high streaming Reynolds numbers. *J. Acoust. Soc. Am.* **135** (1), 27–37.

ROTT, N. 1974 The influence of heat conduction on acoustic streaming. *Z. Angew. Math. Phys.* **25**, 417–421.

SCHLICHTING, H. 1932 Berechnung ebener periodischer Grenzschicht- strommungen [calculation of plane periodic boundary layer streaming]. *Phys. Zcit.* **33**, 327–335.

SCHUSTER, V. K. & MATZ, W. 1940 Über stationare Strömungen im Kundtschen Rohr (on stationary streaming in Kundt tubes). *Akust. Zh.* **5**, 349–352.

SUGIMOTO, N. 2016 Nonlinear theory for thermoacoustic waves in a narrow channel and pore subject to a temperature gradient. *J. Fluid Mech.* **797**, 765–801.

THOMPSON, M. W., ATCHLEY, A. A. & MACCARONE, M. J. 2004 Influences of a temperature gradient and fluid inertia on acoustic streaming in a standing wave. *J. Acoust. Soc. Am.* **117** (4), 1839–1849.

WAXLER, R. 2001 Stationary velocity and pressure gradients in a thermoacoustic stack. *J. Acoust. Soc. Am.* **109** (6), 2739–2750.

WESTERVELT, P. J. 1953 The theory of steady rotational flow generated by a sound field. *J. Acoust. Soc. Am.* **25**, 60–67.

WIKLUND, M. G., GREEN, R. & OHLIN, M. 2012 Acoustofluidics 14: applications of acoustic streaming in microfluidic devices. *Lab on a Chip* **12**, 2438–2451.

WOLFRAM RESEARCH INC. 2018 Mathematica, Version 11.3, Champaign, IL.

ZAREMBO, L. K. 1971 Acoustic streaming. In *High-Intensity Ultrasonic Fields* (ed. L. D. Rozenberg), vol. 175, pp. 135–199. Plenum Press.
O. Hireche, I. Ramadan, C. Weisman, H. Bailliet, Y. Fraigneau, D. Baltean-Carlès, V. Daru

Experimental and numerical investigation of natural convection flows in two horizontal thermoacoustic cavities.

Int. J. Heat Mass Transf. 149 (2020), 119195.
# Experimental and numerical investigation of natural convection flows in two horizontal thermoacoustic cavities

<figure><img src="image_141.png" /></figure>

Check for
updates

Omar Hireche<sup>a</sup>, Islam Ramadan<sup>b</sup>, Catherine Weisman<sup>a,c,*</sup>, Hélène Bailliet<sup>b</sup>, Yann Fraigneau<sup>a</sup>, Diana Baltean-Carlès<sup>a,c</sup>, Virginie Daru<sup>a,d</sup>

## ARTICLE INFO

**Article history:**

Received 1 August 2019

Revised 22 November 2019

Accepted 8 December 2019

**Keywords:**

Natural convection

Thermoacoustics

Porous medium

Particle image velocimetry

Numerical simulation

## ABSTRACT

This study aims at exploring the potential generation of natural convective flows in thermoacoustic devices. To this effect, two thermoacoustic setups are considered without acoustic oscillations, each setup thus constituting a cavity filled with air and comprising a porous medium heated on one of its ends. In the first configuration, the porous medium is a ceramic stack. In the second configuration, an ambient heat exchanger using circulating ambient water is placed on the other end of the porous medium (which is a regenerator). The wall temperature is measured in several points using thermocouples and natural convection flows developing in these cavities are measured in a vertical plane using Particle Image Velocimetry. Moreover numerical calculations are performed based on simplified 3D configurations, using immersed solid bodies at prescribed temperature to represent the heating or cooling elements and describing the porous media as homogeneous and either isotropic or anisotropic. Navier-Stokes equations under Boussinesq assumption are solved in the entire domain with the addition of a Darcy correction term in the porous media. Excellent qualitative and quantitative agreement on velocity fields is obtained between experimental and numerical simulation results for both configurations, which validates the model description. It is shown that the observed natural convection flows are fully three-dimensional and of significant magnitude so that they could potentially interfere with thermoacoustic effects.

© 2019 Elsevier Ltd. All rights reserved.

## 1. Introduction

Thermoacoustic devices use thermal and dynamic interactions between the acoustic oscillations of a working inert gas and conducting solid walls to generate either heat pumping for refrigerating purposes or mechanical work for engine operating purposes [1]. These interactions are called the thermoacoustic effect. They occur nearby the solid walls, in a thermal boundary layer area: the gas experiences acoustic oscillations, while exchanging heat with the solid wall. In order to increase the exchange surface area, a porous medium is used. This porous medium is placed between two heat exchangers, forming the thermoacoustic core (TAC). The heat exchangers enable to either impose a temperature gradient along the porous medium, in order to produce acoustic power (engine), or to benefit from the heat pumping effect steered by an imposed acoustic oscillation (refrigerator). The thermoacoustic core

represents the active part of the device, where the energy conversion takes place. Usually, thermoacoustic devices also include an acoustic resonator for sustaining the acoustic wave, and potentially an acoustic source or a linear electric generator, depending on the desired conversion effect. Thermoacoustic devices represent an interesting ecological alternative (almost no moving parts, the use of inert gas, simplicity) for industrial applications: cryogenics [2], air conditioning [3], domestic refrigeration, solar powered thermoacoustic engines [4], electric energy generators (conversion of the acoustic power into electric power) [5], waste heat recovery (for example, thermoacoustic heat pumps are used to upgrade the waste heat for renewed use in the industrial process) [6]. However thermoacoustic devices still suffer from poor performance generally caused by the existence of several nonlinear phenomena interfering with the process of energy transfer and deteriorating the global performance of the device. These nonlinear phenomena are related to jet-driven streaming, acoustic streaming, harmonic generation, end-effects and natural convection, and are not well described or just partially. Among these phenomena, natural con-

<footnote><sup>a</sup>Laboratoire LIMSI, CNRS, Campus universitaire bât. 507, Rue du Belvedère, Orsay Cedex F-91405, France</footnote>

<footnote><sup>b</sup>Institut Pprime, CNRS - Université de Poitiers - ENSMA, ENSIP, 6 rue Marcel Doré, Bât. B17 - BP 633, Poitiers Cedex 86022, France</footnote>

<footnote><sup>c</sup>Sorbonne Université, Faculté des Sciences et Ingénierie, 4 Place Jussieu, Paris 75005, France</footnote>

<footnote><sup>d</sup>DynFluid Lab., ENSAM, 151 Bd de l'Hôpital, Paris 75013, France</footnote>

<footnote>* Corresponding author.
E-mail address: weisman@limsi.fr (C. Weisman).</footnote>
vection has been the least investigated. The present study focuses on the effect of natural convection flows.

The influence of natural convection on starting conditions for a thermoacoustic engine was pointed out in an early paper of Gardner and Swift [7]. They observed that in the engine under study, the acoustic oscillation did not start although the design parameters were respected. This appeared to be due to a distor- tion of the mean temperature field, which was identified by the authors as originating from gravity-driven convection in the stack. The experimental procedure was modified, in order to reduce the convection flow, by decreasing the mean pressure and start the oscillation at a lower pressure. Once the oscillation was stabilized, the mean pressure was changed to its designed value. This is the option generally chosen by thermoacousticians: try to avoid or to reduce natural convection; but most of the time it is impossible to completely avoid it, since the devices come in a variety of geome- tries. The influence of the tilt angle on the starting condition of a solar powered thermoacoustic engine was analyzed more recently using Thermal Infrared Imager temperature measurements in the TAC [8]. In cousin devices that are pulse tube refrigerators, natural convection was also pointed out as a major source of losses [9,10] but no detailed analysis was performed and no estimate or measurement of the natural convection velocity field is available.

Apart from these few experimental evidence of natural con- vective flows in real devices, some experimental studies aimed at characterizing natural convection in simplified academic configu- rations. Natural convection and heat transfer inside a horizontal empty resonator, taking place between a hot heat exchanger and a cold heat exchanger, was investigated experimentally using ther- mal infrared imaging and particle image velocimetry (PIV) [11], and for different tilt angles [12]. It was shown that heat transfer mode changed from natural convection to forced convection with the addition of acoustic oscillation.

The role of natural convection on the distortion of Rayleigh streaming flow was exhibited in some recent experimental studies [13–15]. Measurement of the sole convection velocity field along a vertical section of the resonator under a temperature gradient was performed in [14] in order to explain the distorted streaming velocity field observed, but the measurement was very localized and there was no theoretical explanation of the complex velocity measurements obtained. Measurements of the mean velocity field in the entire cold resonator part of a simple standing-wave thermoacoustic engine were recently performed by Ramadan et al. [15]. Using both PIV and LDV techniques, they showed the existence of a mean flow with or without thermoacoustic oscil- lations. This flow was identified as natural convection combined, when acoustic waves are present, with Rayleigh streaming. How- ever, in this preliminary study, measurements were performed only in the central horizontal plane, which is limiting for the study of natural convection. More generally, a characterization of the 3D nature of the natural convection velocity field suggested by all above-cited studies is still lacking and its influence on the overall thermoacoustic heat transfer not yet quantified. This is because dedicated adequate measurements of natural convection flow and temperature distribution are difficult in real devices, due to the many coupled intricate physical phenomena.

From the numerical point of view, natural convection is not yet included in thermoacoustic designing tools. It is therefore necessary to develop a numerical tool that can predict the mag- nitude and flow pattern of natural convection flows as well as the associated temperature distribution. So far it is unrealistic to integrate all effects in one complete numerical code, for a given complex device: acoustics, thermoacoustic energy conversion, acoustic streaming, natural convection without powerful comput- ing resources. This is due to the existence of several distinct time scales (acoustic, thermal diffusion) and spatial scales (wavelength,

viscous and thermal boundary layer thickness, device's charac- teristic scales). Such simulations would therefore be very time consuming, since long physical times would have to be discretized in small time steps, and extended physical domains discretized with fine mesh. Morever, for a better understanding of how sev- eral physical processes interact, it seems relevant to study them separately in a first approach. Therefore, this study focuses on the natural convection flow in thermoacoustic device configurations, with no acoustic oscillation. As will be shown further, the geo- metrical configuration of thermoacoustic devices can generally be considered as a series of heated cavities, some of which are filled with a porous medium. Each cavity is submitted to a temperature difference, which can generate a heat-driven flow.

Convection flows through non homogenous and anisotropic porous media in cavities heated from the side are quite well de- scribed, with studies started in the 1980s [16]. The same problem in mixed fluid and porous cavities was treated extensively in the literature [17,18]. The convection flow depends strongly on the geometry of the cavities and the previous general studies cannot be directly applied to the specific case of a thermoacoustic engine. The specific geometry of a thermoacoustic stack was considered in [19], where natural convection in a 2D differentially heated horizontal cavity, filled with a stack of plates, was considered. This study showed that a Darcean porous model gives correct orders of velocity scales, and allows for a qualitative description of the solutions. A preliminary 3D numerical study of the effects of natural convection in a cylindrical guide containing an assembly of a hot heat exchanger and porous medium was conducted recently [20]. The influence of the porous medium's physical characteristics (permeability, thermal conductivity, anisotropy) on the flow and temperature fields was investigated.

The numerical approach of [20] is pursued in this paper on realistic configurations, and confronted with experiments resulting from extending the preliminary study in [15]. In the present anal- ysis, the experimental and numerical models were elaborated and calibrated jointly in order to improve the description of natural convection flows in a thermoacoustic system. The first objective of this paper is to predict natural convection flows in configurations that are representative of such systems. The second objective is to show that the model works well by comparison with possible measurements in parts of the device.

As mentioned above, to isolate natural convection from other phenomena, we chose to work with a thermoacoustic cavity but not with a thermoacoustic device; i.e. no acoustic oscillation was generated in the cavity. Whether or not natural convection should be accounted for in the design of a thermoacoustic device is going to depend highly on the device's geometry and its orientation with respect to gravity. It is also going to depend particularly on the thermo-physical parameters of all constituting elements, and specifically the thermoacoustic core. Therefore, and in order to explore increasingly realistic geometries, we investigate two different configurations: the first one consists in a very simplified thermoacoustic core that is composed of an electric heater at- tached to a ceramic stack; the second one is composed of the same two elements followed by an ambient heat exchanger. These two configurations are further referred to as horizontal cavity without ambient heat exchanger (or 1st configuration), and horizontal cavity with ambient heat exchanger (or 2nd configuration).

In the following, the first configuration is investigated in Section 2. The experimental setup and method are described (Section 2.1) and the numerical method is presented (Section 2.2). Experimental and numerical results are then compared and dis- cussed, and the potentiality of the numerical model for predicting purposes is assessed (Section 2.3). The configuration with ambient heat exchanger is considered in Section 3, following the same approach.
**2. Horizontal cavity without ambient heat exchanger (1st configuration)**

**2.1. Experimental setup, measurement systems and procedure**

For this first comparison between experimental and numerical approaches, a very simple thermoacoustic cavity is used. As shown in Fig. 1a, the resonator is made of a glass tube with inner and outer diameters of 39 mm and 44 mm, respectively. The total length of the resonator is 738 mm and the thermoacoustic core is placed at a distance of 129 mm away from one end of the resonator. The thermoacoustic core consists of an electric heater attached to one side of a 400 CPSI (Cell Per Square Inch) ceramic stack that has length of 60 mm. The axial coordinate is also defined in Fig. 1a, with $x = 0$ corresponding to the right end of the ceramic stack. In such a simple thermoacoustic engine, applying high enough temperature gradient along the stack generates sustainable standing acoustic waves. A more detailed description for the experimental setup can be found in [15]. In order to isolate the natural convection flow, the two small plugs at the ends of the resonator tube were removed to cease any acoustic oscillations and velocity measurements were performed in the central vertical plane. It was checked previously [15] that the removal of the plugs does not affect the natural convection flow.

A PIV system is used to measure the natural convection flow in the central vertical plane. The PIV system consists of a double pulse Nd:YAG laser source that illuminates the vertical plane. Two cameras (Imager Pro X) with spatial resolution of 2048 × 2048 pixels² and maximum repetition rate of 15 Hz are used to record the images. Each camera is equipped with a Nikon lens (Model: Sigma 105 mm). The two cameras are placed next to each other in order to increase the size of the measurement Field Of View (FOV). An atomizer aerosol generator (Model: TOPAS ATM 210) is used to atomize oil to produce seeding particles with a mean particle diameter of 0.4 µm. In order to achieve better correlation between the acquired double-images, the time between laser pulses is adjusted to allow the seeding particles to displace a distance of about 1/4 the size of the interrogation area [21].

The measurement procedure followed in the current study is similar to the procedure we followed in our previous work [15]. In order to have enough seeding particles after the system reaches to the steady-state conditions (i.e. the change in the temperature over the PIV measurement period does not exceed 5%), the preheating technique is followed. A 31.5 W ± 1% electrical power is applied to the hot heat exchanger (coiled resistance). After the preheating of the system, the seeding particles are introduced from one

extremity and flow until reaching the other extremity of the duct so that the seeding particles are dispersed over the whole duct. The system is left for 3 min in order to eliminate the effects of the seeding introduction on the natural convection flow. Finally, the PIV measurements are started and 100 double-images captured. The natural convection flow distribution is measured in the central vertical plane, covering the axial distance from $x = 0$ mm (end of ceramic stack) to $x = 520$ mm. Since wall temperatures are very important for the study of the natural convection flow, several thermocouples are used to measure the wall temperatures and the temperatures at the center of the two side sections of the ceramic stack. As shown by Fig. 1b (large dots), on the tube wall for each axial position there are two temperature measurements corresponding to the top and bottom of the tube wall.

The uncertainty of the measured quantities (i.e. temperature and the velocity) is estimated by repeating the measurements (3 times) under the same input electric power to the electrical heater and calculating the standard deviation of these quantities. The estimated value of the uncertainty in the temperature measurements is thus found to be ± 1 % of the reported mean temperatures. For the velocity measurements, the uncertainty of the measurements is calculated over the whole measurement area and then the spatial-average uncertainty is calculated. The spatial-average uncertainty is found to be ± 13 % of the reported velocities.

**2.2. Numerical model**

**2.2.1. Numerical simulation parameters**

A simplified 3D configuration is introduced in order to model the experimental configuration described in the previous section (Fig. 2a). The domain under study is a horizontal cylindrical tube of length $L$ and radius $R$, filled with air and closed at both ends.

The ceramic stack, constituted of horizontal capillaries of square cross section with side length $a = 1$ mm, is modeled as a homogeneous and anisotropic porous medium of length $L_p$ and radius $R$ (see Fig. 2a and Table 1a). Its porosity $\epsilon$, defined as the ratio of gas volume to total volume is calculated based on the given geometry. The equivalent permeability in the axial direction is $K_{xx} = \frac{a^2\epsilon}{8}$ (calculated using [18]), its permeability in the radial and azimuthal directions are approximately zero (in the code they are set to $10^{-12}$ m²). The equivalent thermal conductivity is assumed isotropic and homogeneous, of scalar value $k_p = \epsilon k_{fluid} + (1 - \epsilon) k_{solid}$ (again based on [18]). The numerical values of these parameters are given by Table 1b.

The hot exchanger (coiled resistance), placed at a distance $L_L$ from the left end, is modeled as 2 solid rings of width $e_h$

<figure><img src="image_142.png" /><figcaption>**Fig. 1.** (a): Schematic for the thermoacoustic resonator and image of the thermoacoustic core. (b): Thermocouples locations and results of measurements (in $K$) in the 1st configuration.</figcaption></figure>
<figure><img src="image_143.png" /><figcaption>**Fig. 2.** (a): Numerical configuration for the 1st configuration, (b): temperature prescribed in the numerical model along the tube wall, as a function of axial position (solid line), and measurements on top/bottom of the outside tube wall (squares) reported from Fig. 1b.</figcaption></figure>

<caption>**Table 1**</caption>

(a): Geometrical and thermal parameters, (b): Thermophysical parameters for air and porous medium, 1st configuration.

<table><thead><tr><th>R</th><th>L</th><th>L<sub>L</sub></th><th>e<sub>h</sub></th><th>L<sub>p</sub></th><th>T<sub>hot</sub></th><th>T<sub>cold</sub></th></tr></thead><tbody><tr><td>mm</td><td>mm</td><td>mm</td><td>mm</td><td>mm</td><td>K</td><td>K</td></tr><tr><td>39</td><td>738</td><td>129</td><td>5.8<br>(a)</td><td>60</td><td>520</td><td>293</td></tr><tr><td>ρ</td><td>μ</td><td>k</td><td>c<sub>p</sub></td><td>ε</td><td>K<sub>xx</sub></td><td>k<sub>p</sub></td></tr><tr><td>kg m<sup>-3</sup></td><td>kg m<sup>-1</sup> s<sup>-1</sup></td><td>W m<sup>-1</sup> K<sup>-1</sup></td><td>J kg<sup>-1</sup> K<sup>-1</sup></td><td></td><td>m<sup>2</sup></td><td>W m<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>1.07</td><td>1.9 × 10<sup>-5</sup></td><td>2.8 × 10<sup>-2</sup></td><td>1004<br>(b)</td><td>0.64</td><td>8.10<sup>-8</sup></td><td>9.18.10<sup>-1</sup></td></tr></tbody></table>

and a central core maintained at certain hot temperature (the central core being hotter than the outside ring). Temperature values prescribed in this hot heat exchanger are: 420 K on the outside ring, 523 K on the middle ring, and 627 K in the central core. Therefore the “average” hot temperature is estimated as $T_{hot} = 520$ K. Along the tube wall, temperature is imposed as a function of the axial coordinate only, based on the measured temperature distribution, as shown by Fig. 2b. Since for each axial position there are two measures of temperature (at the top and bottom of the outside wall) reported in Fig. 1b, the average value between the two measured values was calculated and prescribed for the corresponding axial position. On the left of the HHX, the temperature distribution is then fit with an exponential function. Along the porous medium a linear variation is used. Finally on the right of the porous medium the temperature distribution is fit with another exponential function (Fig. 2b).

The thermophysical parameters for air are calculated for atmospheric pressure and at a reference temperature chosen to be $T_{ref} = 320$ K which is the global mean temperature of the air.

### 2.2.2. Governing equations

In the fluid domain, the Navier-Stokes equations are written as classically under Boussinesq assumptions. The reference air density $\rho$, thermal conductivity $k$, heat capacity $c_p$, dynamic viscosity $\mu$, and volume expansion coefficient at constant pressure $\beta = 1/T_{ref}$ for an ideal gas are constant and equal to their value for a reference temperature $T_{ref}$ and atmospheric pressure.

The thermophysical characteristics of the porous medium are its porosity $\epsilon$, its permeability tensor $K$, and its equivalent thermal conductivity $k_p$. A scalar phase variable $X$ is introduced to identify the nature of the medium ($X = 1$ for gas, $X = 0$ for the porous medium) so that a Darcy term appears as additional term in the
right hand side of the momentum equation in the case of the porous medium zone. Therefore the equations can be written in a unified form as:

$$ \nabla \cdot \mathbf{V} = 0 $$

$$ \frac{\partial(\rho\mathbf{V})}{\partial t} + \nabla \cdot (\rho\mathbf{V} \otimes \mathbf{V}) = -\nabla p + \mu \Delta\mathbf{V} - \rho\beta(T - T_{ref}) $$

$$ \mathbf{g} - (1 - X)\mu\mathbf{K}^{-1}\mathbf{V} $$

$$ \rho c_p \left( \frac{\partial T}{\partial t} + (\mathbf{V} \cdot \nabla) T \right) = \nabla \cdot \left[ (X k + (1 - X) k_p) \nabla T \right] $$

In the above equations **V** is the velocity vector, **p** is the nonhydrostatic pressure, **T** is the temperature, **t** is time and **g** is the acceleration of gravity.

On all solid boundaries (i.e. the cylinder walls and the heat exchanger solid parts), a no-slip condition is imposed, and the temperature is prescribed as detailed in the previous section.

### 2.2.3. Numerical method

Computations are performed with the in-house code SUN-FLUIDH, which has been developed at LIMSI for simulations of unsteady incompressible flows. The previous governing equations are solved in the entire domain using a finite volume approach with a second order discretization in space and time. An implicit treatment of the viscous and thermal conductive terms ensures better numerical stability with regard to the time step which only depends on the CFL criterion. The equation set is solved by means of an ADI method [22]. An incremental projection method is used in order to ensure the divergence-free velocity field [23,24]. This implies to solve a Poisson's equation at each time step which is carried out with a multigrid-SOR algorithm [25]. 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 regular.

The numerical domain is meshed using 64 points regularly spaced along the radial direction, 64 points regularly spaced along the azimuthal direction, and 512 points placed regularly along the axis. The mesh size in the axial direction is $\Delta x = 1.4$ mm. In the radial direction, $\Delta r = 0.61$ mm. It was checked through a mesh-refinement study that this spatial resolution is adequate for the given configuration. The time step is adjusted throughout the simulation in order to satisfy a CFL < 0.5 criterion.

## 2.3. Results and discussion

Fig. 3 shows the experimental and numerical velocity vector maps and colored contours of the axial velocity component for the central vertical plane at the right side of the porous medium, which corresponds to the measurement field of view (close-up views for the figures near the porous medium (i.e. $x = 0$) are also shown in Fig. 3). Fig. 4 shows the numerical temperature field over the same region. The vortex flow observed is indeed that of natural convection of air in a cylindrical fluid cavity, triggered by a horizontal temperature gradient (the temperature of the right side of the porous medium is still much higher than the ambient temperature, as shown on Fig. 4). The numerical temperature at the center of the right end section of the porous medium is found to be equal to 320 K, in good agreement with the measured value of 328 K (see Fig. 1b). Experimental and numerical results agree remarkably well (Fig. 3), considering the numerous simplifications and approximations made in the numerical model. It should be pointed out that the numerical calculation converges to a steady solution, in agreement with experimental results (in which the temperature and thus velocity fields only change very slowly in time). The main difference lies in that for $x > 150$ mm a small flow remains observable along the whole experimental field of

view, while in the numerical simulation it becomes negligible. This could be due to the fact that the temperature imposed on the tube boundary is based on the measurement of the outside temperature and not on the inside, which should be slightly higher, therefore extending the natural convection flow. Also the presence of the openings on the tube ends (open plugs) is not modeled (a closed end with no slip condition is imposed in the numerical domain). In addition, the maximum amplitude and the extent of the convection flow in the field of view are highly dependent on the value of the driving $\Delta T$ imposed through prescribed boundary temperature values. The temperature on the glass tube is not controlled and the ambient temperature varied slightly during the experiment (about 2 degrees). The effect of changing $\Delta T$ by $\pm 2$ K (i.e. representing 12% of $\Delta T = 17$ K for the experimental field of view) could also explain the observed discrepancy between the measurements and the numerical results: A rough model gives the max velocity amplitude as a linear function of $\Delta T$. The corresponding variation of the extent of the convection cell is harder to estimate. Fig. 5 shows the axial profile of the vertical velocity component obtained from experiments and numerical simulations at the central vertical plane, again on the right side of the porous medium, for several values of the radial coordinate $r$ (which on this plane measures the height). Here again the agreement between experimental and numerical results is very good, except for the limit $x \to 0$ (end of the porous medium) where the local geometry changes are not accounted for in the model (the porous medium being considered homogeneous). Also, the porous medium being very anisotropic, the model used for the hot exchanger (and particularly its solid-fluid geometry shown on Fig. 2a) influences the velocity near $x = 0$. As expected, the velocity profiles are not symmetric with respect to the x-axis, and show the development of the convection boundary layer along the “hotter” section $x = 0$, in the bottom half between $r = -15$ mm until $r = 0$. The maximum velocity then decreases because when the flow reaches the top tube wall it is forced to turn. Fig. 6 shows the axial profiles of the axial velocity component obtained from experiments and numerical simulations in the same locations as Fig. 5. Here again the agreement between experimental and numerical results is very good, except for the limit $x \to 0$ (end of the porous medium). As expected, the velocity is negative for the lower part of the tube, and positive for the upper part, forced by the upward flow along the hotter porous medium end $x = 0$. Fig. 7 shows the vertical profile of the axial velocity component obtained from experiments and numerical simulations on the central vertical plane, on the right side of the porous medium, for several values of the axial position. The agreement between experiments and numerical simulations is good, but the maximum values are about 20% off. Also, as noted on Fig. 3, for $x > 150$ mm, there is a small discrepancy between numerical and experimental results.

Fig. 8 a and b show velocity vector maps and colored contours of the axial velocity component obtained from numerical simulations and experimental data from [15] on the central horizontal plane, on the right side of the porous medium. There is again excellent agreement between experiments and numerical simulations. Convective velocities are of course a lot smaller on that plane than on the central vertical plane. Fig. 8c shows the numerical and experimental transverse profile (as a function of the radial coordinate) of the axial velocity component on the central horizontal plane at $x = 60$ mm. On this plane the flow is symmetric with respect to the x-axis.

In order to better understand the convective flow observed in the measurement field of view, contours of the axial velocity component and temperature in the central vertical plane obtained from numerical simulations are shown within the porous medium on Fig. 9a and b, and left of the hot heat exchanger on Fig. 9c and d. The flow inside the porous medium is indeed character-
<figure><img src="image_144.png" /><figcaption>**Fig. 3.** Velocity vector maps and colored contours of the axial velocity component obtained from experiments (a) and numerical simulations (b) on the central vertical plane, on the right side of the porous medium. Zooms on the left part of the figures, for experimental results (c) and numerical simulations (d).</figcaption></figure>

<figure><img src="image_145.png" /><figcaption>**Fig. 4.** Calculated temperature map on the central vertical plane, on the right side of the porous medium.</figcaption></figure>
<figure><img src="image_146.png" /><figcaption>**Fig. 5.** Axial profile of the vertical velocity component on the right of the porous medium (squares: experiments, solid line: numerical simulation) on the central vertical plane at (a): $r = -15$ mm, (b): $r = -10$ mm, (c): $r = -5$ mm, (d): $r = 0$, (e): $r = 5$ mm, (f): $r = 10$ mm, (g): $r = 15$ mm.</figcaption></figure>
<figure><img src="image_147.png" /><figcaption>**Fig. 6.** Axial profile of the axial velocity component on the right of the porous medium (squares: experiments, solid line: numerical simulation) on the central vertical plane at (a): $r = -15$ mm, (b): $r = -10$ mm, (c): $r = -5$ mm, (d): $r = 0$, (e): $r = 5$ mm, (f): $r = 10$ mm, (g): $r = 15$ mm.</figcaption></figure>
<figure><img src="image_148.png" /><figcaption>Fig. 7. Vertical profile of the axial velocity component on the right of the porous medium (squares: experiments, solid line: numerical simulation) on the central vertical plane at (a): x = 0 mm, (b): x = 10 mm, (c): x = 40 mm, (d): x = 70 mm, (e): x = 100 mm, (f): x = 130 mm, (g): x = 150 mm (h): x = 180 mm.</figcaption></figure>
<figure><img src="image_149.png" /><figcaption>**Fig. 8.** Velocity vector maps and colored contours of the axial velocity component (a): experiments, (b): numerical simulations on the central horizontal plane, on the right side of the porous medium. (c): Distribution of the axial velocity component as a function of the radial coordinate (squares: experiments, solid line: numerical simulation) on the central horizontal plane at *x* = 60 mm.</figcaption></figure>

istic of natural convection flow in a cavity filled with a strongly anisotropic porous medium (Fig. 9a). The flow is constrained to remain in the horizontal plane, and the permeability being very small, conduction dominates convection (with mostly horizontal stratification of temperature), and the flow is very slow. Also, because of the porous medium's anisotropy, the model used for the hot heat exchanger influences the flow in the right fluid cavity, i.e. in the measurement field of view: For example the no-slip condition on the solid core and rings and the zero permeability in the radial and azimuthal directions put high constraints on the development of the flow in the porous medium. The flow thus

retains the geometry of the hot heat exchanger (Fig. 9a). In the left fluid cavity the convection vortex turns in the opposite direction, because the cooler boundary is on the left (Fig. 9c).

In order to get an insight in driving effects in the thermoacoustic cavity, the maximum velocity in different regions is examined thereafter in light of classically defined reference Rayleigh numbers and associated velocity scales. On the right side of the hot heat exchanger, the flow is driven by the temperature difference $\Delta T = T_{\text{hot}} - T_{\text{cold}}$. The flow in this region should then be that of a differentially heated cavity partly filled with an anisotropic porous medium. Using the axial permeability $K_{xx}$ and equivalent thermal
<figure><img src="image_150.png" /><figcaption>**Fig. 9.** Contour of the axial velocity component (a) and temperature (b) in the central vertical plane in the porous medium region. Contour of the axial velocity component (c) and temperature (d) in the central vertical plane in the cavity left of the hot heat exchanger.</figcaption></figure>
<figure><img src="image_151.png" /><figcaption>**Fig. 10.** (a): Schematic for the 2nd experimental setup. Photos for (b): The regenerator, (c): The hot heat exchanger, (d): The ambient heat exchanger, (e): Thermocouples locations and results of measurements.</figcaption></figure>

<figure><img src="image_152.png" /><figcaption>**Fig. 11.** Numerical configuration for the 2nd configuration.</figcaption></figure>
conductivity $k_p$, the porous medium Rayleigh number is $Ra_p = \frac{\rho g\beta \Delta THK_{xx}}{\mu \alpha_p}$, with $\alpha_p = \frac{k_p}{\rho c_p}$ the porous medium equivalent thermal diffusivity, and $H$ a reference height [18]. Using all geometrical and thermophysical parameters presented previously (Table 1), and taking the reference height $H = R$, this porous Rayleigh number is $Ra_p = 1.1$. The associated velocity scale is $V_{ref,p} = \frac{\alpha_p}{H} Ra_p = 0.02 \text{ m/s}$ [18]. The maximum velocity in the region on the right of the hot heat exchanger is about $0.04 \text{ m/s}$, which is indeed on the order of this velocity scale. The characteristics of the porous medium therefore determine the strength of the convection flow.

On the left of the hot heat exchanger, the flow is driven by the temperature difference between the left end of the cylindrical tube and the heat exchanger (imposed based on the experiments and equal to about $\Delta T_l = 110 \text{ K}$). The flow in this region should then be that of a differentially heated fluid cavity. The reference fluid Rayleigh number is classically given as $Ra = \frac{\rho g\beta \Delta T_l H^3}{\mu \alpha}$, with $\alpha = \frac{k}{\rho c_p}$ the fluid thermal diffusivity, and $H$ a reference height [17]. Using the same reference height $H = R$, this fluid Rayleigh number is $Ra = 3.36 \times 10^5$. The associated velocity scale is $V_{ref} = \frac{\alpha}{H} Ra^{1/2} = 0.38 \text{ m/s}$ [17]. The maximum velocity in the region on the left of the hot heat exchanger is about $0.15 \text{ m/s}$ (Fig. 9c), which is indeed on the order of this velocity scale. The strength of the convection flow in this region will then be determined by the effectiveness of the isolation of this "hot" part of the tube.

When this simple thermoacoustic engine is functioning (with plugs in place), acoustic velocities of about 2 m/s and streaming velocities of about 1 cm/s have been observed for the same applied electrical power [15]. Therefore the present results confirm

that natural convection is significant and that this phenomenon is important in understanding mean velocity measurements in such a device. This phenomenon, as acoustic streaming and jet streaming, is a source of heat flux that should be taken into account for a realistic thermal balance in such a device. Moreover, the above presented results show that a simple model as the one proposed here gives realistic results, meaning that a thermoacoustic core is satisfactorily described as a homogeneous porous medium in terms of natural convection effects. Finally, it was found that the use of Rayleigh numbers give a correct order of magnitude for the convective velocity, that can be used for a preliminary estimate. In the following section we consider a configuration slightly more complicated and more representative of an optimized system because an ambient heat exchanger is added. Our aims here are the same as above, i.e. investigate the natural convection flows in such a system both experimentally and numerically.

### 3. Natural convection flow in a cavity with ambient heat exchanger

#### 3.1. Experimental setup, measurement systems and procedure

In this second experimental setup, many factors were taken into considerations in the design process to overcome difficulties encountered in the 1st configuration. First, the resonator duct is chosen to have a square cross-sectional area to facilitate the PIV measurements in the vertical plane (x-z plane in Fig. 10a) at different transverse locations (y-direction). Second, an Ambient Heat Exchanger (AHX) is added on the right side of the regen-

<figure><img src="image_153.png" /><figcaption>**Fig. 12.** Velocity vector map overlapped on the color map of the axial velocity component at the exit of the AHX on vertical plane at $y = -19 \text{ mm}$ (top: experiments, bottom: numerical simulations).</figcaption></figure>
erator to keep its temperature constant. At the left side of the regenerator ends is placed an electric resistance as in the 1st configuration. Third, a thermal insulation is used to cover the hot part of the resonator and hence to reduce the heat loss through the walls of the resonator tube. As shown in Fig. 10a, the setup consists of a square duct with inner side length of 44 mm and the total length of the duct is 2655 mm. The left part of the duct is made of aluminum whereas the right part of the duct is made of acrylic to grant an optical access for the laser measurements. Two loudspeakers (Model: PHL - B20/2410) are attached to the ends of the duct. As before, in order to isolate the natural convection flow, no acoustic wave was generated in the resonator, that is the loudspeakers were turned off. The thermoacoustic core consists of a regenerator placed between the Hot Heat Exchanger (HHX) and the AHX. The coordinate system is also defined in Fig. 10a, with x = 0 corresponding to the right end of the AHX. As shown in Fig. 10b, the regenerator is made of stainless steel wire meshes stacked in a steel canister with a total length of 40 mm. The volumetric porosity (i.e. the ratio of the gas volume to the total volume) of the regenerator is 80%. The HHX is made of an insulated electric heater wire with total length of 450 mm and an outer diameter of 5 mm. The electrical wire is formed as shown in Fig. 10c, yielding a porosity of the hot heat exchanger of 43% and a length in the axial direction (i.e. x-axis) of 10 mm. As shown in Fig. 10d, the AHX is made of three parts. The central part, which is responsible for the heat exchange between the gas and the water, is printed (using 3D printing technology) from steel powder. The other two parts (white in Fig. 10d) are printed from plastic to avoid the heat transfer through the walls. The three parts are

glued together and assembled with the resonator duct. The central part of the AHX is made of fins with thickness of 1 mm and plate spacing of 1 mm. Also, there are two water paths where the water flow direction in one of them is in the opposite direction to the other to achieve a temperature uniformity over the cross-section of the AHX. The inner and outer diameters of the water tubes are 3 mm and 5 mm, respectively. The total length of the AHX in the axial direction (i.e. x-axis) is 20 mm and its porosity is 40.5%.

The same PIV system as for the first setup was used and a similar experimental procedure was followed. Here the input electrical power was 30 W ± 1% and the temperatures of the walls of the aluminum duct are measured (see Fig. 10e). Thermocouples are placed beneath the thermal insulation and in contact with the walls of the duct to accurately measure the wall temperatures. It was checked that thanks to the AHX, the temperature of the right part of the resonator duct (i.e. acrylic duct) remains at the ambient temperature therefore, no temperature measurements were conducted on the walls of the acrylic duct. The PIV measurements performed in this system focus on measuring the natural convection flow at the end of the AHX in the vertical plane (i.e. x-z plane) and at 7 different transverse locations (i.e. y = 0, ±5 mm, ±15 mm and ±19 mm).

The uncertainty in the temperature measurements is ± 1 % of the reported mean temperatures. The value of the spatial-average uncertainty in the velocity measurements is reduced to ± 6 % (compared to ± 13 % in the first experimental setup in which a circular duct is used) because the laser light reflections are less when the duct walls are flat (i.e. square duct).

<figure><img src="image_154.png" /><figcaption>**Fig. 13.** Velocity vector map overlapped on the color map of the axial velocity component at the exit of the AHX on vertical plane at y = -5 mm (top: experiments, bottom: numerical simulations).</figcaption></figure>
## 3.2. Numerical model

As previously, a simplified 3D configuration is introduced in order to model the experimental configuration (Fig. 11). The horizontal duct has a square cross-section (with side length equal to the inner side length of the experiment $l = 44$ mm).

The regenerator is modeled as a homogeneous parallelepipedic isotropic porous medium block of length $L_p = 40$ mm and cross-section width 38 mm, inserted in a steel case (as in the experiment). Its thermophysical characteristics are its porosity $\epsilon = 0.8$, its permeability $K_{xx} = K_{yy} = K_{zz} = K = 10^{-8}$ m$^2$ (based on experimental work performed on similar materials [26]), and its equivalent thermal conductivity $k_p = 3.16 \text{ W m}^{-1} \text{ K}^{-1}$ (calculated as $k_p = \epsilon k_{fluid} + (1 - \epsilon) k_{solid}$).

The hot heat exchanger (electric wire) is modeled as 6 vertical solid plates of axial width $e_h = 10$ mm (as in the experiment) and transverse width equal to 4 mm, maintained at given hot temperature $T_{hot} = 438$ K.

The ambient heat exchanger, placed on the right side of the porous medium, is modeled as 6 horizontal solid plates of axial width $e_c = 20$ mm (as in the experiment) and transverse width equal to 4 mm maintained at given cold (water) temperature $T_{cold} = 291$ K. In between the solid plates, the vertical fins are modeled as a homogeneous anisotropic porous medium (i.e. stack of thin plates). The characteristics of that porous medium are its porosity $\epsilon = 0.4$ (as in the experiment), its permeability assumed anisotropic $K_{xx} = K_{zz} = d_1^3/(12d_1 + d_2) = 10^{-8}$ m$^2$ (where $d_1$ is the plate's interspace and $d_2$ the plate's width, based on [18]),

$K_{yy} \approx 0$, and its thermal conductivity $k_p = 27.4 \text{ W m}^{-1} \text{ K}^{-1}$ (calculated from the porosity and thermal conductivities of gas and plates).

On the left of the HHX, in order to model the thermal insulation, the duct walls are adiabatic. Therefore the steady temperature in that region should be homogeneously equal to the hot temperature $T_{hot}$, and the natural convection flow should not extend very far. For that reason and in order to lower computational costs, the length of the duct left of the HHX is reduced to $L_l = 160$ mm.

Since the thermal insulation also covers part of the regenerator, the duct walls are adiabatic along the insulation (with axial extension identical to that used in the experiment). The temperature of the remainder of the duct wall along the regenerator is prescribed equal to $T_{amb} = 294$ K. The ambient temperature $T_{amb}$ is prescribed along the duct walls on the right side of the AHX.

The temperature difference $T_{amb} - T_{cold} = 3$ K being quite small, the natural convection flow should not extend very far. For that reason and in order to lower computational costs, the length of the duct right of the AHX is reduced so that the total computational length is smaller than that of the experiment and equal to $L = 430$ mm. The reference ambient temperature $T_{amb}$ is used to calculate the thermophysical parameters of the fluid.

The numerical domain is meshed using 64 points regularly spaced along the cross sectional sides in the y and z directions. The mesh size in the axial direction is $\Delta x = 1.6 \times 10^{-3}$ 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.5 criterion.

<figure><img src="image_155.png" /><figcaption>**Fig. 14.** Velocity vector map overlapped on the color map of the axial velocity component at the exit of the AHX on vertical plane at $y = 0$ (top: experiments, bottom: numerical simulations).</figcaption></figure>
### 3.3. Results and discussion

As for the previous configuration both simulation and experiments reach a steady state. Figs. 12–14 show the velocity vector map overlapped on the color map of the axial velocity component of the natural convection flow obtained from experiments and numerical simulations at the exit of the AHX on vertical plane for 3 different transverse locations $y = -19$ mm (Fig. 12), $y = -5$ mm (Fig. 13), and $y = 0$ (center plane, Fig. 14). Experimental and numerical results agree remarkably well, both qualitatively and quantitatively, considering the numerous simplifications and approximations made in the numerical model. As in the case of the first experiment -without AHX-, the natural convection vortices extend for a longer axial distance in the experiments than in numerical calculation. Again, this could be due to the fact that the temperatures imposed on the tube boundary are based on the measurement of the outside temperature instead of the inside. The convection vortex is in the opposite direction than in the configuration without ambient heat exchanger, and the maximum axial velocity component is about half. The direction of the flow comes from the fact that the ambient heat exchanger is very effective, so that the tubes containing the cold water are cooler than the ambient temperature on the right fluid cavity. Therefore the fluid tends to move downward in the vicinity of the AHX and then is forced to turn right at the bottom of the duct. The order of magnitude of the velocity can be compared to the natural convection reference velocity in a fluid cavity calculated using $\Delta T_r = T_{amb} - T_{cold} = 3$ K and a reference height $H = 44$ mm. The calculation gives the Rayleigh number $Ra = 1.63 \times 10^4$ and $V_{refp} = 0.07 \times 10^{-3}$ m/s. The

maximum velocity in that region is about 0.025 m/s (Figs. 12–14), which is indeed on the order of this velocity scale. In order to better understand the convection flow, contours of the calculated axial velocity component and temperature in the central vertical plane around the regenerator (the measurement field of view does not extend left of $x = 0$) are shown on Fig. 15. Analysis of this figure and of Fig. 12 shows that in this experiment the AHX prevents the development of a strong natural convection flow on the right part of the duct. Therefore the duct is clearly separated in two parts, left of the AHX and right of the AHX.

Since the regenerator is placed between the heat exchangers, the flow is similar to the natural convection flow in a porous differentially heated cavity (Fig. 15). However the regenerator's permeability being very small, so are velocities (on the order of 0.002 m/s) and isothermal lines are almost vertical, i.e. the regenerator acts like a conducting solid. The order of magnitude of the velocity can be compared to the porous natural convection reference calculated using $\Delta T = T_{hot} - T_{cold} = 147$ K, $H = 44$ mm. The calculation gives the porous Rayleigh number $Ra_p = 1.26 \times 10^{-1}$ and $V_{refp} = 2.2 \times 10^{-3}$ m/s. The maximum velocity in the regenerator is about 0.002 m/s (Fig. 15), which is indeed on the order of this velocity scale. Both the ambient heat exchanger and the isolation of the hot duct contribute to the concentration of temperature gradients in the regenerator only. The convection flow inside the regenerator depends then mostly on its permeability.

On the left of the HHX the flow is clearly driven by the flow between the mesh wires of the regenerator (Fig. 15). The temperature in that region is almost homogeneously hot.

<figure><img src="image_156.png" /><figcaption>Fig. 15. Contours of the axial velocity component (top) and temperature (bottom) in the vertical central plane (numerical simulation results).</figcaption></figure>
In sum, the simple model proposed here predicts the natural convection flow inside the entire device, with good agreement with measurements in the field of view.

## 4. Conclusion

Natural convection flows were analyzed both experimentally and numerically for two different horizontal thermoacoustic cavities, when no acoustic oscillations were present. The first cavity consisted of a cylindrical resonator filled with air and containing a thermoacoustic core made of a ceramic stack placed next to a wire resistance (hot heat exchanger). The second cavity consisted of a duct with square cross-section filled with air and containing a thermoacoustic core made of a regenerator placed between a wire resistance (hot heat exchanger) and an ambient heat exchanger (using a water stream for the heat exchange). Temperature was measured by thermocouples in several points along the cavities, and velocity fields were obtained from PIV measurements in several vertical planes (only the central vertical plane in the first experiment). The numerical model consisted in solving Navier-Stokes equations under Boussinesq assumptions inside a simplified cavity based on each experimental setup. In the first experiment, the cylindrical cavity was of similar geometry and dimensions as the experiment, the hot heat exchanger was modeled as several immersed solid bodies at given temperature, and the ceramic stack was modeled as a homogeneous anisotropic porous medium. The temperature along the cylindrical tube was prescribed as a boundary condition based on the experimental measurements. In the second experiment, the cavity was also modeled based on experiments, the regenerator being modeled as a homogeneous isotropic porous medium and the ambient heat exchanger as an immersed solid body at given temperature inside a homogeneous anisotropic porous medium. The temperature boundary conditions were also based on experimental measurements.

It was shown that the flow and temperature distributions are fully three-dimensional. They depend on the thermophysical properties of the working gas and of the porous media constituting the thermoacoustic stack or regenerator and heat exchangers, as well as on the geometrical characteristics of the device and the temperature conditions. There was excellent qualitative and quantitative agreement on the velocity fields between experimental and numerical simulation results for both configurations. Also the velocity scales observed in the different regions of the thermoacoustic cavities were in good agreement with classical scales for natural convection flows within a porous or fluid cavity. It should therefore be possible to predict the strength and patterns of natural convection flows inside functioning thermoacoustic devices using an approach similar to the one proposed here. It was noted that the observed natural convection flows could be quite significant in such devices and potentially interfere with thermoacoustic effects. We hope the fairly simple models used here will be further used by others to assess potential natural convection effects and contribute to optimal design of thermoacoustic devices. However, possible interaction and coupling between thermoacoustic and natural convection effects was not studied here and is the topic of ongoing work.

### CRediT authorship contribution statement

**Omar Hireche:** Formal analysis, Writing - original draft. **Islam Ramadan:** Conceptualization, Visualization, Writing - original draft. **Catherine Weisman:** Supervision, Validation, Writing - review & editing. **Hélène Bailliet:** Project administration, Writing - original draft. **Yann Fraigneau:** Software. **Diana Baltean-Carlès:** Project administration, Writing - original draft. **Virginie Daru:** Validation, Supervision.

### Acknowledgments

This research is financially supported by Agence Nationale de la Recherche (Projet ANR-17-CE06-0007-01). The authors also appreciate the technical assistance provided by Pascal Biais, Laurent Philippon and Philippe Szeger.

### Supplementary material

Supplementary material associated with this article can be found, in the online version, at 10.1016/j.ijheatmasstransfer.2019.119195

### References

[1] G.W. Swift, Thermoacoustics: a unifying perspective for some engines and refrigerators, ASA Press, Springer, 2002.

[2] W. Dai, E. Luo, J. Hu, H. Ling, A heat driven thermoacoustic cooler capable of reaching liquid nitrogen temperature, Appl. Phys. Lett. 86 (2005) 224103.

[3] L. Zontjens, C. Howard, A. Zander, B. Cazzolato, Feasibility study of an automotive thermo acoustics refrigerator, in: Proceedings of Acoustics, Busselton, Australia, 2005.

[4] J.A. Adeff, T.J. Hofler, Design and construction of a solar-powered, thermoacoustically driven, thermoacoustic refrigerator, J. Acoust. Soc. Am. 107 (6) (2005) L37-L42.

[5] S. Backhaus, E. Tward, M. Petach, Traveling wave thermoacoustic electric generator, Appl. Phys. Lett. 85 (2004) 1085.

[6] S. Spoelstra, M.E.H. Tijani, Thermoacoustic heat pumps for energy savings, in seminar "boundary crossing acoustics", in: Proceedings of the Acoustical Society of the Netherlands, 2005.

[7] D.L. Gardner, G.W. Swift, A cascade thermoacoustic engine, J. Acoust. Soc. Am. 114 (4) (2003) 1905-1919.

[8] Y.-L. He, H.-B. Ke, F.Q. Cui, W.Q. Tao, Explanations on the onset and damping behaviors in a standing-wave thermoacoustic engine, Appl. Therm. Eng. 58 (2013) 298-304.

[9] G. Thummes, M. Schreiber, R. Landgraf, C. Heiden, Convective heat losses in pulse tube coolers: effect of pulse tube inclination, in: R.G. Ross (Ed.), Cryocoolers 9, Springer, Boston, 1997.

[10] M.E. Poese, R.M. Keolian, R.W.M. Smith, S.L. Garrett, E.C. Mitchell, Performance estimates of a helium-based thermoacoustic-Stirling chiller, Int. J. Low-Carbon Technol. 8 (2013) i79-i84.

[11] N. Pan, S. Wang, C. Shen, Visualization investigation of the flow and heat transfer in thermoacoustic engine driven by loudspeaker, Int. J. Heat Mass Transf. 55 (2012) 7737-7746.

[12] N. Pan, S. Wang, C. Shen, A fundamental study on characteristic of thermoacoustic engine with different tilt angles, Int. J. Heat Mass Transf. 74 (2014) 228-237.

[13] P. Debesse, D. Baltean-Carlès, F. Lusseyran, M.X. François, Oscillating and streaming flow identification in a thermoacoustic resonator from under-sampled PIV measurements, Meas. Sci. Technol. 25 (16pp) (2014) 25005.

[14] E. Saint Ellier, Y. Bailly, L. Girardot, D. Ramel, P. Nika, Temperature gradient effects on acoustic and streaming velocities in standing acoustic waves resonator, Exp. Therm. Fluid Sci. 66 (2015) 1-6.

[15] I. Ramadan, H. Bailliet, J.C. Valière, Experimental investigation of the influence of natural convection and end-effects on Rayleigh streaming in a thermoacoustic engine, J. Acoust. Soc. Am. 143 (1) (2018) 361-372.

[16] D. Poulikakos, A. Bejan, Natural convection in vertically and horizontally layered porous media heated from the side, Int. J. Heat Mass Transf. 26 (1983) 1805-1814.

[17] A. Bejan, Convection Heat Transfer, Mac Graw Hill, 1985.

[18] D.A. Nield, A. Bejan, Convection in Porous Media, Springer and Business Media, New York, 2006.

[19] C. Weisman, D. Baltean-Carlès, P. Duthil, P.L. Quéré, Natural convection in a stack of horizontal plates in a differentially heated cavity, in: Proceedings of the 19th International Symposium on Transport Phenomena (ISTP-19), Reykjavik, Iceland, 2008.

[20] O. Hireche, C. Weisman, D. Baltean-Carlès, V. Daru, Y. Fraigneau, Numerical study of the effects of natural convection in a thermoacoustic device, in: Proceedings of the 24th Congres Francais de Mecanique, Brest, France, 2019.

[21] A.K. Prasad, Particle image velocimetry, Curr. Sci. 79 (1) (2000) 51-60.

[22] D.W. Peaceman, H.H. Rachford, The numerical solution of parabolic and elliptic differential equations, J. Soc. Indust. Appl. Math. 3 (1) (1955) 28-41.

[23] K. Goda, A multistep technique with implicit difference schemes for calculating two or three-dimensional cavity flows, J. Comput. Phys. 30 (1979) 76-95.

[24] J.L. Guermond, P.D. Minev, J. Shen, An overview of projection methods for incompressible flows, Comp. Meth. Appl. Mech. Eng. 195 (2006) 6011-6045.

[25] Y. Saad, Iterative methods for sparse linear systems, Society for Industrial and Applied Mathematics, second ed., 2003.

[26] S.H. Tasnim, Porous Media Thermoacoustic Stacks: Measurements and Models, University of Waterloo, Ontario, Canada, 2011 Thesis.