1. Introduction
Standing acoustic waves in homogeneous fluids have long been known to generate Eulerian mean flows. In the 19th century, the observation of air flows in a Kundt's tube by Dvorák (Reference Dvorák1876) was one of the three phenomena that led Rayleigh (Reference Rayleigh1884) to develop his pioneering analysis of acoustic streaming. A general framework, reviewed by Riley (Reference Riley2001), has since been achieved: the Reynolds stress divergence induced by a standing acoustic wave generates a streaming flow provided this stress divergence cannot be balanced by a mean pressure gradient; that is, provided the acoustic wave has non-zero vorticity. If acoustic wave attenuation takes place over a much longer time scale than the acoustic period, fluctuating vorticity is localized in thin oscillatory boundary layers, and the characteristic streaming velocity scale 
$U_{s_*}=U_*^2/a_*$, with 
$U_*$ the maximum acoustic-wave velocity and 
$a_*$ the speed of sound. This streaming velocity remains modest and limits practical applications to microfluidics, where acoustic forcing is used to mix dilute chemicals along the direction transverse to the microchannel (Bengtsson & Laurell Reference Bengtsson and Laurell2004). Recently, it has been shown that larger streaming velocities can be achieved if the solid boundaries include sharp edges having radii of curvature smaller than the thickness of the oscillating boundary layer, a promising finding albeit one that is difficult to realize experimentally (Huang et al. Reference Huang, Xie, Ahmed, Rufo, Nama, Chen, Chan and Huang2013; Zhang et al. Reference Zhang, Guo, Brunet, Costalonga and Royon2019).
Standing acoustic waves in inhomogeneous fluids drive streaming flows with radically different features. This essential outcome of the last two decades stems from experimental evidence presented by Loh et al. (Reference Loh, Hyun, Ro and Kleinstreuer2002), Hyun, Lee & Loh (Reference Hyun, Lee and Loh2005) and Stockwald, Kaestle & Ernst (Reference Stockwald, Kaestle and Ernst2014), who reported streaming velocities in stratified gases two orders of magnitude larger than the estimate 
$U_*^2/a_*$ and streaming patterns qualitatively different from that predicted by Rayleigh's theory (Dreeben & Chini Reference Dreeben and Chini2011). This change of phenomenology is also apparent in direct numerical simulations (DNS) of the compressible Navier–Stokes equations (Lin & Farouk Reference Lin and Farouk2008; Aktas & Ozgumus Reference Aktas and Ozgumus2010). Based on these results, Chini, Malecha & Dreeben (Reference Chini, Malecha and Dreeben2014) derived a new set of wave–mean-flow interaction equations that are able to capture the dynamics of streaming flows in an inhomogeneous gas. These authors demonstrated that the physical origin of the enhanced streaming flows is the baroclinic production of fluctuating vorticity, as is evident by taking the curl of the linearized Euler equation, viz.
\begin{equation} \boldsymbol{\nabla} \times \left( \rho_0 \frac{\partial \boldsymbol{u}'}{\partial t} ={-} \boldsymbol{\nabla} p' \right) \Rightarrow \frac{\partial (\boldsymbol{\nabla} \times \boldsymbol{u}')}{\partial t} = \frac{(\boldsymbol{\nabla} \rho_0) \times (\boldsymbol{\nabla} p')}{\rho_0^2},\end{equation}
where 
$\rho _0$ is the non-uniform background density, 
$\boldsymbol {u}'$ is the acoustic wave velocity, 
$p'$ is the acoustic wave pressure, 
$\boldsymbol {\nabla }$ is the spatial gradient operator and 
$t$ is the time variable. Whereas the fluctuating vorticity required to drive streaming flows in homogeneous fluids results from viscous torques and is confined to thin Stokes boundary layers, in (stably) density stratified fluids wave vorticity can readily fill the entire domain owing to its generation via an inviscid process. Accounting for this fundamental change leads to a characteristic streaming velocity 
$U_{s_*}=U_*$ in strongly inhomogeneous fluids. Given that the acoustic Mach number 
$U_*/a_*$ typically is very small compared with unity, this result provides a rationale for the large-amplitude streaming flows previously reported. Karlsen, Augustsson & Bruus (Reference Karlsen, Augustsson and Bruus2016) extended this procedure to inhomogeneities in compressibility (negligible in gases) and demonstrated that, in an inhomogeneous gas, the wave-induced Reynolds stress divergence (i.e. the acoustic force density 
$\boldsymbol {f}_{ac}$) can be expressed as
\begin{equation} \boldsymbol{f}_{ac} ={-} \tfrac{1}{2} \,\overline{|\boldsymbol{u}'|^2} \,\boldsymbol{\nabla} \rho_0,\end{equation}
where the overbar indicates a time average. This concise representation confirms that the forcing is inviscid and associated with baroclinicity (also proportional to 
$\boldsymbol {\nabla } \rho _0$, as shown in (1.1)). Moreover, this expression correctly suggests that two-way coupling between the fast acoustic waves and the slowly evolving streaming flow can be realized. Specifically, in the framework of Chini et al. (Reference Chini, Malecha and Dreeben2014) and Karlsen et al. (Reference Karlsen, Augustsson and Bruus2016), the waves drive a streaming flow that advects inhomogeneities in density that, in turn, feed back on the wave velocity field 
$\boldsymbol {u}'$. Furthermore, given that the velocity field of the standing acoustic wave depends on the background density distribution over the entire domain, this two-way coupling is spatially non-local.
In addition to manifesting more complex dynamics, this new regime of streaming is practically important. In high-intensity discharge lamps, a temperature difference of several thousand degrees exists between the arc and the tube wall, enabling baroclinic acoustic streaming to be used to improve the efficiency of these lamps dramatically (Dreeben & Chini Reference Dreeben and Chini2011; Chini et al. Reference Chini, Malecha and Dreeben2014). Similarly, streaming flows in thermoacoustic devices strongly depend on the inhomogeneous temperature distribution (Daru et al. Reference Daru, Weisman, Baltean-Carlès and Bailliet2021a, Reference Daru, Weisman, Baltean-Carlès and Baillietb). In microfluidics, properly accounting for inhomogeneities in both density and compressibility is crucial for obtaining accurate predictions of the acoustically driven mixing of different fluids (Karlsen et al. Reference Karlsen, Qiu, Augustsson and Bruus2018; Pothuri, Azharudeen & Subramani Reference Pothuri, Azharudeen and Subramani2019; Qiu et al. Reference Qiu, Joergensen, Corato, Bruus and Augustsson2021). Acoustic wave forcing also has been proposed as an experimental means to mimic and tune gravity in regions where 
$\boldsymbol{\nabla}(\overline{|\boldsymbol{u}'|^2})$ is spatially uniform, 
$\boldsymbol {f}_{ac}$ then being locally similar to gravity in the Boussinesq approximation up to a gradient term (Koulakis & Putterman Reference Koulakis and Putterman2021; Koulakis et al. Reference Koulakis, Ofek, Pree and Putterman2023), although we stress that this connection no longer holds once two-way coupling sets in. Perhaps the most promising application involves the use of acoustics to enhance the rate of heat transfer from heated objects immersed in a cooler fluid medium, especially in scenarios in which forced convection is difficult to establish or natural convection does not occur (e.g. as for cooling electronic components aboard spacecraft).
Acoustic streaming in a straight channel with a temperature differential imposed between the top and bottom boundaries is a simple yet practical configuration in which to study the effects of inhomogeneity. Many prior studies have focused on the use of liquids as the working fluid, for which variations in compressibility also contribute to the acoustic force density (Karlsen et al. Reference Karlsen, Augustsson and Bruus2016). These studies include the large set of numerical simulations performed by Kumar et al. (Reference Kumar, Azharudeen, Pothuri and Subramani2021) and Rajendran, Solomon & Subramani (Reference Rajendran, Solomon and Subramani2022) in which the effects of two-way coupling were neglected. Their results demonstrate that acoustics can either enhance or reduce the heat flux associated with natural convection depending on whether the acoustic wave displacements are, respectively, predominantly parallel or perpendicular to the channel walls. Here, we consider the second class of fluids, i.e. gases, typically air in experiments, approximated as an ideal gas in theoretical analyses. The experiments of Nabavi, Siddiqui & Dargahi (Reference Nabavi, Siddiqui and Dargahi2008) clearly demonstrate that the streaming pattern derived in the homogeneous limit by Rayleigh (Reference Rayleigh1884) is modified for temperature differences as low as a fraction of a degree. Results from DNS of the compressible Navier–Stokes equations for larger temperature differentials have been reported in Lin & Farouk (Reference Lin and Farouk2008), Aktas & Ozgumus (Reference Aktas and Ozgumus2010) and Baran et al. (Reference Baran, Machaj, Malecha and Tomczuk2022). The modest acoustic amplitudes considered in these simulations resulted in only minor modifications of the background density profile, allowing Michel & Chini (Reference Michel and Chini2019) to obtain an explicit expression for these baroclinic streaming flows by neglecting the two-way coupling. Larger acoustic amplitudes lead to two-way coupling, as analysed theoretically by Chini et al. (Reference Chini, Malecha and Dreeben2014), Karlsen et al. (Reference Karlsen, Qiu, Augustsson and Bruus2018) and Michel & Chini (Reference Michel and Chini2019) and realized experimentally by Michel & Gissinger (Reference Michel and Gissinger2021).
In the present study, the dynamics of an ideal gas undergoing standing acoustic-wave oscillations in a differentially heated channel is investigated using an extension of the wave–mean-flow interaction equations derived in Michel & Chini (Reference Michel and Chini2019). Our main objective is to characterize the strongly nonlinear regimes attained at large acoustic amplitudes that manifest two-way coupling. More specifically, we quantify the dependence of the acoustically enhanced heat flux on the geometry of the flow, i.e. on the aspect ratio 
$\delta = k_*H_*$, where 
$k_*$ is the acoustic-wave wavenumber and 
$H_*$ the height of the channel. Previous theoretical investigations have shown that the heat-transfer enhancement becomes negligible in the limit 
$\delta \rightarrow 0$ (Michel & Chini Reference Michel and Chini2019), while recent DNS performed for modest acoustic amplitudes suggest that maximum enhancement is reached for 
$\delta = O(1)$, although these simulations did not achieve strong two-way coupling (Malecha Reference Malecha2023). Moreover, the scaling of the heat flux with the wave amplitude is largely undocumented in this regime. Presuming the background temperature field is homogenized in the interior of the domain as a result of the streaming-induced mixing, baroclinicity (and, hence, the acoustic force density) will localize near the walls. The consequence of this qualitative change in behaviour on the generation of steady mechanical and thermal boundary layers and on the overall heat transfer rate is also addressed in the present investigation.
The remainder of this article is organized as follows. The two time scale wave–mean-flow system governing the two-way coupled dynamics for arbitrary (but fixed) values of 
$\delta$ is introduced in § 2. The numerical algorithm for simulating the wave–mean-flow interaction equations is described in § 3, where detailed results are also presented. A summary of our key findings and their potential implications is given in § 4.
2. Problem formulation
2.1. Flow configuration
The flow configuration is similar to that investigated by Michel & Chini (Reference Michel and Chini2019). Dimensional variables are denoted with tildes and dimensional parameters with asterisks. These variables and parameters and their definitions are summarized in table 1.
Table 1. Dimensional variables and parameters.
As illustrated in figure 1, we investigate the two-dimensional (2-D) dynamics of an ideal gas in an infinitely long channel of height 
$H_*$. No-slip and zero normal-flow boundary conditions on the velocity along with Dirichlet conditions on the temperature are imposed along each horizontal wall, with the temperature set to 
$T_*$ at 
$\tilde {y}=0$ and to 
$T_*+\Delta \varTheta _*$ at 
$\tilde {y}=H_*$. This environment is assumed to be gravity-free and, thus, the imposed temperature differential 
$\Delta \varTheta _*>0$ neither generates natural convection nor restoring forces. In the absence of acoustic forcing, a linear temperature profile is established owing to the strictly diffusive heat flux across the channel. An unspecified external agency then drives a standing acoustic wave of angular frequency 
$\omega _*$ and wavenumber 
$k_*$ along the horizontal (
$\tilde {x}$) direction to generate a streaming flow that enhances this heat flux. The horizontal (spatial) periodicity 
$2{\rm \pi} /k_*$ of the acoustic wave is presumed constant, being fixed in laboratory experiments, for example, by the finite horizontal length of the channel (or, more properly, cavity). Consequently, the angular frequency 
$\omega _*$ may vary with time as the temperature field slowly evolves. Moreover, the spatial phase of the acoustic wave is fixed by setting to zero the horizontal velocity component at 
$\tilde {x}=0$, i.e. 
$\tilde {u}(\tilde {x}=0, \tilde {y},\tilde {t})=0$, where 
$(\tilde {x},\tilde {y})$ are the horizontal and vertical coordinates and 
$\tilde {t}$ is the time variable.
Figure 1. Schematic of the 2-D system configuration, similar to Michel & Chini (Reference Michel and Chini2019). An ideal gas is confined between two horizontal, no-slip and impermeable walls separated by a distance 
$H_*$. The temperatures of the lower and upper walls are fixed at 
$T_*$ and 
$T_* + \Delta \varTheta _*$, respectively (but note that gravity is not included). A standing acoustic wave of horizontal wavenumber 
$k_*$ generates a counter-rotating cellular streaming flow spanning the channel that enhances the initially diffusive heat flux.
The flow is governed by the compressible Navier–Stokes equations, supplemented by the energy equation and the ideal gas equation of state,
$$\begin{gather} \tilde{\rho} \left[ \partial_{\tilde{t}} \tilde{\boldsymbol{u}} + ( \tilde{\boldsymbol{u}} \boldsymbol{\cdot} \tilde{\boldsymbol{\nabla}} )\tilde{\boldsymbol{u}} \right] ={-} \tilde{\boldsymbol{\nabla}} \tilde{p} + \mu_* \left[ \tilde{{\nabla}}^2 \tilde{\boldsymbol{u}} + \tfrac{1}{3}\tilde{\boldsymbol{\nabla}} ( \tilde{\boldsymbol{\nabla}} \boldsymbol{\cdot} \tilde{\boldsymbol{u}} ) \right], \end{gather}$$
$$\begin{gather}\partial_{\tilde{t}} \tilde{\rho} + \tilde{\boldsymbol{\nabla}} \boldsymbol{\cdot} \left( \tilde{\rho} \tilde{\boldsymbol{u}} \right) = 0, \end{gather}$$
$$\begin{gather}\tilde{\rho} c_{v_*} \left[ \partial_{\tilde{t}} \tilde{T} + ( \tilde{\boldsymbol{u}} \boldsymbol{\cdot} \tilde{\boldsymbol{\nabla}} )\tilde{T} \right] ={-} \tilde{p} \left( \tilde{\boldsymbol{\nabla}} \boldsymbol{\cdot} \tilde{u} \right) + \kappa_* \tilde{\nabla}^2 \tilde{T}, \end{gather}$$
$$\begin{gather}\tilde{p} = \tilde{\rho} R_* \tilde{T}, \end{gather}$$
where 
$\tilde {\rho }$, 
$\tilde {p}$, 
$\tilde {T}$ and 
$\tilde {\boldsymbol {u}}=(\tilde {u},\tilde {v})$ denote the density, pressure, temperature and velocity fields, respectively, and 
$\tilde {\boldsymbol {\nabla }} = (\partial _{\tilde {x}}, \partial _{\tilde {y}})$. As discussed in Michel & Chini (Reference Michel and Chini2019), bulk viscosity and viscous heating are neglected and the variation of dynamic viscosity with temperature is ignored in (2.1)–(2.4).
2.2. Scaling and non-dimensionalization
Dimensionless variables and parameters are introduced using the scalings given in table 2. The small parameter subsequently used in the asymptotic analysis is the acoustic Mach number 
$\epsilon = U_*/a_*$, also referred to as the inverse Strouhal number of the oscillating flow. The dimensionless temperature differential 
$\varGamma = \Delta \varTheta _* / T_*$ is taken to be fixed and order unity as 
$\epsilon \to 0$ to characterize a stratified gas in which baroclinicity plays a crucial role. This distinguished limit is similar to that adopted in Chini et al. (Reference Chini, Malecha and Dreeben2014) and Michel & Chini (Reference Michel and Chini2019). In contrast with these previous investigations, which focused on long thin channels for which 
$\delta = O (\sqrt {\epsilon })$, in this work we consider acoustic waves with horizontal wavelengths comparable to the height of the channel, i.e. 
$\delta = O(1)$. Note that the Reynolds and Péclet numbers based on the acoustic-wave oscillatory velocity and wavelength also characterize the streaming flow, which has typical velocities of order 
$U_{S_*} = U_*$ and develops structures having a spatial periodicity comparable to that of the waves. (In Chini et al. (Reference Chini, Malecha and Dreeben2014) and Michel & Chini (Reference Michel and Chini2019), cross-channel, i.e. vertical, diffusion is enhanced owing to the thinness of the channel, resulting in a disparity between the acoustic and streaming Reynolds numbers.)
Table 2. Dimensionless variables and parameters, as in the previous analyses of Chini et al. (Reference Chini, Malecha and Dreeben2014) and Michel & Chini (Reference Michel and Chini2019) except for the aspect ratio, the Reynolds number and the Péclet number, which in the present work are 
$O(1)$ quantities (asymptotically).
It is instructive to compare these scalings with typical experimental values. Using the definitions given in table 2, the acoustic cavity considered by Michel & Gissinger (Reference Michel and Gissinger2021), for instance, would be characterized by the following parameter values: 
$\epsilon = 1.8\times 10^{-4}$, 
$\delta = 1.5$, 
$\varGamma \in [0.1,0.3]$, 
$Re=134$, 
$Pe=86$ and 
$A \in [0.2,2]$ (where the wave amplitude 
$A = O(1)$ is introduced in § 2.4). Thus, the scaling requirements 
$\epsilon \ll 1$, 
$\delta = O(1)$, 
$\varGamma = O(1)$, 
$Re = O(1)$ and 
$Pe = O(1)$ encompass this set-up.
2.3. Asymptotic analysis
A multiple time scale analysis of the governing equations is performed (i) to disentangle the dynamics of the fast acoustic waves from the comparatively slowly evolving streaming flow and (ii) to consistently suppress negligible terms to simplify the numerical implementation and facilitate appropriate physical interpretation. To this end, we expand the various fields in powers of the small dimensionless parameter 
$\epsilon$:
$$\begin{gather} (u,v) = \epsilon (u_1,v_1) + O(\epsilon^2), \end{gather}$$
$$\begin{gather}P = 1 + \epsilon {\rm \pi}_1 + \epsilon^2 {\rm \pi}_2 + O(\epsilon^3), \end{gather}$$
$$\begin{gather}\varTheta = 1 + \varGamma y + \varTheta_0 + \epsilon \varTheta_1 + O(\epsilon^2), \end{gather}$$
$$\begin{gather}\rho = \rho_0 + \epsilon \rho_1 + O(\epsilon^2). \end{gather}$$
Note that the dimensionless temperature 
$\tilde {T}/T_*$ is denoted by 
$\varTheta$, cf. (2.7), rather than by 
$T$, the latter notation being reserved for the ‘slow’ time variable (see below).
In the absence of acoustic waves, 
$u=v=0$ and, in steady state, the temperature profile reduces to the conduction profile, 
$\varTheta =1+\varGamma y$, with the pressure being uniform since gravity is neglected. When excited, the leading-order acoustic velocity is 
$O(\epsilon )$, by definition since 
$U_* = \epsilon a_*$, and similarly small perturbations in the pressure and temperature fields are generated. The 
$O(1)$ temperature disturbance 
$\varTheta _0$ arises from the reorganization of temperature inhomogeneities by the strong streaming flow, as shall be made explicit in the wave–mean-flow interaction equations given subsequently.
The distinction between the waves and the streaming flow can be made based on a separation of time scales. While the acoustic wave fields exhibit fast oscillations with zero mean value, the streaming flow effectively remains constant over this time scale, evolving instead on a slow time 
$T \equiv \epsilon t$. The formal separation is captured via a Wentzel–Kramers–Brillouin–Jeffreys approximation through the introduction of a rapidly varying phase 
$\phi (t) = \varPhi (T)/\epsilon$. In this framework, any field 
$f(x,y,t)$ is expressed as 
$f(x,y,\phi,T)$, where 
$\phi$ and 
$T$ are taken to be independent variables, and thus can be decomposed as
\begin{equation} f(x,y,\phi,T ) = \bar{f}(x,y,T) + f'(x,y,\phi;T). \end{equation}Here, the streaming flow component is represented by
\begin{equation} \bar{f}(x,y,T) = \frac{1}{2{\rm \pi} n} \int_\phi^{\phi + 2 n {\rm \pi}} f(x,y,s,T) \,{\rm d}s \end{equation}
for sufficiently large positive integer 
$n$ (and arbitrary 
$\phi$). In contrast, 
$f'(x,y,\phi ;T)$ describes the acoustic wave of zero mean value (
$\overline {f'}(x,y,\phi ;T) = 0$). Finally, the instantaneous angular frequency 
$\omega (T)$ satisfies
\begin{equation} \omega(T) = \frac{\mathrm{d}\phi}{\mathrm{d}t} = \frac{\mathrm{d}\varPhi}{\mathrm{d}T} = \omega_0(T) + O(\epsilon). \end{equation}2.4. Leading-order multiscale wave–mean-flow interaction equations
This multiscale analytical approach has been used by Chini et al. (Reference Chini, Malecha and Dreeben2014) to derive a closed set of equations describing the coupled dynamics of the acoustic wave and the streaming flow. The present formulation differs only in the scaling of the aspect ratio and diffusive terms, which can be easily traced in the derivation. Here, we simply state the resulting leading-order equations.
The dynamics of the streaming flow is governed by the following mean-flow system:
\begin{gather} \bar{\rho}_0 \left( \partial_T \bar{u}_1 + \bar{u}_1 \partial_x \bar{u}_1 + \bar{v}_1 \partial_y \bar{u}_1 \right) ={-} \frac{\partial_x \bar{\rm \pi}_2 }{\gamma} - \left[ \partial_x \left(\bar{\rho}_0 \overline{u_1^{'2}} \right) + \partial_y \left(\bar{\rho}_0 \overline{u_1'v_1'} \right) \right] \nonumber\\ + \frac{1}{Re} \left[ \partial_{xx} \bar{u}_1 + \frac{1}{\delta^2} \partial_{yy} \bar{u}_1 + \frac{1}{3} \left( \partial_{xx} \bar{u}_1 + \partial_{xy} \bar{v}_1 \right) \right], \end{gather}
\begin{gather} \bar{\rho}_0 \left( \partial_T \bar{v}_1 + \bar{u}_1 \partial_x \bar{v}_1 + \bar{v}_1 \partial_y \bar{v}_1 \right) ={-} \frac{\partial_y \bar{\rm \pi}_2 }{\gamma \delta^2} - \left[ \partial_x \left(\bar{\rho}_0 \overline{u_1' v_1'} \right) + \partial_y \left(\bar{\rho}_0 \overline{v_1^{'2}} \right) \right]\nonumber\\ + \frac{1}{Re} \left[ \partial_{xx} \bar{v}_1 + \frac{1}{\delta^2} \partial_{yy} \bar{v}_1 + \frac{1}{3\delta^2} \left( \partial_{xy} \bar{u}_1 + \partial_{yy} \bar{v}_1 \right) \right], \end{gather}
\begin{gather} \partial_T \bar{\rho}_0 + \partial_x \left( \bar{\rho}_0 \bar{u}_1 \right) + \partial_y \left( \bar{\rho}_0 \bar{v}_1 \right) = 0, \end{gather}
\begin{gather} \bar{\rho}_0 \left[ \partial_T \bar{\varTheta}_0 + \bar{u}_1 \partial_x \bar{\varTheta}_0 + \bar{v}_1 \left(\varGamma + \partial_y \bar{\varTheta}_0 \right) \right] = (1-\gamma) \left(\partial_x \bar{u}_1 + \partial_y \bar{v}_1 \right) \nonumber\\ + \frac{\gamma}{Pe} \left( \partial_{xx} \bar{\varTheta}_0 + \frac{1}{\delta^2} \partial_{yy} \bar{\varTheta}_0 \right), \end{gather}
\begin{gather} \bar{\rho}_0 = \frac{1}{1+\varGamma y + \bar{\varTheta}_0}. \end{gather}
This system describes the wave-averaged motion of a compressible ideal gas driven by an acoustic force density 
${\boldsymbol {f}}_{ac}=- \boldsymbol {\nabla } \boldsymbol {\cdot } (\bar {\rho }_0 \overline {\boldsymbol {u}_1'\boldsymbol {u}_1'})$ that, for an ideal gas, can be reduced to 
$-(1/2) \overline {|\boldsymbol {u}_1'|^2} \boldsymbol {\nabla } \bar {\rho }_0$ (up to a pressure gradient), as reported in (1.2) (Karlsen et al. Reference Karlsen, Augustsson and Bruus2016). If the acoustic force density is fixed, the problem reduces to that of forced convection.
In baroclinic acoustic streaming, however, the situation is considerably more subtle, as the slowly evolving temperature disturbance 
$\bar {\varTheta }_0$, and hence the density field 
$\bar {\rho }_0$, feeds back on the acoustic wave. This coupling is clearly evident in the corresponding equations for the waves, viz.
$$\begin{gather} \omega_0 \bar{\rho}_0 \partial_\phi u_1' + \frac{1}{\gamma} \partial_x {\rm \pi}_1' = 0, \end{gather}$$
$$\begin{gather}\omega_0 \bar{\rho}_0 \partial_\phi v_1' + \frac{1}{\gamma \delta^2} \partial_y {\rm \pi}_1' = 0, \end{gather}$$
$$\begin{gather}\omega_0 \partial_\phi \rho_1' + \partial_x \left( \bar{\rho}_0 u_1' \right) + \partial_y \left( \bar{\rho}_0 v_1'\right) = 0, \end{gather}$$
$$\begin{gather}\omega_0 \partial_\phi \varTheta_1' + u_1' \partial_x \bar{\varTheta}_0 + v_1' \left(\varGamma + \partial_y \bar{\varTheta}_0 \right) + \frac{\gamma -1}{\bar{\rho}_0 } \left( \partial_x u_1' + \partial_y v_1' \right)=0, \end{gather}$$
$$\begin{gather}{\rm \pi}_1' = \frac{\rho_1'}{\bar{\rho}_0 } + \bar{\rho}_0 \varTheta_1'. \end{gather}$$
On the fast time scale, the dynamics of the acoustic waves is governed by a linear homogeneous system, the solution of which consists of a sum of modes of various amplitudes and angular frequencies. Here, however, we assume that only one acoustic mode is externally forced, i.e. only one wave has finite amplitude. This acoustic mode has dimensional horizontal wavenumber 
$k_*$ and is the solution of the eigenvalue problem for the wave associated with the lowest eigenvalue (i.e. the smallest angular frequency) that also has a non-uniform horizontal structure. Thus, any oscillatory field 
$f_1'$, representing 
$(u_1', v_1', {\rm \pi}_1', \varTheta _1', \rho _1')$, can be expressed as
\begin{equation}
f_1'(x,y,\phi;T) = \frac{A(T)}{2}\left[\hat{f}_1(x,y;T) {\rm e}^{{\rm i}\phi} + \mathrm{c.c.} \right],
\end{equation}
where 
$\mathrm {c.c.}$ denotes the complex conjugate, 
$A(T)$ is the amplitude of the mode and 
$\hat {f}_1$ is a complex function that characterizes its spatial structure. For 
$A(T)$ to be uniquely specified, a normalization condition must be imposed; we opt to require the eigenmodes to satisfy
\begin{equation} \underset{x\in[0,2{\rm \pi}], y\in [0,1]}{\mathrm{max}} \left \vert \hat{\rm \pi}_1(x,y,T)\right \vert = 1. \end{equation}
In practice, this amplitude 
$A$ can be readily inferred experimentally from the (dimensional) steady-state amplitude of the pressure oscillations measured by a sensor at 
$(x_s,y_s)$, which yields 
$\epsilon A \vert \hat {{\rm \pi} }(x_s,y_s)\vert p_*$. In the limit of very narrow channels (Michel & Chini Reference Michel and Chini2019), the eigenvalue problem for the lowest mode can be reduced to a one-dimensional problem (in 
$x$) that can be solved analytically for a linear temperature profile. Moreover, in that limit, the multiple scale analysis can be carried out to next order to obtain an equation governing the temporal evolution of the modal amplitude 
$A(T)$. Unfortunately, the corresponding higher-order analysis cannot be easily executed for the fully 2-D eigenvalue problem obtained here. The amplitude of the acoustic mode therefore will be assumed to be constant, i.e. an external control parameter, in this investigation for simplicity. Experimentally, this specification could be achieved via slowly varying acoustic forcing mechanisms controlled by the feedback of a pressure sensor.
3. Numerical simulations of the wave–mean-flow equations
3.1. Methods
The two time scale quasilinear wave–mean-flow interaction system obtained in § 2.4 is integrated numerically over the slow time variable 
$T$. Specifically, we solve the initial-value problem (2.12)–(2.16) governing the slow evolution of the streaming fields using a third-order four-stage Runge–Kutta scheme implemented in the spectral computational framework Dedalus (Burns et al. Reference Burns, Vasil, Geoffrey, Oishi, Lecoanet and Brown2020). Equations (2.14) and (2.15) are reformulated to yield a divergence condition on the streaming velocity field and a modified energy equation, viz.
$$\begin{gather} \partial_x \bar{u}_1 + \partial_y \bar{v}_1 = \frac{1}{Pe} \left( \partial_{xx} \bar{\varTheta}_0 + \frac{1}{\delta^2} \partial_{yy} \bar{\varTheta}_0 \right), \end{gather}$$
$$\begin{gather}\bar{\rho}_0 \left[ \partial_T \bar{\varTheta}_0 + \bar{u}_1 \partial_x \bar{\varTheta}_0 + \bar{v}_1 \left(\varGamma + \partial_y \bar{\varTheta}_0 \right) \right] = \frac{1}{Pe} \left( \partial_{xx} \bar{\varTheta}_0 + \frac{1}{\delta^2} \partial_{yy} \bar{\varTheta}_0 \right). \end{gather}$$
At each (slow) time step, the acoustic fields 
$\hat {u}_1$ and 
$\hat {v}_1$ are required to evaluate the wave-induced Reynolds stress divergence. Accordingly, the set of equations (2.17)–(2.21) governing the fast wave dynamics is consolidated to form an eigenvalue problem for the fluctuating pressure field:
\begin{equation}
\partial_{x}\left[\frac{1}{\bar{\rho}_0} \partial_{x} {\rm \pi}_1'\right] + \delta^{{-}2}\partial_{y}
\left[\frac{1}{\bar{\rho}_0} \partial_{y} {\rm \pi}_1'\right] = \omega_{0}^2\,\partial_{\phi \phi}{\rm \pi}_1'.
\end{equation}
We solve, at each time step, the 2-D eigenvalue problem (3.3) in MATLAB using a Fourier–Chebyshev collocation discretization (Trefethen Reference Trefethen2000; Driscoll, Hale & Trefethen Reference Driscoll, Hale and Trefethen2014). In both the Dedalus and MATLAB codes, Fourier series in the periodic 
$x$-direction and Chebyshev polynomials in the wall-normal 
$y$-direction are used to represent all field variables. The number of modes is chosen to correspond to a 
$122 \times 122$ spatial grid resolution, and the time step 
$\Delta T=0.01$. Numerical convergence in space and time has been verified.
To model a physical experiment in which a specific acoustic mode is forced for time 
$T \geqslant 0$, the steady basic conduction state (
$\bar {\varTheta }_0 = \bar {v}_0=\bar {u}_0=0$) is taken as the initial condition. The solution of the eigenvalue problem targets the horizontal standing wave of smallest non-zero eigenvalue. For particular values of 
$\delta$, standing waves along the vertical direction may share the same angular frequency as the horizontal mode of interest, but in practice the nature of the forcing (e.g. oscillating horizontal or vertical walls) would ensure that only one of these modes is excited. Both the resulting eigenfunction and eigenfrequency can vary with 
$T$. The results discussed in this section are obtained for various acoustic wave amplitudes 
$A$ and aspect ratios 
$\delta$, keeping the following parameters fixed:
\begin{equation} \varGamma = 0.3, \quad\gamma = 1.4, \quad Re = 2500, \quad Pe = 1775.\end{equation}
This set of parameters corresponds to using air as the working fluid, for which 
$\gamma = 1.4$ and the Prandtl number 
$Pe/Re=0.71$, inside an acoustic cavity similar to the one employed by Michel & Gissinger (Reference Michel and Gissinger2021) but with larger Reynolds and Péclet numbers to clearly reveal strongly nonlinear dynamical behaviour. We monitor the top- and bottom-wall Nusselt numbers, defined as the ratio of the respective heat fluxes at the upper and lower channel walls to the strictly conductive heat flux and evaluated as
\begin{equation} Nu_{t} (T)=1+ \frac{1}{2{\rm \pi} \varGamma} \int_0^{2{\rm \pi}} \partial_y \bar{\varTheta}_0(x,1,T) \,{{\rm d}x}, \quad Nu_{b} (T)=1+ \frac{1}{2{\rm \pi} \varGamma} \int_0^{2{\rm \pi}} \partial_y \bar{\varTheta}_0(x,0,T)\,{{\rm d}x}, \end{equation}
where the subscripts ‘t’ and ‘b’ denote ‘top’ and ‘bottom’. In a steady state, integrating (3.1) over the entire domain yields the expected equality 
$Nu_{t}(\infty )=Nu_{b}(\infty )\equiv Nu$.
This numerical approach, in which an eigenvalue problem is solved at each time step to account for the slow evolution of the acoustic mode, has been introduced in Karlsen et al. (Reference Karlsen, Qiu, Augustsson and Bruus2018) and Michel & Chini (Reference Michel and Chini2019) and is necessary to consistently capture the acoustic force density 
$-(1/2) \overline {|\boldsymbol {u}_1'|^2} \boldsymbol {\nabla } \bar {\rho }_0$, since 
$\boldsymbol {u}_1'$ depends non-locally on the entire background density field 
$\bar {\rho }_0(x,y,T)$. Nevertheless, the relevance of regularly updating the acoustic wave fields remains insufficiently documented, leading to its arbitrary omission in most subsequent investigations. To quantitatively assess the importance of two-way coupling, we therefore also performed a set of numerical simulations, referred to as ‘one-way coupled’, in which the temporal integration is carried out without solving the eigenvalue problem at each time step. Instead, 
$\boldsymbol {f}_{ac}(x,y,T)$ is set to 
$-(1/2) \overline {|\boldsymbol {u}_1'|^2} \boldsymbol {\nabla } \bar {\rho }_0$ evaluated at the initial time 
$T=0$. We acknowledge that several alternative ‘one-way coupled’ approximations can be proposed, such as (i) only fixing 
$\overline {|\boldsymbol {u}_1'|^2}$ in time and evaluating 
$\boldsymbol {f}_{ac}$ using the actual density gradient 
$\boldsymbol {\nabla } \bar {\rho }_0$ or (ii) setting the fluctuating pressure 
${\rm \pi} '(x,y,T,\phi ) = A \cos (x) \cos (\phi )$ (i.e. a single Fourier mode in 
$x$ that is independent of 
$y$) and evaluating the acoustic force density with the fluctuating acoustic velocity obtained from (2.17) and the current density gradient. The latter option may be appropriate for well-mixed steady states realized at large dimensionless amplitude 
$A$ (see § 3.4).
3.2. Convergence to a steady state
For small and intermediate acoustic wave amplitudes 
$A$, our numerical simulations reveal that the system converges to parameter-dependent steady states. A typical run for 
$\delta =4$ and 
$A=4$ shows that a few hundred slow time units are required to reach this steady state; see figure 2 and the animation included in the supplemental material available at https://doi.org/10.1017/jfm.2024.744. This requirement highlights the importance of using a multiple-scale algorithm, since a single time unit corresponds in practice to hundreds or thousands of acoustic cycles (cf. 
$\epsilon \approx 10^{-4}$ for the set-up of Michel & Gissinger (Reference Michel and Gissinger2021) discussed in § 2.2). During the transient regime, the top and bottom heat fluxes evolve non-monotonically, ultimately converging to a common value significantly larger than the conductive heat flux. The variation of 
$Nu$ with the amplitude of the wave 
$A$ and with the aspect ratio 
$\delta$ will be discussed in § 3.4.
Figure 2. Time series for 
$A=4$ and 
$\delta =4$ of (a) the top and bottom Nusselt numbers 
$(Nu_{t},Nu_{b})$ and (b) of the acoustic-wave angular frequency 
$\omega _{0}$. The total steady-state temperature field 
$1+ \varGamma y + \bar {\varTheta }_{0}$ shown in (c) exhibits strong variations in 
$x$ associated with localized jets at 
$x=\lbrace 0,{\rm \pi} /2, {\rm \pi}, 3{\rm \pi} /2, 2{\rm \pi} \rbrace$ and boundary layers in 
$y$ close to both walls.
This steady state exhibits both thermal boundary layers and vertical jets. Low-amplitude waves (
$A\ll 1$) generate smooth cellular structures (not shown here, for brevity), similar to the ones computed theoretically and numerically in Michel & Chini (Reference Michel and Chini2019) for 
$\delta \ll 1$ and in the DNS of Lin & Farouk (Reference Lin and Farouk2008), Aktas & Ozgumus (Reference Aktas and Ozgumus2010), Baran et al. (Reference Baran, Machaj, Malecha and Tomczuk2022) and Malecha (Reference Malecha2023).
To gain a more detailed understanding of these steady states, we further analyse the Reynolds stress divergence, also referred to as the acoustic force density 
$\boldsymbol {f}_{ac}$, resulting from the acoustic wave. This force density, corresponding to the second terms on the right-hand sides of (2.12)–(2.13), is balanced by mean inertia, mean viscous forces and the mean pressure gradient. Crucially, the part of this wave forcing that actually sustains the streaming flow, i.e. the part that is not balanced by a mean pressure gradient, is given by the curl 
$\boldsymbol {\nabla } \times \boldsymbol {f}_{ac}$. In figure 3, we plot this (scalar) field at the initial time, when the mean temperature distribution corresponds to the linear conduction profile, and in the steady state. The striking difference confirms the two-way coupling between the waves and the streaming flow in this system, in contrast to the usual forced convection configuration. A streaming flow generated by the curl of the acoustic force modifies the mean density distribution and, hence, the properties of the acoustic waves and therefore also 
$\boldsymbol {\nabla }\times \boldsymbol {f}_{ac}$.
Figure 3. Evolution for 
$A=4$ and 
$\delta = 4$ of the curl of the acoustic force density 
$\boldsymbol {\nabla } \times \boldsymbol {f}_{ac}$: the initial condition (a) and the steady state (b). Here 
$\boldsymbol {e}_{z} \equiv \boldsymbol {e}_{x}\times \boldsymbol {e}_{y}$, where 
$\boldsymbol {e}_{x}$ and 
$\boldsymbol {e}_{y}$ are unit vectors in the 
$x$ and 
$y$ directions, respectively. The feedback from the evolving streaming density field 
$\overline {\rho }_{0}$ leads to the localization of 
$\boldsymbol {\nabla }\times \boldsymbol {f}_{ac}$ near the upper and lower walls.
As noted in § 1, the generation of fluctuating vorticity in a gas is key to understanding streaming flows. In a homogeneous gas, a non-zero acoustic-wave vorticity implies that the acoustic force density
\begin{equation} \boldsymbol{f}_{ac} \equiv{-} \boldsymbol{\nabla}\boldsymbol{\cdot} \left(\bar{\rho}_0\,\overline{\boldsymbol{u}_1'\boldsymbol{u}_1'}\right)={-}\bar{\rho}_0\, \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\overline{\boldsymbol{u}_1'\boldsymbol{u}_1'}\right)= \text{gradient term} - \bar{\rho}_0\,\overline{\left(\boldsymbol{\nabla} \times \boldsymbol{u}_1'\right) \times \boldsymbol{u}_1'} \end{equation}generically cannot be balanced by a gradient term and, thus, may drive a streaming flow. Acoustic-wave vorticity in that case is usually generated by viscous and/or thermal diffusion in thin oscillatory boundary layers. In an inhomogeneous and inviscid gas, the acoustic force density instead reduces to (Karlsen et al. Reference Karlsen, Augustsson and Bruus2016)
\begin{align} \boldsymbol{f}_{ac} &\equiv{-} \boldsymbol{\nabla} \boldsymbol{\cdot} \left(\bar{\rho}_0\,\overline{\boldsymbol{u}_1'\boldsymbol{u}_1'}\right)= \text{gradient term} - \tfrac{1}{2}\,\overline{\left\vert \boldsymbol{u}_1' \right\vert^2}\,\boldsymbol{\nabla} \bar{\rho}_{0} \end{align}
\begin{align} &= \text{gradient term} - \tfrac{1}{2}\left(\bar{\rho}_0\,\overline{\left(\boldsymbol{\nabla} \times \boldsymbol{u}_1'\right) \times \boldsymbol{u}_1'} + \overline{( \boldsymbol{u}_1' \boldsymbol{\cdot}\boldsymbol{\nabla} \overline{\rho}_0) \boldsymbol{u}_1'} \right), \end{align}
where (3.8) can be derived directly from (3.7) and the leading-order fluctuating equations. In the absence of acoustic-wave vorticity, i.e. in the absence of fluctuating baroclinicity 
$\boldsymbol {\nabla }\overline {\rho }_0 \times \boldsymbol {\nabla } {\rm \pi}_1' = 0$, the steady-state contribution to 
$\boldsymbol {\nabla }\times \boldsymbol {f}_{ac}$ from the term 
$\overline {(\boldsymbol {u}_1' \boldsymbol {\cdot } \boldsymbol {\nabla }\bar {\rho }_0) \boldsymbol {u}_1'}$ in (3.8) is negligible. Thus, even in inhomogeneous gases, acoustic-wave vorticity plays a crucial role in driving streaming flows. Physical insights follow from and qualitative predictions can be made based on this understanding. First, it accounts for the localization of 
$\boldsymbol {\nabla }\times \boldsymbol {f}_{ac}$ close to the walls, as evident in figure 3. As shown in figure 4, 
$\boldsymbol {\nabla } \times \boldsymbol {f}_{ac}$ is significant only in regions where acoustic-wave vorticity also is large in magnitude; and wave vorticity is confined to thermal boundary layers, since the fluctuating isobars and the mean isopycnals are almost orthogonal there (see (1.1)). Second, this understanding explains why the orientation of the standing acoustic wave is crucial in such systems. Consider, for example, a configuration in which the initial conduction state of the present set-up is perturbed with a standing wave oscillating along the vertical rather than horizontal direction: the fluctuating isobars and the mean isopycnals would then be approximately aligned, resulting in zero acoustic vorticity and, consequently, negligible baroclinic acoustic streaming (see also Kumar et al. (Reference Kumar, Azharudeen, Pothuri and Subramani2021)).
Figure 4. Comparison for 
$A=4$ and 
$\delta = 4$ of the normalized steady-state amplitudes of (a) the curl of the Reynolds stress divergence 
$\boldsymbol {\nabla }\times \boldsymbol {f}_{ac}$ and (b) the acoustic-wave vorticity 
$\boldsymbol {\nabla }\times \boldsymbol {u}_1'$: (a) 
$|\boldsymbol {\nabla }\times\boldsymbol {f}_{ac}|/|\boldsymbol {\nabla }\times \boldsymbol {f}_{ac}|_{max}$; (b) 
$\overline {|\boldsymbol {\nabla }\times \boldsymbol {u}_1^{\prime }|}/\overline {|\boldsymbol {\nabla }\times \boldsymbol {u}_1^{\prime }|}_{max}$. Since a horizontal standing acoustic mode is considered, the fluctuating isobars (isovalues of 
$p'$) are essentially vertical – see inset in (b) – and acoustic vorticity is therefore localized where vertical gradients of density exist (see (1.1)), i.e. in the top and bottom boundary layers.
3.3. Dependence on wave amplitude
The influence of the acoustic wave amplitude 
$A$ is investigated through a suite of numerical simulations performed for the same fixed aspect ratio 
$\delta = 4$ but for various values of 
$A \in [2 \times 10^{-3}, 4]$. All simulations converge to steady states, with the Nusselt number 
$Nu$ plotted in figure 5(a). The results are compared with one-way coupled simulations (defined in § 3.1) that also converge to steady states for 
$A \in [2 \times 10^{-3}, 0.2]$. In the range 
$A \in [1, 4]$, the one-way coupled dynamics evolves to statistically stationary states, and the corresponding Nusselt numbers are obtained by time-averaging over more than 
$10^3$ slow time units.
Figure 5. (a) Steady-state Nusselt number 
$Nu$ (minus one) as a function of the acoustic wave amplitude 
$A$ for 
$\delta = 4$, along with the asymptotic prediction 
$Nu-1 \propto A^4$ that can be derived in the limit 
$A \ll 1$. (b) Steady-state Nusselt number 
$Nu$ as a function of the aspect ratio 
$\delta$ for both 
$A=4$ and 
$A=0.01$. The one-way coupled simulations, in which the evolution of the acoustic waves is neglected, provide accurate results only in the small amplitude limit 
$A \ll 1$. The inset (b.1) shows 
$Nu-1$ (rather than 
$Nu$) versus 
$\delta$.
In the limit of small acoustic-wave amplitude 
$A \ll 1$, the acoustic force density 
$\boldsymbol {f}_{ac}\propto A^2$ generates weak streaming flows and correspondingly small changes in the background temperature fields. Although of limited practical interest since the resulting Nusselt numbers are close to unity, this regime can be easily analysed theoretically by neglecting the evolution of the acoustic wave properties, i.e. by assuming that the one-way coupled dynamics provides an accurate approximation. Similarly to the analysis of Michel & Chini (Reference Michel and Chini2019), which assumes 
$\delta \ll 1$, we obtain for 
$\delta = O(1)$ that
\begin{equation} Nu -1 \underset{A \ll 1}{ = } A^4 Re^2 Pe^2 \gamma^{{-}4} F(\delta,\varGamma), \end{equation}
with 
$F(\delta,\varGamma )$ a function that has to be computed numerically; in the limit 
$\delta \ll 1$, 
$F(\delta,\varGamma ) = \delta ^8 G(\varGamma )$ for some function 
$G(\varGamma )$, as predicted by Michel & Chini (Reference Michel and Chini2019). The results of the numerical simulations reported in figure 5 support the 
$A^4$ scaling for 
$A \leqslant 0.02$, a range for which one-way and two-way coupled simulations yield nearly indistinguishable results.
Figure 5 also reveals that the one-way coupled numerical simulations and resulting theoretical prediction (3.9) fail to capture the two-way coupling that emerges for 
$A \geqslant 0.1$. In this regime, 
$Nu$ is significantly larger than unity, and accurate results can be obtained only by using the full machinery of the multiple-scale numerical algorithm implemented in our numerical simulations.
3.4. Dependence on aspect ratio
The impact of varying the aspect ratio 
$\delta = k_* H_*$ is investigated with a set of numerical simulations performed for fixed 
$A = 0.01$ or 4 and varying 
$\delta \in [ 0.25, 10 ]$. This variation can be interpreted as being achieved, for example, by adjusting the channel height 
$H_*$, since 
$\delta$ is the only parameter involving 
$H_*$ (see table 2).
The Nusselt numbers reported in figure 5(b) for small acoustic amplitude 
$A=0.01$ reach a maximum for 
$\delta \simeq 3$. Since the acoustic amplitude 
$A \ll 1$, the evolution of the acoustic wave properties can be neglected, and the one-way coupled simulations accurately describe this regime. The Nusselt number is found to converge to unity as 
$\delta \rightarrow 0$, following the asymptotic scaling law 
$Nu-1 \propto \delta ^8$ derived in Michel & Chini (Reference Michel and Chini2019), which evidently here holds for 
$\delta \leqslant 1$ (see the inset in figure 5b). For 
$A=4$, the feedback on the acoustic waves is essential, and two-way coupled simulations are required to reach a consistent steady state. In this case, the Nusselt number still exhibits significant variation with the aspect ratio, with a maximum value of 
$Nu = 56.4$ at 
$\delta = 7$. In this two-way coupled regime, the aspect ratio that maximizes the Nusselt number is expected to depend on all other parameters (
$\varGamma, Re, Pe, \gamma$) in a complicated fashion. Given that, in practice, the height of the acoustic cavity is not easily modified once an experimental apparatus is built, these simulations are particularly valuable for identifying optimal parameters beforehand.
The evolution of the streaming fields with 
$\delta$ is depicted in figures 6 and 7. For the smallest value of 
$\delta$ (
${=}0.25$), diffusion in the wall-normal (
$y$) direction prevails over inertia (see the factor 
$1/\delta ^2$ in (2.12) and (2.13)) and smooth cellular structures are observed, similar to the ones reported in Lin & Farouk (Reference Lin and Farouk2008), Aktas & Ozgumus (Reference Aktas and Ozgumus2010), Baran et al. (Reference Baran, Machaj, Malecha and Tomczuk2022) and Malecha (Reference Malecha2023). For larger values of 
$\delta$, top and bottom boundary layers accompanied by vertical jet-like structures are generated. For such flows, the 
$x$-averaged temperature 
$1+\varGamma y + \langle \bar{\varTheta}_0(x,y) \rangle _x$, where 
$\langle ({\cdot })\rangle _x$ denotes an 
$x$-average, does not monotonically vary from 1 at 
$y=0$ to 
$1+ \varGamma$ at 
$y=1$, but instead displays thin regions of reversed temperature gradient. The same feature is observed in quasilinear models of buoyancy-driven convection (Herring Reference Herring1968; O'Connor, Lecoanet & Anders Reference O'Connor, Lecoanet and Anders2021). In these reduced models, convection is driven by the interaction between hydrodynamic modes, assumed to be linearized solutions that evolve on a fast time scale determined by an eigenvalue problem in which the more slowly evolving temperature arises as a non-constant coefficient, and the slowly evolving temperature field itself, which is forced by flux divergences analogous to a wave-induced Reynolds stress divergence. Note that these emerging viscous and thermal boundary layers have a thickness that remains large compared with that of the oscillatory Stokes layers, induced by the no-slip condition at each wall, exhibited by the acoustic velocity field. These Stokes layers, which are passive in our analysis – only generating a higher-order (i.e. weaker) streaming flow – and, thus, self-consistently not resolved in our model, are of dimensional thickness 
$\delta _{BL}^* = \sqrt {\mu _*/(\rho _*a_*k_*})$. The corresponding dimensionless expression is 
$\delta _{BL}^* /H_* = \delta ^{-1}\sqrt {\epsilon /Re}$.
Figure 6. Steady-state streaming velocity fields for 
$A=4$ and various aspect ratios 
$\delta$. The white lines are isovalues of the mass current potential 
$\phi _\rho$, defined such that 
$\bar {\rho }_0 \bar {\boldsymbol {u}}_1 = \boldsymbol {\nabla } \phi _\rho$, and the colours correspond to the normalized velocity magnitude 
$(|\bar {\boldsymbol {u}}_1|/|\bar {\boldsymbol {u}}_1|_{max})$. The ratio of the height to width of each panel is set to the respective value of 
$\delta$ to facilitate qualitative comparison.
Figure 7. Wall-normal profiles of steady-state streaming temperature and horizontal velocity for 
$A=4$ and various aspect ratios 
$\delta$. (a) The total streaming temperature averaged over the horizontal 
$x$ direction, 
$1 + \varGamma y + \langle \bar{\varTheta}_0(x,y)\rangle _x$. (b) The streaming 
$x$-velocity component at a fixed location 
$x = {\rm \pi}/4$, i.e. 
$\bar{u}_1(x={\rm \pi} /4,y)$. The smooth profiles observed for 
$\delta = 0.25$ develop viscous and thermal boundary layers as 
$\delta$ increases.
4. Conclusion
The interaction between a thermally stratified gas and a standing acoustic wave has been investigated in a channel with top and bottom walls of fixed but differing temperatures. As shown in Chini et al. (Reference Chini, Malecha and Dreeben2014), Karlsen et al. (Reference Karlsen, Augustsson and Bruus2016), Karlsen et al. (Reference Karlsen, Qiu, Augustsson and Bruus2018) and Michel & Chini (Reference Michel and Chini2019), a two-way coupling develops in such inhomogeneous systems between the streaming flow, which is forced by the standing acoustic wave, and the wave, whose properties are modified by the varying background temperature. The present study extends the work of Michel & Chini (Reference Michel and Chini2019) to characterize the resulting heat fluxes in channels having heights comparable to the acoustic-wave wavelength. Moreover, the numerical simulations reported here highlight several features of baroclinic streaming that are not clearly apparent in the literature.
For the given geometry, the evolution of the Nusselt number 
$Nu$, i.e. the heat flux normalized by that realized in the absence of the acoustic wave, is of primary interest. In contrast with scaling laws derived in the limit of small acoustic amplitudes, see (3.9), and with previous investigations in narrow channels (Michel & Chini Reference Michel and Chini2019), we report values of 
$Nu$ significantly larger than unity. This finding supports the practical potential of acoustics to enhance heat transfer in situations where forced convection is difficult to establish. The numerical algorithm that has been developed, coupling an initial value problem to an eigenvalue problem, can be used to determine the geometry that maximizes the Nusselt number for given acoustic amplitudes and gas properties, as shown for two specific sets of dimensionless parameters in figure 5(b). Moreover, the results reported here constitute a reference case for future investigations seeking to further increase 
$Nu$ by tuning the boundary conditions (e.g. enforcing a constant heat flux instead of a constant temperature along the upper boundary to better model a heat source, or allowing for a non-zero wall-normal velocity to model a porous substrate).
We also document several features of baroclinic streaming that are expected to be generic and to apply, for instance, to the interaction of a standing acoustic wave with air surrounding a hot (or cold) cylinder or sphere or filling a half-plane above a warmer horizontal wall. In each scenario, these features can be traced to the baroclinic generation of acoustic-wave vorticity. In strongly stratified gases, streaming flows are forced where the curl of the acoustic force density is non-zero; this only occurs where the acoustic-wave vorticity does not vanish, i.e. where the fluctuating isobars and mean isopycnals are approximately orthogonal. The following sequence of events therefore can be anticipated for an inhomogeneous gas suddenly forced by a standing acoustic wave. (i) Given a smooth initial distribution of the thermal inhomogeneity, acoustic-wave vorticity will spread throughout the domain and drive a streaming flow (figure 3a). (ii) The streaming flow will mix the interior of the fluid, causing inhomogeneities in temperature to largely be confined to thin boundary layers (figures 2c and 7a). (iii) This evolution in the temperature field will modify the characteristics of the acoustic wave; in particular, acoustic vorticity, and hence the component of the acoustic force density effectively driving the streaming flow, also will localize in these boundary layers (figure 4). (iv) The streaming flow will develop viscous boundary layers close to the hot (or cold) solid boundary (figure 7b) that results in jets sustaining the mixing in the interior of the domain (figure 6b–e). We conjecture that a transition to turbulence also may be expected to occur as the acoustic amplitude or Reynolds number is further increased.
This dynamical picture highlights the two-way coupling between the acoustic wave and streaming flow. To stress this point, we carried out one-way coupled simulations in which the acoustic wave fields are not updated as the streaming temperature field evolves. The resulting statistically stationary state is no longer strictly steady and the associated heat flux is not accurately estimated; see figure 5. Holding the acoustic velocities fixed in the evaluation of the acoustic force density (1.2), either to simplify the numerical algorithm by not regularly solving an eigenvalue problem or to draw an analogy with an effective gravitational field (Koulakis & Putterman Reference Koulakis and Putterman2021; Koulakis et al. Reference Koulakis, Ofek, Pree and Putterman2023), therefore restricts the validity of the approach to waves of vanishing amplitude that do not drive any significant streaming flow and for which the advection of temperature inhomogeneities can be neglected. Indeed, figure 8, which shows the acoustic wave kinetic energy for both one-way (figure 8
$a$) and properly two-way (figure 8
$b$) coupled simulations, confirms that merely updating the mean density distribution in the acoustic force density (1.2) is not sufficient for reliable and robust simulation of these flows. Multiple scale analysis is the appropriate approach to treat such disparate time scale coupling. Finally, this multiscale analysis can be extended to identify the scaling of the temperature differential across the domain 
$\varGamma = \Delta \varTheta _*/T_*$ for which the feedback of the streaming flow on the wave must be included, i.e. to identify the crossover between pure Rayleigh streaming (
$\varGamma = 0$) and baroclinic streaming (
$\varGamma = O(1)$).
Figure 8. Evolution for 
$A=4$ and 
$\delta = 4$ of the acoustic-wave kinetic energy 
$\frac {1}{2}\overline {|\boldsymbol {u}_1'|^2}$: the initial condition (a) and the steady state (b). The evolution of the acoustic-wave kinetic energy, in addition to the streaming density field, also results in significant modifications to the acoustic force density (1.2).
Supplementary movie
Supplementary movie is available at https://doi.org/10.1017/jfm.2024.744.
Funding
This work is supported by Centre national d’études spatiales (CNES). We also acknowledge support from the Woods Hole Geophysical Fluid Dynamics Summer School (NSF 1829864), where this work was initiated.
Declaration of interests
The authors report no conflict of interest.
Open access
