<footnote n="29"><sup>29</sup>C. H. Lees. "The free periods of a composite elastic column or composite stretched wire," Proc. Phys. Soc. **41**, 204–213 (1929). Cited in H. Gué-noche, "Flame propagation in tubes and in closed vessels," in *Non-steady Flame Propagation*, edited by G. H. Markstein (Pergamon Press, New York, 1964), pp. 107–181.</footnote>

<footnote n="30"><sup>30</sup>G. Mozurkewich, "Time-average temperature distribution in a thermo-acoustic stack," J. Acoust. Soc. Am. **103**, 380–388 (1998).</footnote>

<footnote n="31"><sup>31</sup>V. Gusev, "Relaxation-time approximation for analytical evaluation of temperature field in thermoacoustic stack," J. Sound Vib. **235**, 711–726 (2000).</footnote>

<footnote n="32"><sup>32</sup>H. Ishikawa and D. J. Mee, "Numerical investigations of flow and energy fields near a thermoacoustic couple," J. Acoust. Soc. Am. **111**, 831–839 (2002).</footnote>

<footnote n="33"><sup>33</sup>A. Berson and P. Blanc-Benon, "Nonperiodicity of the flow within the gap of a thermoacoustic couple at high amplitudes," J. Acoust. Soc. Am. **122**, 122–127 (2007).</footnote>

<footnote n="34"><sup>34</sup>L. Zoontjens, C. Q. Howard, A. C. Zander, and B. S. Cazzolato, "Numerical study of flow and energy fields in thermoacoustic couples of non-zero thickness," Int. J. Therm. Sci. **48**, 733–746 (2008).</footnote>
L. Ma, C. Weisman, D. Baltean-Carlès, I. Delbende, L. Bauwens

Effect of a resistive load on the starting performance of a standing wave thermoacoustic engine : A numerical study

J. Acoust. Soc. Am. 138 (2), 847-857, 2015
# Effect of a resistive load on the starting performance of a standing wave thermoacoustic engine: A numerical study

Lin Ma, Catherine Weisman,ª Diana Baltean-Carlès, and Ivan Delbende

Université Pierre et Marie Curie UPMC Sorbonne Universités, UFR 919 - Laboratoire d'Informatique pour la Mécanique et les Sciences de l'Ingénieur, 91405 Orsay Cedex, France

Luc Bauwens

Department of Mechanical and Manufacturing Engineering, University of Calgary, Calgary, Alberta T2N 1N4, Canada

(Received 4 August 2014; revised 23 April 2015; accepted 8 July 2015; published online 17 August 2015)

The influence of a resistive load on the starting performance of a standing-wave thermoacoustic engine is investigated numerically. The model used is based upon a low Mach number assumption; it couples the two-dimensional nonlinear flow and heat exchange within the thermoacoustic active cell with one-dimensional linear acoustics in the loaded resonator. For a given engine geometry, prescribed temperatures at the heat exchangers, prescribed mean pressure, and prescribed load, results from a simulation in the time domain include the evolution of the acoustic pressure in the active cell. That signal is then analyzed, extracting growth rate and frequency of the dominant modes. For a given load, the temperature difference between the two sides is then varied; the most unstable mode is identified and so is the corresponding critical temperature ratio between heater and cooler. Next, varying the load, a stability diagram is obtained, potentially with a predictive value. Results are compared with those derived from Rott's linear theory as well as with experimental results found in the literature. © 2015 Acoustical Society of America. [http://dx.doi.org/10.1121/1.4927092]

[JDM]

Pages: 847–857

## I. INTRODUCTION

A thermoacoustic engine is a device absorbing heat at the hot heat exchanger and releasing heat at the cold heat exchanger while producing acoustic work as an output. Engines are meant to carry a load, which in the case of an acoustic device might include a thermoacoustic refrigerator,¹ a pulse tube refrigerator,² a (possibly linear) electric generator,³ or a piezoelectric transducer.⁴ The strong coupling that will exist between the engine and the load influences both starting performance and steady (periodic) operation.

Thermoacoustic devices have been studied extensively for more than 30 years. Most of the work in the literature is based upon the linear theory developed by Rott⁵ and later extended by Swift *et al.*¹,⁶ The theory can handle different geometries for the various parts of the system, such as stack, heat exchangers, and ducts. Numerical implementation is based upon a lumped element approximation.⁷ The linear theory will then yield the acoustic field and, as a result, the performance of a given system when in the periodic regime.¹ Linear theory can also predict the onset conditions (critical starting temperature gradient) as a function of the mean pressure by performing a stability analysis for a fixed geometry and working fluid. Results from linear theory are in good agreement with experiments⁸⁻¹⁰ for the onset of a standing wave thermoacoustic engine inside a resonator closed at both ends. Recently, a completely analytical boundary layer theory extending linear theory was developed by Sugimoto,¹¹,¹²

reducing the problem to a one-dimensional model. That theory was then used in an analytical stability analysis for systems with a stack under given temperature gradient inside straight or looped resonators, for given acoustic frequency, without resorting to lumped elements. Linear theory has important limitations: the stability analysis is performed in the frequency domain for each mode, the transient regime is not described, and nonlinear and multidimensional effects are not considered or are only partially taken into account.

Transient regimes were addressed by de Waele¹³ who developed a one-dimensional time-dependent model of a thermoacoustic Stirling engine taking into account the dynamics of individual components, again using a lumped-element approximation. de Waele studied engine start phase and transient effects while acoustics is still considered linear. A mixed empirical and theoretical one-dimensional model was developed by Guedra and co-workers¹⁴,¹⁵ based upon a transmission matrix formulation, also assuming linear acoustics. The transmission matrix was obtained experimentally from impedance measurements. Using that model, they studied the onset of oscillations for both a basic standing wave engine and a traveling wave-engine. None of these studies included the influence of a load on onset conditions.

In contrast, the influence of the load has been accounted for in studies of thermoacoustic devices in the periodic regime. The performance of a loaded thermoacoustic engine was investigated by Olson *et al.*,¹⁶ and also others,¹⁷,¹⁸ using linear theory¹ and the lumped elements approximation. More recently, Hatori *et al.*¹⁹ developed an experimental method based on measuring the acoustic impedance of the engine and of the load separately, and predicted the operating point

<footnote>ªElectronic mail: catherine.weisman@limsi.fr</footnote>
of the combination engine plus load. Coming to nonlinear models, some numerical simulations of the full compressible Navier-Stokes equations in loaded thermoacoustic devices with complicated geometry have been performed.²⁰,²¹ These simulations were performed using commercial codes and gave a good description of thermoacoustic engines in the transient and periodic regimes. However, each test case required extremely large and time consuming simulations; therefore, these are not appropriate for a parametric study of a variable load coupled with an engine.

The aim of the current work is to study the influence of a variable load on the marginal stability conditions of a standing-wave thermoacoustic engine. The model used implements a hybrid low Mach number approximation. It is based on coupling the two-dimensional nonlinear flow and heat exchange in the thermoacoustic active cell with one-dimensional linear acoustics in the loaded resonator. For a detailed description, see Ref. 22. In the load model used in the present study, instantaneous acoustic pressure and velocity are assumed to be in phase at the loaded end, and the load is thus reduced to a real scalar value measuring the ratio of velocity over acoustic pressure. For a given engine geometry and a given temperature ratio between heater and cooler, numerical simulations of the active cell show that depending on the load value, the engine will or will not start. From the numerical results, for any given configuration, a critical temperature ratio is thus determined, corresponding to the onset of thermoacoustic oscillations, together with the frequency of the corresponding, most unstable mode. The analysis is performed using a time-dependent nonlinear model that takes into account two-dimensional effects as well as multi-frequency content. Finally, the influence of the load on the starting phase of the engine is discussed and the results are compared, both with results derived from the linear theory and with Atchley's experiments,⁸ bringing new insight into the engine-load coupling.

## II. PHYSICAL MODEL

### A. Low Mach number approximation

The geometry consists of a long tube with length $\tilde{L}_{res}$, within which an active cell is placed at a distance $\tilde{L}_L$ from the left end. The active cell consists of a stack of horizontal solid plates with length $\tilde{L}_{stack}$, placed between two heat exchangers also consisting of solid horizontal plates that have the same periodicity as the stack (see Fig. 1). The geometry can then be reduced to a domain consisting of two half-plates plus the gap between them, plus a region that extends away from the solid plates.

First, a dimensionless coordinate system $\hat{x}$ is introduced, that resolves the acoustic aspect of the device, which is thus scaled by the resonator length $\tilde{L}_{res}$. The inert resonator end located at $\hat{x} = -\hat{l}_L = -\tilde{L}_L/\tilde{L}_{res}$ is a closed end. The second resonator end consists of a load placed at a fixed location $\hat{x} = \hat{l}_R = 1 - \hat{l}_L$, described in more detail below. The active cell is positioned at $\hat{x} = 0$. It has negligible size when scaled by $\tilde{L}_{res}$, as discussed below [see Fig. 2(a)].

Next, a different dimensionless coordinate system $(x, y)$ is introduced, that characterizes the active cell geometry,

<figure><img src="image_81.png" /><figcaption>FIG. 1. Simplified geometry of a thermoacoustic engine.</figcaption></figure>

where lengths are scaled by the stack length $\tilde{L}_{stack}$, in both $x$ and $y$. The corresponding dimensionless geometrical parameters include the vertical gap or flow passage between two stack plates $h$, the stack porosity $h/H$, where $H$ is the domain height, equal to the vertical gap plus the thickness of one stack plate, the vertical gap between two heat exchanger plates $h_x$, the horizontal gap between heat-exchangers and stack $L_{xs}$, and the heat-exchanger lengths $L_{hx}$ and $L_{cx}$ [see Fig. 2(b)].

The low Mach number model is described in detail in previous work.²² The model assumes that the particle displacement spans a length on the order of $\tilde{L}_{stack}$, hence typical velocities on the order of $\tilde{U}_{ref} = \tilde{L}_{stack}/\tilde{\tau}$, where $\tilde{\tau}$ is on the order of the fundamental acoustic period of the empty resonator, $\tilde{\tau} = \tilde{L}_{res}/\tilde{c}$, and $\tilde{c}$ is the speed of sound. A reference Mach number $M = \tilde{U}_{ref}/\tilde{c}$ can then be introduced which also equals $M = \tilde{L}_{stack}/\tilde{L}_{res}$. In the following, the Mach number is assumed to be small or equivalently the stack is then taken to be much shorter than the resonator. The active cell size is then negligible in the $\hat{x}$ coordinate system.

Flow in the resonator is then described by linear acoustics. The solutions are expressed as a pair of planar traveling waves that move respectively left and right at the speed of sound. The two Riemann invariants $\mathcal{L}(\hat{x}, t)$ and $\mathcal{R}(\hat{x}, t)$ are used, defined by

$$ \mathcal{L} = \gamma u - \sqrt{T} p', \quad \mathcal{R} = \gamma u + \sqrt{T} p' \quad (1) $$

with $\gamma$ the ratio of specific heats, $u(\hat{x}, t)$ the dimensionless gas velocity, $T$ the dimensionless temperature (using the cold temperature as a reference) which is assumed to be constant in each resonator part ($T=T_h$ on the left side of the

<figure><img src="image_82.png" /><figcaption>FIG. 2. Geometry of resonator (a) and of active cell geometry (b) in dimensionless coordinates.</figcaption></figure>
active cell, and $T = T_c = 1$ on the right side), and $p'(\hat{x}, t)$ the dimensionless order $M$ (acoustic) pressure (using the mean gas pressure as a reference). At each tube end, the boundary condition results in a relationship that determines the outgoing wave as a function of the incoming wave, as explained in Ref. 22. The order $M$ pressure and acoustic velocity on both sides of the active cell ($\hat{x} = 0^{\pm}$) are expressed using these same quantities at previous times accounting for propagation on each side and reflections at the respective ends. The specific treatment of the loaded end will be described in detail in Sec. II B.

Flow in the active cell is described by two-dimensional Navier–Stokes equations, with order $M^2$ pressure gradients superimposed to a pressure spatially uniform up to order $M$, and varying density and temperature. The solid stack plate problem is described by the heat conduction equation. In the stack and heat exchanger plates, physical properties such as density, thermal conductivity, and heat capacity are taken to be constant. Thermal conductivity and heat capacity of the gas are also constant, with values corresponding to cold temperature, chosen as the reference temperature. Reference density of the gas is that corresponding to the reference (cold) temperature and the mean gas pressure.

At the solid-gas boundaries, continuity of temperature and heat flux are imposed, as well as a no slip condition. In the heat exchangers, uniform, constant temperatures are imposed with dimensionless values $T_c=1$ in the cooler plates, $T_h > 1$ in the heater plates. Given the periodicity, open active cell boundaries are effectively adiabatic.

Boundary conditions at the open ends of the active cell are provided by appropriate matching of the inner (active cell) and outer (resonator) solutions. The latter are expressed using the Riemann invariants $\mathcal{L}$ and $\mathcal{R}$ defined by Eq. (1), traveling respectively left and right in the resonator, at the speed of sound. At each extremity of the active cell, the outgoing invariant is determined as a function of the incoming one, taking the corresponding resonator boundary conditions into account. These consist respectively of a zero velocity at the closed end, and the load, described in Sec. II B, at the loaded end.²² Finally, a global energy balance over the active cell relates velocities at the two cell extremities to heat transfer.

From the standpoint of resonator acoustics, the active cell is transparent to pressure but provides a source of volume. From the standpoint of the active cell, matching with the resonator occurs at infinity so that, in principle, the active cell has infinite length on both sides. However the numerical domain has to be finite. Its length is chosen equal to 5 times the stack length, which ensures proper matching, as confirmed by validation. The dimensionless parameters characterizing the problem are then the reference Mach number $M$, the relative position of the active cell in the resonator $\hat{l}_L$, parameters describing the geometry of the active cell ($H, h_s, h_x, L_{xs}, L_{hx}, L_{cx}$), a Péclet number $P_e = \tilde{U}_{ref}\tilde{H}/\tilde{\alpha}_{ref}$, where $\tilde{\alpha}_{ref}$ is the reference thermal diffusivity, a Reynolds number $R_e = \tilde{\rho}_{ref}\tilde{U}_{ref}\tilde{H}/\tilde{\mu}_{ref}$, $\tilde{\rho}_{ref}$ and $\tilde{\mu}_{ref}$ being the reference density and viscosity, and the Péclet number $P_{es} = P_e \tilde{\alpha}_{ref}/\tilde{\alpha}_s$ in solid walls, where $\tilde{\alpha}_s$ is the thermal diffusivity of the solid material.

## B. Load model

The present study considers a purely resistive load, so that instantaneous acoustic pressure $p'(\hat{l}_R, t)$ and velocity $u(\hat{l}_R, t)$ at the loaded end are assumed to be in phase. The load is thus characterized by a scalar value equal to the ratio of dimensionless acoustic pressure to velocity, i.e., $f = p'(\hat{l}_R, t)/u(\hat{l}_R, t)$, with $f$ a positive real. Such a load behaves as a pure resistance in the electrical analogy. Actual loads may include inertial and capacitance effects; these will depend not only on mechanical design, but also on the driven devices. A detailed study of these issues is beyond the scope of the present work.

The limit value $f \to 0$ corresponds to an open end, with $p'(\hat{l}_R, t) = 0$, while the limit value $f \to \infty$ corresponds to a closed end, with $u(\hat{l}_R, t) = 0$.

The boundary condition at the load location $\hat{x} = \hat{l}_R$ is carried along the returning (left-moving) characteristic, resulting in the corresponding boundary condition at the right end of the active cell, at $\hat{x} = 0^+$. Specifically, the boundary condition at $\hat{x} = 0^+$ relates the current local value of the left-moving Riemann invariant and the value of the right-moving Riemann invariant at the same location but at a previous time that differs by the round trip time between the active cell and the load at the speed of sound $\sqrt{T_c}$,

$$ \mathcal{L}(0^+, t) = Z\mathcal{R}(0^+, t - 2\hat{l}_R / \sqrt{T_c}), \quad (2) $$

in which a real coefficient $Z$ was introduced, defined as

$$ Z = \frac{\gamma - f}{\gamma + f}. \quad (3) $$

While $f$ goes from zero for an open end to infinity for a closed end, $Z$ has the finite span from $-1$ for a closed end to $+1$ for an open end. Furthermore, the expression emphasizes the particular value $f = \gamma$, for which $Z$ vanishes, corresponding to a left-moving Riemann variable with zero value, regardless of the right-moving variable, hence to a non-reflecting load. If $Z=0$, the solution no longer depends upon the length of the right side of the resonator; such a load may be difficult to implement, and of little practical value, as results below will confirm; however, that case still has some interest as it separates different regimes. In the framework of the harmonic approximation, it is possible to relate $f$ or $Z$ to the reflection coefficient at the loaded end hence to the standing wave ratio.²³ In that context the absolute value of $Z$ is equal to the modulus of the reflection coefficient. The limit values $Z = \pm 1$ correspond to a pure standing wave in the resonator, while the limit case $Z=0$ corresponds to a pure traveling wave moving from the active cell towards the load. Between those limit cases, the wave is a combination of two waves traveling in opposite directions with unequal amplitudes, i.e., neither standing or purely traveling. In the current study, however, the effect of the load parameter on starting performance will be investigated in the time domain with no a priori on the nature of the wave developing in the engine. Critical temperature ratio between heater and cooler and starting frequencies will be analyzed by varying the value of
f in the range 0 to +∞. Results will be more conveniently shown as functions of the coefficient Z in the range [-1, 1].

### III. NUMERICAL SIMULATION

#### A. Algorithm

The problem in the active cell is solved numerically using a finite volume code.²² Diffusion is dealt with implicitly while advection is explicit. The scheme is second-order accurate in space and in time. A fractional time step projection method adapted for variable density is used to enforce continuity. The ADI algorithm is used to solve the Helmholtz equation for time advance of temperature and velocities. A multigrid algorithm determines the pressure correction. Solutions of the coupled equations providing velocity boundary conditions and acoustic pressure at the active cell location from resonator acoustics and energy balance are appropriately integrated in the solution sequence.²²

#### B. Procedure

Simulations account for a specific engine geometry, characterized by the reference Mach number M, the relative distance of the active cell to the left end $\hat{L}$, and the dimensions of the active cell H, h, $h_x$, $L_x$, and $L_{xs}$. Most of the current results were obtained for $L_{xs}$ equal to h. Reference Péclet number $P_e$, solid Péclet number $P_{es}$, and Reynolds number $R_e$ are maintained constant. The evolution in time of the velocity and temperature fields is computed for given heater and cooler temperatures and load. The required initial conditions consider fluid at rest and temperature profiles, obtained numerically, corresponding to steady conduction in the walls and in fluid at rest. Initial conditions also need to be specified for resonator acoustics, and these need to be nonzero since the steady solution to the global problem consists of fluid at rest, together with a stationary conduction solution. To that effect, and in order not to favor specific modes, an initial random noise is imposed in the resonators. Small $O(10^{-4})$ random values are thus assigned to the Riemann invariants at the active cell entrance and exit for a time window prior to zero, of length equal to the round trip time between active cell and the respective resonator end.

Starting from initial conditions specified above, the velocity $u(x, y, t)$, dynamic pressure $p(x, y, t)$, temperature $T(x, y, t)$, and density $\rho(x, y, t)$ fields are calculated inside the active cell at each time step. Coupling at the active cell entrance and exit provides the velocity at the entrance $u(0^-, t)$ and exit $u(0^+, t)$ of the active cell, as well as the value of the acoustic pressure at the active cell location $p'(0, t)$ at each time step. The velocity values at the entrance and exit of the active cell serve as boundary conditions for velocity inside the active cell. The order M acoustic pressure $p'(0, t)$, which is included in the active cell model only through boundary conditions, is an outcome of the simulation. That quantity is analyzed here to characterize the device's regime: motion damping or amplification.

#### C. Post-processing

Acoustic modes are characterized by amplitude and complex frequency, i.e., frequency and growth rate. If the temperature ratio between heater and cooler is small, growth rates are all negative, i.e., the amplitudes of all modes decrease in time and motion stops. A critical temperature ratio is obtained when one of the modes exhibits positive growth: instability occurs and the device starts. In order to find the critical temperature ratio $T_{crit} = \tilde{T}_{hcrit}/\tilde{T}_c$ for a given load, repeated simulations are performed for increasing values of the hot temperature, until a first growth rate becomes positive.

The pressure signal $p'(0, t)$ is analyzed to extract the growth rate and frequency of the dominant modes, by fitting the signal to an approximation of the form

$$p'(0, t) = \sum_{i=1}^{n} A_i e^{\sigma_i t} \cos(\omega_i t) + B_i e^{\sigma_i t} \sin(\omega_i t), \quad (4)$$

where n is the desired number of modes, with typical value 4. The fitting uses a window that includes m consecutive data points, collected after several reference periods have elapsed so that a dominant exponential mode has emerged from the initial phase affected by initial conditions. In Eq. (4), $\sigma_i$ is the growth rate of mode i, $\omega_i$ is its angular frequency, and $A_i, B_i$ are amplitudes.

Fitting the ($A_i, B_i, \sigma_i, \omega_i$) coefficients is performed using an iterative two stage procedure. Starting with an initial guess ($\sigma_i, \omega_i$), the ($A_i, B_i$) coefficients are determined by least mean square error minimization, using the m data points. Then the ($\sigma_i, \omega_i$) coefficients are obtained using a simplex algorithm. Iterations continue until the fit satisfies a given criterion, with the quadratic error below some adequate lower bound. The whole procedure is iterated until an error below $10^{-4}$ on ($\sigma_i, \omega_i$) is reached.

### IV. RESULTS

#### A. Conditions

Numerical simulations were carried out on an existing prototype of thermoacoustic engine,²² with a “short” thermoacoustic active cell inserted in a long resonator closed at both ends. The engine is modeled as the simplified device introduced in Fig. 2, assuming acoustics is lossless in the resonators and concentrating all dissipation at the loaded end.

The working fluid is helium. The dimensions (resonator length, stack length, distance between the left end and the active cell, plate spacing) and the conditions of the experiments, such as mean pressure and cold temperature, are kept constant. The relevant dimensionless parameters are summarized in Table I, and the relevant fluid properties are summarized in Table II.

Numerical simulations were carried out using a regular two-dimensional mesh for the active cell, with either 2048 × 32 or 4096 × 64 grid points. The simulation domain is long enough to ensure that the flow velocity near the entrance and exit is parallel to the resonators. The time step is adjusted in order to satisfy the numerical stability
<table><caption>TABLE I. Dimensionless parameters characterizing the reference experimental device (Ref. 22).</caption><thead><tr><th>M</th><th>P<sub>e</sub></th><th>P<sub>es</sub></th><th>R<sub>e</sub></th><th>l&#x0302;<sub>L</sub></th><th>H</th><th>h/H</th><th>h<sub>s</sub>/H</th><th>L<sub>hx</sub></th><th>L<sub>ex</sub></th><th>L<sub>xx</sub></th></tr></thead><tbody><tr><td>0.02</td><td>1100</td><td>4900</td><td>1600</td><td>0.09</td><td>0.007</td><td>0.73</td><td>0.73</td><td>0.05</td><td>0.05</td><td>0.005</td></tr></tbody></table>

condition, with about 200 (coarse grid) to 2000 (fine grid) time steps per reference acoustic period. Typically only 60 min of CPU time are necessary for each numerical run (about 100 periods), such as the one shown in Fig. 3.

For given load value and given heater temperature, the time history of acoustic pressure in the active cell was recorded from the computation. An example of the time signal is shown in Fig. 3(a). The first oscillations are rather erratic [details are not shown in Fig. 3(a)] and they depend upon initial conditions. The signal eventually becomes smoother and amplitudes appear to correspond to a superimposition of exponentially varying modes. The signal appears to remain multi-frequency. This is confirmed in Fig. 3(b) which shows the corresponding spectrogram. In the figure, the angular frequency is scaled by a reference angular frequency $\omega_{ref}$ equal to the fundamental angular frequency of the empty resonator closed at both ends, filled with gas at a mean pressure $\bar{P}_m$ and a temperature equal to the cold temperature $\bar{T}_c$. At any given time, the power spectrum exhibits successive peaks corresponding to the resonant frequencies of the engine for the specific load value. As time evolves, most modes are damped while some are amplified. A more accurate determination of the frequencies and growth rates of the dominant modes is presented in Sec. IV B.

## B. Results at a given active cell location

First, the effect of the load on the critical temperature ratio is studied, for a single active cell location.

Analyzing the pressure signal using the approach described in Sec. III C, the growth rate and frequency of the four dominant modes are calculated, for a range of increasing values of the temperature of the heater $T_h$, for a given load $f$, until unstable modes appear. The most unstable mode is identified as the first mode for which the growth rate $\sigma$ goes from negative to positive, as $T_h$ is increased. Plotting $\sigma$ as a function of the dimensionless temperature $T_h$ helps one to identify the precise value of the critical temperature ratio, as Fig. 4 shows. In the vicinity of the critical value, $\sigma$ varies nearly linearly with $T_h$, and a precise determination of the critical temperature ratio is obtained by interpolation. For

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

<figure><img src="image_84.png" /><figcaption>FIG. 3. (Color online) (a) Typical time signal of acoustic pressure $p'(0, t)$ and (b) associated spectrogram (the magnitude is the square root of the power density spectrum). $Z = 0.988$. $T_h = 1.08$.</figcaption></figure>

example in Fig. 4, the critical temperature ratio is found to be $T_{crit} = 1.845$.

Next, varying the load and reproducing the same procedure for each load value results in a stability diagram (Fig. 5) in the $(Z, T_h)$ plane showing stable and unstable zones. The zones are separated by a stability curve connecting all the critical temperature ratios causing the engine to start. The region above the stability curve corresponds to instability. Figure 5 shows a stability domain roughly delimited by a symmetrical triangle peaking at $Z=0$, when a semi-log scale is used.

<figure><img src="image_85.png" /><figcaption>FIG. 4. Growth rate of most unstable mode as a function of the hot exchanger temperature, $Z = -0.963$.</figcaption></figure>

<table><caption>TABLE II. Fluid properties in reference experiment (Ref. 22).</caption><thead><tr><th colspan="8">Helium gas</th></tr><tr><th>P&#x0302;<sub>m</sub><br>MPa</th><th>T&#x0302;<sub>c</sub><br>K</th><th>ρ&#x0302;<br>kg m<sup>-3</sup></th><th>μ&#x0302;<br>Pa s</th><th>ζ&#x0302;<sub>p</sub><br>J kg<sup>-1</sup> K<sup>-1</sup></th><th>κ&#x0302;<br>W m<sup>-1</sup> K<sup>-1</sup></th></tr></thead><tbody><tr><td>1.0</td><td>293</td><td>1.6</td><td>2 × 10<sup>-5</sup></td><td>5193</td><td>0.15</td></tr></tbody></table>
<figure><img src="image_86.png" /><figcaption>FIG. 5. Stability diagram in the (Z, T_h) plane, showing the stability curve (●-●). The region above the curve corresponds to instability.</figcaption></figure>

Figure 6 shows the angular frequency ω of the first mode to become unstable as a function of Z. As |Z| decreases, the frequency of the most unstable mode is observed to increase. Although the general shape of the curve is approximately symmetrical as it is in Fig. 5 for the critical temperature, the values of ω do not display any symmetry. Indeed, critical modes near Z = -1 are of the closed-end type (multiples of λ/2) whereas critical modes near Z = 1 are of the open-end type (odd multiples of λ/4). This results in distinct sets of onset frequencies, as is further discussed and analyzed in Secs. IVC and V when commenting on Figs. 7, 8, and 11.

For Z close to zero, Fig. 5 shows a very large critical temperature ratio, too large to be realistic. Likewise, frequencies become so high that resolution, both of the wavelength in space and of the period in time, becomes quite poor. Furthermore, these curves were obtained assuming constant viscosity and thermal conductivity conditions, which is unrealistic at high temperature. These results are thus somewhat dubious and, therefore, Fig. 5 merely gives a qualitative picture for Z close to zero. While such numerical critical temperature ratios may be of limited value, the conclusion is valid that in the vicinity of Z=0 the device will not start for any realistic temperature ratio (T_h less than about 10), all modes remaining stable.

<figure><img src="image_87.png" /><figcaption>FIG. 6. Angular frequency of the most unstable mode scaled by ω ref, as a function of Z.</figcaption></figure>

<figure><img src="image_88.png" /><figcaption>FIG. 7. Growth rates as a function of the active cell location for 0 < $\hat{I}_L$ < 0.5, Z = -0.963, T_h = 1.2. (●-●) fundamental mode (the λ/2 mode), (—*) first harmonic mode (the λ mode).</figcaption></figure>

### C. Results varying the active cell location or configuration

The stability diagrams shown in Figs. 5 and 6 can be used, for a given configuration, to predict the critical temperature ratio and the dominant frequency as a function of the load.

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

<figure><img src="image_90.png" /><figcaption>FIG. 8. Stability diagram showing the critical temperature ratio as a function of the active cell location for (a) (□-□) Z = -0.99, (-*-) Z = -0.92, (●-●) Z = -0.85, the critical mode is the fundamental mode (λ/2), and (b) (●-●) Z = 0.79, (-*-) Z = 0.90, (□-□) Z = 0.99, the critical mode is the first harmonic (3λ/4).</figcaption></figure>
These diagrams are modified if the location of the active cell measured by $\hat{l}_L$ is varied. The stability curve becomes a stability surface, which can be used to predict the critical temperature ratio for a given setup. For given load and heater temperature values, the position of the active cell has to be within a certain region for the engine to start. For an active cell outside this region, all modes are dampened and the engine will not start.

Figure 7 shows the growth rate of the fundamental mode and of the first harmonic as a function of the active cell location for $Z = -0.963$ and $T_h = 1.2$. It is found that if the active cell is positioned at $0.083 < \hat{l}_L < 0.12$, the first harmonic is the most unstable mode. In contrast, if it is placed at $0.12 < \hat{l}_L < 0.43$, the fundamental mode is dominant. For that parameter setting, the overall optimal position for starting performance is $\hat{l}_L = 0.27$, for which the growth rate is maximum. Another way to address the issue is to determine the active cell location for which the critical temperature ratio is minimum. Figure 8 shows the value of the critical temperature ratio for varying active cell locations and for selected load values. Figure 8(a) shows results for negative values of $Z$ near $Z = -1$, in which case both ends are closed. For each value of $Z$, there is one stack position that minimizes the critical temperature ratio. The optimal position is approximately 0.25 for $Z = -0.99$. Since the fundamental mode (the $\lambda/2$ mode) is the critical one in that region, this corresponds to a stack located at $\lambda/8 \approx \hat{L}_{res}/4$ as expected ($\hat{L}_{res} = 1$). As $Z$ departs from $Z = -1$, the wavelength increases so that the optimal position shifts towards the center of the resonator.

Figure 8(b) corresponds to positive values of $Z$, near $Z=1$, for which the cold end is open. The optimal position is approximately 0.16 for $Z=0.99$. Now the first harmonic (the $3\lambda/4$ mode) is the critical mode in that range of stack locations. This corresponds to a stack located at $\lambda/8 \approx \hat{L}_{res}/6 \approx 0.17$. For decreasing values of $Z$, below $Z=1$, the wavelength increases so that the optimal position again shifts towards the resonator center. Note that here the active cell is located at $\hat{l}_L \le 0.25$. If the stack were placed further away, the fundamental mode would become the critical one and the optimal position would shift towards 0.5.

The effect of varying the distance between the heat exchangers and the stack $L_{xs}$ was also studied. In all previous cases the horizontal spacing $L_{xs}$ was equal to the vertical stack plate spacing $h$. Figure 9 shows the stability diagram for two other situations: for heat-exchangers and stack in direct contact ($L_{xs}=0$) and doubling the horizontal spacing ($L_{xs} = 2h$). While placing the heat-exchangers and stack incontact together results in significant reduction of the critical temperature ratio for all values of the load (potentially leading to a significant increase in longitudinal conduction losses), the effect of doubling the spacing is insignificant. The nature of the critical modes will be discussed in Sec. V (see comments of Fig. 11).

## V. COMPARISON WITH THE LINEAR THEORY

In this section the critical temperature ratio associated to a given acoustic mode is estimated following Rott's linear theory,⁵,²⁴ which is outlined here using dimensional values

<figure><img src="image_91.png" /><figcaption>FIG. 9. Stability diagram in the $(Z, T_h)$ plane, showing the stability curve (critical temperature ratios) for three distances between heat-exchangers and stack, $(-\Box-) L_{xs} = 0$, $(-\bullet-) L_{xs} = h$, $(-\ast-) L_{xs} = 2h$. (a) $-1 < Z < -0.6$, (b) $0.6 < Z < 1$.</figcaption></figure>

for the various parameters. Linear results were obtained for a slightly different configuration, with the same total resonator length, stack length, and position of the stack center within the resonator as above. However, the detailed heat exchanger geometry used above was replaced by imposing a linear wall temperature profile in the stack, constant mean temperatures $\tilde{T}_h$ on the left side of the stack, and constant mean cold temperature $\tilde{T}_c$ on the right side.

The flow in the resonator tubes left and right of the stack was taken as one-dimensional isentropic linear acoustics with different speeds of sound corresponding to the respective temperatures $\tilde{T}_h$ and $\tilde{T}_c$. The boundary conditions on the left and right ends were the same as previously, with a load modeled as a real impedance $\tilde{R}$ relating acoustic velocity and acoustic pressure. The acoustic pressure and velocity amplitudes $[\tilde{p}'_h(\tilde{x}), \tilde{u}_h(\tilde{x})]$ on the left resonator part and $[\tilde{p}'_c(\tilde{x}), \tilde{u}_c(\tilde{x})]$ on the right resonator part are functions of position which can be expressed simply in terms of two unknowns, namely, the complex maximum amplitudes $\tilde{U}_h$ and $\tilde{U}_c$.

The stack region is described as a two-dimensional channel of height $2\tilde{y}_0$ equal to the vertical distance between stack plates, extending from $\tilde{x} = 0$ to $\tilde{x} = \tilde{L}_{stack}$. Following Rott's approach, solutions of the viscous flow equations are sought for as harmonic functions in time with angular frequency $\tilde{\omega}$,
using a boundary layer approximation and integrating vertically over the channel height. Also, the stack is assumed to be acoustically compact, so that the acoustic pressure amplitude $\tilde{p}'$ is homogeneous through the stack. Since temperature varies linearly through the stack, $d\tilde{T}/d\tilde{x}$ is a constant.

The resulting dimensional equation for the acoustic velocity in Fourier space is²⁵

$$ \frac{d\tilde{u}_1}{d\tilde{x}} - \frac{\tilde{f}_\kappa - \tilde{f}_\nu}{(1-\tilde{f}_\nu)(1-P_r)} \frac{1}{\tilde{T}} \frac{d\tilde{T}}{d\tilde{x}} \tilde{u}_1 = -\frac{i\tilde{\omega}}{\gamma \tilde{P}_m} [1+(\gamma-1)\tilde{f}_\kappa]\tilde{p}', \quad (5) $$

where $\tilde{P}_m$ is the gas mean pressure, Rott's factors are defined as

$$ \begin{aligned} \tilde{f}_{\nu} &= \frac{\tanh\left[\frac{(1+i)\tilde{y}_0}{\sqrt{2\tilde{\mu}/(\tilde{\rho}\tilde{\omega})}}\right]}{(1+i)\tilde{y}_0/\sqrt{2\tilde{\mu}/(\tilde{\rho}\tilde{\omega})}}, \\ \tilde{f}_{\kappa} &= \frac{\tanh\left[\frac{(1+i)\tilde{y}_0}{\sqrt{2\tilde{k}/(\tilde{\rho}\tilde{c}_p\tilde{\omega})}}\right]}{(1+i)\tilde{y}_0/\sqrt{2\tilde{k}/(\tilde{\rho}\tilde{c}_p\tilde{\omega})}} \end{aligned} \quad (6) $$

with constant gas properties $\tilde{\mu}$, $\tilde{c}_p$, $\tilde{k}$ (therefore the Prandtl number $P_r$ is constant also), and density $\tilde{\rho}$ depending on position through variation of temperature. The $\tilde{f}_\nu$ and $\tilde{f}_\kappa$ factors are functions of $\tilde{x}$. Equation (5) is a first order linear inhomogeneous equation, so the solution is expressed as

$$ \tilde{u}_1(\tilde{x}) = \exp \left[ \int_0^{\tilde{x}} \frac{\tilde{f}_\kappa(\tilde{x}_1) - \tilde{f}_\nu(\tilde{x}_1)}{(1-\tilde{f}_\nu(\tilde{x}_1))(1-P_r) \tilde{T}(\tilde{x}_1)} \frac{1}{d\tilde{x}_1} d\tilde{x}_1 \right] \\ \times \left\{ C_1 - \int_0^{\tilde{x}} \frac{i\tilde{\omega}}{\gamma \tilde{P}_m} [1+(\gamma-1)\tilde{f}_\kappa(\tilde{x}_1)] \tilde{p}' \exp \left[ -\int_0^{\tilde{x}_1} \frac{\tilde{f}_\kappa(\tilde{x}_2) - \tilde{f}_\nu(\tilde{x}_2)}{[1-\tilde{f}_\nu(\tilde{x}_2)](1-P_r) \tilde{T}(\tilde{x}_2)} \frac{1}{d\tilde{x}_2} d\tilde{x}_2 \right] d\tilde{x}_1 \right\}, \quad (7) $$

in which $C_1$ is an integration constant.

That expression can be rewritten using the relationship between position and temperature as

$$ \tilde{u}_1(\tilde{x}) = \exp \left[ - \int_{\tilde{T}(\tilde{x})}^{\tilde{T}_h} \frac{\tilde{f}_\kappa(\theta) - \tilde{f}_\nu(\theta)}{[1 - \tilde{f}_\nu(\theta)](1 - P_r)\theta} d\theta \right] \\ \times \left\{ C_1 + \int_{\tilde{T}(\tilde{x})}^{\tilde{T}_h} \frac{i\tilde{\omega}}{\gamma \tilde{P}_m} [1 + (\gamma - 1)\tilde{f}_\kappa(\theta)] \tilde{p}' \right. \\ \times \exp \left[ \int_\theta^{\tilde{T}_h} \frac{\tilde{f}_\kappa(\theta^*) - \tilde{f}_\nu(\theta^*)}{[1 - \tilde{f}_\nu(\theta^*)](1 - P_r)\theta^*} d\theta^* \right] \frac{\tilde{L}_{stack}}{\tilde{T}_c - \tilde{T}_h} d\theta \right\}. \quad (8) $$

Boundary conditions are provided by matching with resonator acoustics. While the resonator acoustics yield velocities uniform transversally, the boundary layer model incorporated in the linear theory used in the stack results in a specific nonuniform velocity profile. Matching the mean velocities values is justified based upon a multiple length approach similar to that used above. Indeed the transition region is small compared with acoustic lengths, so that mass and pressure corrections due to the transition zone occur at higher order. For known boundary conditions at the outer resonator ends, acoustics yield a relationship between acoustic pressure in the stack, $\tilde{p}'$, and the uniform acoustic velocities, $\tilde{U}_h$ and $\tilde{U}_c$, respectively, on the left and right side. Continuity conditions require mean velocities and acoustic pressure to match at both stack ends.

Requiring the resulting globally homogeneous problem to have non-trivial solutions results in a dispersion relation that cannot be solved in closed form for the complex

frequency $\omega$. A numerical solution was thus implemented using a shooting method, integrating over the stack length, and iterating until the relationship between pressure and velocity at the end match the acoustic boundary condition. The method requires an initial guess for the complex frequency, conveniently provided by the numerical results of Sec. IV B.

Figures 10 and 11 show a comparison between numerical results and results obtained using the linear theory. Only the range $|Z| > 0.6$ is shown, which contains all cases of practical interest.

In Fig. 10, the critical temperature ratio is shown as a function of Z. There is very good agreement between the two approaches, with the ratio obtained through the linear theory slightly above the numerical values, which may be ascribed to a shorter length for lossless acoustics in the linear approach. As $Z \to -1$, the critical temperature ratio becomes very close to 1. For instance, for $Z = -0.997$, the critical temperature ratio is found to be $T_{crit} = 1.024$ numerically while linear theory gives $T_{crit} = 1.14$. Similarly, for an almost open-end resonator, $Z = 0.998$, numerical simulations yield a value $T_{crit} = 1.025$ vs $T_{crit} = 1.13$ from linear theory.

In Fig. 11, the onset frequency normalized by $\omega_{ref}$ is shown as a function of Z for both approaches, again showing excellent agreement. Also plotted (continuous lines) are the resonant frequencies of a lossless resonator with a temperature discontinuity positioned at $\hat{x} = 0$, for hot temperature $T_h$ set to the critical value $T_{crit}(Z)$ determined by the simulation, for a closed right end in Fig. 11(a) or an open right end in Fig. 11(b). The small frequency increase as $T_h$ increases is barely noticeable on the continuous (nearly horizontal) lines. The load itself has a small effect since the points from both
<figure><img src="image_92.png" /></figure>

<figure><img src="image_93.png" /><figcaption>FIG. 10. Stability diagrams in the (Z, T<sub>h</sub>) plane showing critical temperature ratios, (•) numerical simulation, (— —) linear theory. (a) −1 < Z < −0.6; (b) 0.6 < Z < 1.</figcaption></figure>

numerical results and linear theory, and the points on the
horizontal lines, which only take into consideration the
changes in temperature due to Z, but not the values of Z itself
only differ slightly. However, the most visible effect of the
load is the switch in the first mode to become unstable (and
thus of the onset frequency), represented by vertical lines in
Fig. 11. In Fig. 11(a), each vertical range corresponds to a
multiple of λ/2, and in Fig. 11(b), each vertical range corre-
sponds to an odd multiple of λ/4. As |Z| diminishes, progres-
sively higher frequency modes are the first to become
unstable. For the current configuration, |Z| has to be
extremely close to 1 for a fundamental mode to be the first to
become unstable. The results, both numerical and from lin-
ear theory, show that as |Z| decreases, the wavelengths of
each mode increase. Thus the position within the resonator
maximizing power shifts away from the closed end for each
mode. Hence, the optimal position of the active cell for
destabilizing any mode (λ<sub>n</sub>/8, n = 1, 2, 3, ... ) shifts to the
right, as was already shown in Fig. 8. Mode n will remain
the most unstable until λ<sub>n</sub>/8 of the n-th mode shifts too far
away from the active cell location and λ<sub>n+1</sub>/8 of the next
mode becomes more favorably located, resulting in step
increases in the frequency of the most unstable mode as
shown in Fig. 11.

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

<figure><img src="image_95.png" /><figcaption>FIG. 11. Frequency of the most unstable modes vs Z. (□) Linear theory, (×) simulations, (a) −1 < Z < −0.6; (—) resonant frequencies of a closed-closed lossless resonator with temperature discontinuity at $\dot{x} = 0$ and $T_h = T_{crit}(Z)$ determined by numerical simulations, (— —) approximate vertical separation between modes (multiples of $\lambda/2$). (b) $0.6 < Z < 1$, (—) resonant frequencies of an open-closed lossless resonator with temperature discontinuity at $\dot{x} = 0$ and $T_h = T_{crit}(Z)$, (— —) approximate vertical separation between modes (odd multiples of $\lambda/4$).</figcaption></figure>

VI. COMPARISON WITH EXPERIMENTAL RESULTS

In this section, numerical results from the current model
are compared with the experimental results of Atchley et al.⁸
The experimental setup consists of a thermoacoustic cell that
is about 6.5 cm long, with a 3.5 cm long stack inserted inside
a 1 m long resonator closed at both ends, at about one tenth
of its length, 0.87 m from the right end.²⁶ The inner radius of
the resonator is 1.91 cm. Here its cross-sectional area is
denoted by A. The only “load” on the engine is dissipation in
the gas contained in the resonator (helium under pressure,
initially at ambient temperature). The experiments⁸ show
that at low mean pressure, the most unstable mode is the first
harmonic, while at higher pressure, the fundamental mode
becomes the most unstable. There is also a region where
both modes are unstable.

For comparison, simulations were performed for several
mean pressure values. One of the model assumptions is that
acoustics is lossless in the resonators so that all dissipation is
due to the load, while the experimental results⁸ were
obtained with no load but with losses. For comparison, the
load value in the present model was identified such that the
critical temperature ratio obtained by simulation matches the experiment. For $\tilde{P}_m = 0.44$ MPa, it is found that a dimensionless load value of $f=45$ in the present model yields a critical heater temperature that coincides with the value obtained experimentally, $\tilde{T}_h = 660$ K. This value corresponds to a dimensional resistive load value $\tilde{R} = fM\tilde{P}_m/(A\tilde{U}_{ref}) = 17.1 \text{ MPa s m}^{-3}$. The mode that becomes unstable is the fundamental mode, which is also consistent with the experiments. As the heater temperature is increased in the simulation while keeping the load value constant, the first harmonic also becomes unstable when $\tilde{T}_h \approx 821$ K, which is again in agreement with the value $\tilde{T}_h \approx 800$ K found in the experiment.

For $\tilde{P}_m = 0.15$ MPa, a load value $f=16$ is necessary to reproduce the measured critical heater temperature $\tilde{T}_h = 718$ K, corresponding to dimensional value $\tilde{R} = 2.1 \text{ MPa s m}^{-3}$. For that mean pressure value, the mode that first becomes unstable is the first harmonic, again in the experiments as well as in the present simulations. The load value that matches the experimental stability diagram of Atchley et al.,⁸ shown in Fig. 12, is found to depend only on the mean pressure. In the electro-acoustic analogy, the cold resonator in the experimental configuration is a simple right- ended duct filled with helium prone to viscous and thermal dissipation. In the present model, the lossless resonator is connected to a purely resistive load, which is found to mainly play the role of thermal resistance. Indeed the ther- mal resistance in the experiment may be estimated based on the classical expression,<fnref n="1" /> $R_κ = 2γ\tilde{P}_m/[ω(γ-1)Sδ_κ]$, where S is the surface area of the resonator side and $δ_κ$ the thermal boundary-layer thickness. If $\tilde{P}_m = 0.44$ MPa (respectively, 0.15 MPa), $\tilde{T}_{cold} = 293$ K, the thermal resistance is found to be on the order of 40 MPa s m⁻<fnref n="3" /> (respectively, 5 MPa s m⁻<fnref n="3" />). For $\tilde{P}_m = 0.15$ MPa, the angular frequency cor- responds to the first harmonic, since this is the mode that destabilizes first. These values are found to be in qualitative agreement with the abovementioned values $\tilde{R}$ obtained numerically, thus lossless acoustics with a loaded end can adequately model a dissipative resonator. Moreover, the value of the load strongly influences the ability of the ther- moacoustic engine to start on its own.

<figure><img src="image_96.png" /><figcaption>FIG. 12. Dimensional load $\tilde{R}$ yielding the critical temperatures of Atchley et al. (Ref. 8) vs mean pressure $\tilde{P}_m$.</figcaption></figure>

## VII. CONCLUSION

The dynamics of a thermoacoustic engine was investi- gated numerically for varying mean pressure and tempera- ture ratio between heater and cooler. The model is based upon a low Mach number assumption, leading to a combina- tion between a numerical solution of the flow and heat transfer in the active cell and an analytical solution of linear acoustics in the resonator. Numerical simulations were carried out on existing thermoacoustic engine configura- tions,²²,²⁶ showing how a resistive load influences the onset of thermoacoustic oscillations.

The resistive load was varied to encompass a full range of situations from a resonator with a closed end to an open end, hence values range from infinity to zero. To that effect, for a large number of load values, the critical temperature ratio was determined, as well as the angular frequency of the most unstable mode. Results show that for a wide range of intermediate load values, the critical temperature ratio is so high that, in practice, the engine will not start on its own.

The active cell consisted of a stack of parallel plates placed between heat exchangers also made of plates aligned with the stack plates. For this simple configuration, compari- son was possible between the numerical results and a theoretical analysis based on Rott’s linear theory. Good agreement was found for both the critical temperature ratio and the dominant oscillation frequency. The hybrid numeri- cal model will be most useful whenever a realistic active cell with complex geometry must be dealt with, for example, for misaligned plates or even fully three-dimensional configurations.

The current model of a lossless resonator associated with a resistance placed at the resonator end can be extended to describe real loads to which the thermoacoustic engine is coupled. The approach used also allows one to account for dissipation in a simple experimental resonator with a closed end,⁸ provided the resistive load value is adequately chosen. In that case, the load is mainly associated with the thermal resistance of the resonator.

## ACKNOWLEDGMENTS

This work was financed by a Ph.D. grant from Université Pierre et Marie Curie. Support of the Natural Sciences and Engineering Research Council of Canada is gratefully acknowledged (L.B.).

<footnote n="1">¹G. W. Swift, “Thermoacoustic engines,” J. Acoust. Soc. Am. **84**, 1145–1180 (1988).</footnote>

<footnote n="2">²R. Radebaugh, “A review of pulse tube refrigeration,” Adv. Cryog. Eng. **35B**, 1191–1205 (1990).</footnote>

<footnote n="3">³A. Miglori and G. W. Swift, “Liquid sodium thermoacoustic engine,” Appl. Phys. Lett. **53**, 355–357 (1988).</footnote>

<footnote n="4">⁴C. Jensen and R. Raspet, “Thermoacoustic power conversion using piezo-electric transducer,” J. Acoust. Soc. Am. **128**(1), 98–103 (2010).</footnote>

<footnote n="5">⁵N. Rott, “Damped and thermally driven acoustic oscillations in wide and narrow tubes,” Z. Angew. Math. Phys. **20**, 230–243 (1969).</footnote>

<footnote n="6">⁶W. P. Arnott, H. E. Bass, and R. Raspet, “General formulation of thermoacoustics for stacks having arbitrarily shaped pore cross section,” J. Acoust. Soc. Am. **90**, 3228–3237 (1991).</footnote>

<footnote n="7">⁷B. Ward, J. Clark, and G. W. Swift, “Design environment for low-amplitude thermoacoustic energy conversion (DeltaEC) user’s guide (version 6.2),” Los Alamos National Laboratory (2008).</footnote>
<footnote n="8">⁸A. A. Atchley and F. M. Kuo, "Stability curves for a thermoacoustic prime mover," *J. Acoust. Soc. Am.* **95**, 1401–1404 (1994).</footnote>

<footnote n="9">⁹W. P. Arnott, J. R. Belcher, R. Raspet, and H. E. Bass, "Stability analysis of a helium-filled thermoacoustic engine," *J. Acoust. Soc. Am.* **96**, 370–375 (1994).</footnote>

<footnote n="10">¹⁰T. Yazaki, A. Iwata, T. Maekawa, and T. Tominaga, "Travelling wave thermoacoustic engine in a looped tube T," *Phys. Rev. Lett.* **81**, 3128–3131 (1998).</footnote>

<footnote n="11">¹¹N. Sugimoto, "Thermoacoustic-wave equations for gas in a channel and a tube subject to temperature gradient," *J. Fluid Mech.* **658**, 89–116 (2010).</footnote>

<footnote n="12">¹²H. Hyodo and N. Sugimoto, "Stability analysis for the onset of thermoacoustic oscillations in a gas-filled looped tube," *J. Fluid Mech.* **741**, 585–618 (2014).</footnote>

<footnote n="13">¹³A. T. A. M. de Waele, "Basic treatment of onset conditions and transient effects in thermoacoustic Stirling engines," *J. Sound Vib.* **325**, 974–988 (2009).</footnote>

<footnote n="14">¹⁴M. Guedra, G. Penelet, P. Lotton, and J.-P. Dalmont, "Theoretical prediction of the onset of thermoacoustic instability from the experimental transfer matrix of a thermoacoustic core," *J. Acoust. Soc. Am.* **130**, 145–152 (2011).</footnote>

<footnote n="15">¹⁵M. Guedra, G. Penelet, and P. Lotton, "Experimental and theoretical study of the dynamics of self-sustained oscillations in a standing wave thermoacoustic engine," *J. Appl. Phys.* **115**, 024504 (2014).</footnote>

<footnote n="16">¹⁶J. R. Olson and G. W. Swift, "A loaded thermoacoustic engine," *J. Acoust. Soc. Am.* **98**, 2690–2693 (1995).</footnote>

<footnote n="17">¹⁷R. Bao, G. B. Chen, K. Tang, Z. Z. Jia, and W. H. Cao, "Effect of RC load on performance of thermoacoustic engine," *Cryogenics* **46**, 666–671 (2006).</footnote>

<footnote n="18">¹⁸P. Duthil, D. Baltean Carlès, A. Bétrancourt, M.-X. François, Z. B. Yu, and J. P. Thermeau, "Experiments and valve modeling in thermoacoustic device," *Adv. Cryog. Eng.* **51**(B), 1095–1102 (2006).</footnote>

<footnote n="19">¹⁹H. Hatori, T. Biwa, and T. Yazaki, "How to build a loaded thermoacoustic engine," *J. Appl. Phys.* **111**, 074905 (2012).</footnote>

<footnote n="20">²⁰G. Y. Yu, E. C. Luo, W. Dai, and J. Y. Hu, "Study of nonlinear processes of a large experimental thermoacoustic-Stirling heat engine by using computational fluid dynamics," *J. Appl. Phys.* **102**, 074901 (2007).</footnote>

<footnote n="21">²¹J. A. Lycklama à Nijeholt, M. E. H. Tijani, and S. Spoelstra, "Simulation of a traveling-wave thermoacoustic engine using computational fluid dynamics," *J. Acoust. Soc. Am.* **118**, 2265–2270 (2005).</footnote>

<footnote n="22">²²O. Hireche, C. Weisman, D. Baltean-Carlès, P. Le Quéré, and L. Bauwens, "Low Mach number analysis of idealized thermoacoustic engines with numerical solution," *J. Acoust. Soc. Am.* **128**(6), 3438–3448 (2010).</footnote>

<footnote n="23">²³R. Raspet, J. Bruster, and H. E. Bass, "New approximation method for thermoacoustic calculations," *J. Acoust. Soc. Am.* **103**(5), 2395–2402 (1998).</footnote>

<footnote n="24">²⁴N. Rott, "Thermoacoustics," *Adv. Appl. Mech.* **20**, 135–243 (1980).</footnote>

<footnote n="25">²⁵G. W. Swift, *A Unifying Perspective for Some Engines and Refrigerators* (Acoustical Society of America, Melville, NY, 2002), p. 90, Eq. 4.70.</footnote>

<footnote n="26">²⁶A. A. Atchley, H. E. Bass, T. J. Hofler, and H. T. Lin, "Study of a thermoacoustic prime mover below onset of self-oscillation," *J. Acoust. Soc. Am.* **91**, 734–743 (1992).</footnote>
V. Daru, D. Baltean-Carlès, C. Weisman, P. Debesse, G. Gandikota

Two-dimensional numerical simulations of nonlinear acoustic streaming in standing waves
Wave Motion 50, 955-963 (2013).
# Two-dimensional numerical simulations of nonlinear acoustic streaming in standing waves

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

CrossMark

Virginie Daru<sup>a,b,*</sup>, Diana Baltean-Carlès<sup>a,c</sup>, Catherine Weisman<sup>a,c</sup>, Philippe Debesse<sup>a</sup>, Gurunath Gandikota<sup>a</sup>

## HIGHLIGHTS

* Numerical simulations of acoustic streaming in standing wave 2D channels are shown.

* Compressible Navier-Stokes equations are solved with high-order differential schemes.

* The transition from regular to irregular streaming flow regimes is described.

* Simulations are in good alignment with experiments of nonlinear streaming regimes.

* Intricate coupling exists between the mean temperature field and the streaming flow.

## ARTICLE INFO

**Article history:**

Received 4 August 2012
Received in revised form 17 March 2013
Accepted 21 March 2013
Available online 1 April 2013

**Keywords:**

Acoustic streaming
Standing wave
Numerical simulation
Nonlinear streaming regime

## ABSTRACT

Numerical simulations of compressible Navier-Stokes equations in closed two-dimensional channels are performed. A plane standing wave is excited inside the channel and the associated acoustic streaming is investigated for high intensity waves, in the nonlinear streaming regime. Significant distortion of streaming cells is observed, with the centers of streaming cells pushed toward the end-walls. The mean temperature evolution associated with the streaming motion is also investigated.

© 2013 Elsevier B.V. All rights reserved.

## 1. Introduction

Acoustic streaming is generated inside a two-dimensional channel as a consequence of the interaction between a plane standing wave and the solid boundaries. It consists of a mean second order flow produced mainly by shear forces within the viscous boundary layer along the solid walls. This motion was initially studied by Rayleigh [1] in the case of wide channels, in which the boundary layer thickness is negligible in comparison with the channel width. This streaming flow is characterized by four steady counter-rotating vortices in one-half wavelength of the standing wave outside the boundary layer, nowadays referred to as Rayleigh streaming. The vortices develop along the half wavelength of the standing wave. Along the central axis of the channel, the streaming motion is oriented from axial particle velocity nodes to antinodes. Inside the boundary layer four additional vortices are created simultaneously, with the streaming motion oriented from axial particle velocity antinodes to nodes along the inner walls of the tube [2,3].

<footnote><sup>a</sup> LIMSI-UPR CNRS 3251, BP133, 91403 Orsay Cedex, France</footnote>

<footnote><sup>b</sup> Arts et Métiers ParisTech, Lab. DynFluid, 151 bd de l'hôpital, 75013 Paris, France</footnote>

<footnote><sup>c</sup> Université Pierre et Marie Curie, 4 Place Jussieu, 75252 Paris Cedex 05, France</footnote>

<footnote>* Corresponding author at: LIMSI-UPR CNRS 3251, BP133, 91403 Orsay Cedex, France. Tel.: +33 169858078.
E-mail address: virginie.daru@limsi.fr (V. Daru).</footnote>
In the case of wide channels, Menguy and Gilbert [4] showed that streaming itself can be linear (case of slow streaming) or nonlinear (case of fast streaming), and both regimes are characterized by a reference nonlinear Reynolds number $Re_{NL} = (M \times y_0/\delta_v)^2$ reflecting the influence of inertial effects on the streaming flow (M is the acoustic Mach number, $M = U_{max}/c_0$, with $U_{max}$ the maximum axial particle velocity inside the channel and $c_0$ the initial speed of sound, $y_0$ is the half width of the channel and $\delta_v$ the viscous boundary layer thickness). The value of $Re_{NL}$ determines the magnitude of the nonlinear hydrodynamic effect of the fluid inertia in the time-averaged momentum equation for the streaming velocity field [4]. Most analytical streaming models have been established in the case of slow streaming, characterized by $Re_{NL} \ll 1$. They are based on successive approximations of the nonlinear hydrodynamic equations and have been derived for arbitrary values of the ratio $y_0/\delta_v$, taking into account the variations of heat conduction and viscosity with temperature [5], and the existence of a longitudinal temperature gradient [6]. In the case $Re_{NL} = O(1)$, Menguy and Gilbert [4] derived an asymptotic model for streaming flow inside wide cylindrical resonators, with no mean temperature gradient, and showed a distortion of streaming patterns due to inertia effects. However, this model does not cover the strongly nonlinear streaming regime ($Re_{NL} \gg 1$), and does not explain the nonlinear effects on acoustic streaming recently observed in several experimental works [7–9], where the temperature gradient along the resonator wall has a significant influence.

Numerical simulations in the linear regime, yielding results for non idealized geometries, were performed in the specific cases of thermoacoustic refrigerators [10] or in annular resonators [11] and solved the dynamics of the flow without taking heat transfer into account.

Simulations in the nonlinear regime were first performed by Yano [12], who studied the acoustic streaming associated with resonant oscillations with periodic shock waves in tubes with aspect ratio (width over length) very large (0.1). He solved the full 2D Navier–Stokes equations with an upwind finite-difference TVD scheme and showed the existence of irregular vortex structures and even turbulent streaming for high streaming Reynolds numbers (based on a characteristic streaming velocity, the tube length, and the kinematic viscosity, $R_s = U_sL/\nu$). This is a different configuration than our configuration, since it considers low frequency acoustic waves in wide tubes with respect to their length and focuses on turbulent streaming.

Simulations of acoustic streaming in the linear and nonlinear regime, taking heat transfer into account, in a two-dimensional rectangular enclosure, were performed by Aktas and Farouk [13]. In their study, the wave is created by vibrating the left wall of the enclosure and the full compressible Navier–Stokes equations are solved, with an explicit time-marching algorithm (a fourth order flux-corrected transport algorithm) to track the acoustic waves. Their numerical results are in agreement with theoretical results in the linear regime and show irregular streaming motion in the nonlinear regime, but they show the existence of irregular streaming at small values of $Re_{NL}$, in contradiction with experiments cited above. Moreover, these simulations do not analyze the deformation of the streaming cells until they split onto several cells.

We propose in this work to conduct numerical 2D compressible simulations for studying the origin of the distortion of streaming cells (of Rayleigh type) that were experimentally observed. Calculations are performed for channels with aspect ratios ranging from 0.01 to 0.07, and the coupling between streaming effects and thermal effects in the channel (existence of a mean temperature gradient) is also investigated.

## 2. Problem description and numerical model

We consider a rectangular channel of length $L$ and half width $y_0$, initially filled with the working gas. In order to initiate an acoustic standing wave in the channel, it is shaken in the longitudinal direction ($x$), so that an harmonic velocity law is imposed, $V(t) = (V(t), 0)^T$, with $V(t) = x_p\omega \cos(\omega t)$, $\omega$ being the angular frequency and $x_p$ the amplitude of the channel displacement. The channel being undeformable, the flow can be modeled by the compressible Navier–Stokes equations expressed in the moving frame attached to the channel, so that a forcing source term is added. The model reads

$$ \left\{ \begin{aligned} \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) &= 0 \\ \frac{\partial (\rho \mathbf{v})}{\partial t} + \nabla \cdot (\rho \mathbf{v} \otimes \mathbf{v}) + \nabla p &= \nabla \cdot (\bar{\tau}) - \rho \frac{d\mathbf{V}}{dt} \\ \frac{\partial \rho E}{\partial t} + \nabla \cdot (\rho E \mathbf{v} + p \mathbf{v}) &= \nabla \cdot (k \nabla T) + \nabla \cdot (\bar{\tau} \mathbf{v}) - \rho \mathbf{v} \cdot \frac{d\mathbf{V}}{dt} \end{aligned} \right. \quad (1) $$

where $\mathbf{v} = (u, v)^T$ is the flow velocity, $E = e + \frac{1}{2}\mathbf{v} \cdot \mathbf{v}$ is the total energy, with $e = \frac{p}{(\gamma-1)\rho}$ the internal energy, $\gamma$ the specific heat ratio, $\bar{\tau} = -\frac{2}{3}\mu(\nabla \cdot \mathbf{v})\bar{I} + 2\mu\bar{D}$ the viscous stress tensor of a Newtonian fluid, $\bar{D}$ the strain tensor, $\mu$ the dynamic viscosity, $k$ the thermal conductivity. The thermo-physical properties $\mu$ and $k$ are supposed to be constant. The gas is considered as a perfect gas obeying the state law $p = r\rho T$, where $T$ is the temperature and $r$ is the perfect gas constant corresponding to the working gas. The physical boundary conditions employed in the moving frame are no slip and isothermal walls.

The model is numerically solved by using high order finite difference schemes, developed in Daru and Tenaud [14]. An upwind scheme, third order accurate in time and space, is used for convective terms, and a centered scheme, second order, is used for diffusion terms. More detail about the scheme and computations showing its good qualities can be found in Daru and Gloerfelt [15] and Daru and Tenaud [16]. This scheme can be derived up to an arbitrary order of accuracy for convective terms in the case of a scalar equation. Here the third order scheme is selected, after having done several comparisons using
higher order schemes (up to the 11th order), that have shown that third order gives sufficient accuracy for a reasonable CPU cost. In cases where shock waves are present, the scheme can be equipped with a flux limiter (MP), preserving monotonicity, intended for suppressing the parasitic numerical oscillations generated in the shock region, while preserving the accuracy of the scheme in smooth regions. However, the flows considered here are always low Mach number flows. Although traveling shock waves are a main feature of the flow for high acoustics levels, as noticed by several authors [17], they are of weak intensity and the numerical oscillations are very small and do not spoil the solution. Thus the MP limiter, which is expensive in terms of CPU cost, was not activated in these calculations. For solving the 2D Navier–Stokes equations, the scheme is implemented using Strang splitting. This reduces the formal accuracy of the scheme to second order. However, numerical experiments have shown that a very low level of error is still achieved.

Let us describe our numerical procedure. The system (1) can be written in vector form:

$$ \frac{\partial w}{\partial t} + \frac{\partial}{\partial x}(f - f^v) + \frac{\partial}{\partial y}(g - g^v) = h \quad (2) $$

where w is the vector of conservative variables ($\rho$, $\rho u$, $\rho v$, $\rho E)^T$, f and g are the inviscid fluxes $f = (\rho u, \rho u^2+p, \rho uv, \rho Eu+pu)^T$ and $g = (\rho v, \rho uv, \rho v^2+p, \rho Ev + pv)^T$, $f^v$ and $g^v$ being the viscous fluxes $f^v = (0, \tau_{xx}, \tau_{xy}, k \frac{\partial T}{\partial x} + u \tau_{xx} + v \tau_{xy})^T$, $g^v = (0, \tau_{xy}, \tau_{yy}, k \frac{\partial T}{\partial y} + u \tau_{xy} + v \tau_{yy})^T$. The source term reads $h = (0, -\rho \frac{dV}{dt}, 0, -\rho u \frac{dV}{dt})^T$. Denoting by $w_{ij}^n$ the numerical solution at time $t = n\delta t$ and grid point $(x, y) = (i\delta x, j\delta y)$, we use the following Strang splitting procedure to obtain second order of accuracy every two time steps:

$$ w_{i,j}^{n+2} = L_{\delta x} L_{\delta y} L_{\delta x} w_{ij}^n \quad (3) $$

where $L_{\delta x}$ (resp. $L_{\delta y}$) is a discrete approximation of $L_x(w) = w + \delta t(-f_x + f_x^v + h)$ (resp. $L_y(w) = w + \delta t(-g_y + g_y^v)$). The 1D operators being similar in the two directions, we only describe the x operator. The scheme is implemented as a correction to the second order MacCormack scheme. It consists of three steps, as follows:

$$ \begin{aligned} w_{i,j}^* &= w_{i,j}^n - \frac{\delta t}{\delta x}(f_{i+1,j} - f_{i,j} - f_{i+1/2,j}^v + f_{i-1/2,j}^v)^n + \delta t h_{i,j}^n \\ w_{i,j}^{**} &= w_{i,j}^* - \frac{\delta t}{\delta x}(f_{i,j} - f_{i-1,j} - f_{i+1/2,j}^v + f_{i-1/2,j}^v)^* \end{aligned} \quad (4) $$

$$ w_{i,j}^{n+1} = \frac{1}{2}(w_{i,j}^n + w_{i,j}^{**}) + C_{i+1/2,j}^x - C_{i-1/2,j}^x. $$

The viscous fluxes are discretized at each interface using centered second order finite differences formulas. The corrective term $C_{i+1/2,j}^x - C_{i-1/2,j}^x$ provides the third order accuracy and the upwinding for the inviscid terms. Let us define $\psi_{i+1/2,j} = \frac{1}{6} \sum_{l=1}^4 \{ |v_{i+1/2,j}^l| (1 - v_{i+1/2,j}^l) (1 + v_{i+1/2,j}^l) \delta \alpha_{i+1/2,j}^l \cdot d_{i+1/2,j}^l \}$, where $v^l = \frac{\delta t}{\delta x} \lambda^l$, $\lambda^l$ and $d^l$ are the eigenvalues and eigenvectors of the Roe-averaged Jacobian matrix $A = \frac{df}{d\bar{w}}$ [18], and $\delta \alpha^l$ is the contribution of the l-wave to the variation $(w_{i+1/2,j}^n - w_{i-1/2,j}^n)$. Using the function $\psi$, the corrective term reads

$$ C_{i+1/2,j}^x = \begin{cases} -\psi_{i+1/2,j}^n + \psi_{i-1/2,j}^n & \text{if } v_{i+1/2,j} \ge 0 \\ \psi_{i+3/2,j}^n - \psi_{i+1/2,j}^n & \text{if } v_{i+1/2,j} < 0. \end{cases} \quad (5) $$

This completes the description of the numerical method.

We are interested in the acoustic streaming generated by the interaction of the imposed plane standing wave and the channel wall. Resonant conditions are imposed, for which $L = \lambda/2$, $\lambda = c_0/f$ being the wavelength, $c_0$ the speed of sound for initial state and f the vibration frequency of the channel. It is known [5] that boundary layers develop along the walls, with thickness $\delta_v = \sqrt{2\nu/\omega}$, $\nu$ being the kinematic viscosity $\nu = \mu/\rho_0$, and $\rho_0$ the density at initial state. Depending on the value of the ratio $y_0/\delta_v$, several patterns of streaming can appear: Rayleigh-type streaming in the central region, and boundary layer type streaming near the longitudinal walls. The boundary layer is of small thickness and must be correctly resolved by the discretization mesh. After several trials, we have determined that a value of 5 points per boundary layer thickness is sufficient for reasonable accuracy of the simulations. The results obtained using 10 points per boundary layer thickness show very small differences with the former, the maximum value of the differences being less than 3%. All results presented below are thus obtained using a cartesian mesh of rectangular cells of constant size $\delta x$ and $\delta y$, composed of 500 points in the axial direction x, and of $5 \times y_0/\delta_v$ points in the y direction normal to the axis. In the considered geometry, this leads to cells such that $\delta x = 5.8 \delta y$, which is sufficient to correctly solve the equations in the region close to the vertical walls. The flow being symmetrical with respect to the x axis (at least in the range of parameters treated), only the upper half of the channel was considered. Also, the scheme being fully explicit, the time step $\delta t$ is fixed such as to satisfy the stability condition of the scheme which can be written as

$$ \delta t \le \frac{1}{2} \min(\delta y^2/\nu, \delta y^2/(k/\rho_0 c_0), \delta y/c_0). \quad (6) $$
<table><caption>**Table 1**
Values of the parameters of the simulations.</caption><tr><td>x<sub>p</sub> (μm)</td><td>5</td><td>10</td><td>50</td><td>80</td><td>100</td></tr><tr><td>U<sub>max</sub> (m/s)</td><td>7.0</td><td>26.45</td><td>61.11</td><td>70.94</td><td>89.96</td></tr><tr><td>y<sub>0</sub>/δ<sub>v</sub></td><td>10</td><td>40</td><td>20</td><td>10</td><td>20</td></tr><tr><td>y<sub>0</sub>/L</td><td>0.0180</td><td>0.0718</td><td>0.0359</td><td>0.018</td><td>0.0359</td></tr><tr><td>Re<sub>NL</sub></td><td>0.041</td><td>9.469</td><td>12.636</td><td>4.257</td><td>27.384</td></tr></table>

As shown in Eq. (6), the first two limiting values $\delta y^2/\nu$ and $\delta y^2/(k/\rho_0 c_0)$ are related to the viscous and thermal conduction terms, and the third one $\delta y/c_0$ is related to the acoustic propagation. In all cases considered here, the time step limitation is acoustic, i.e. $\delta t \le \frac{1}{2}\delta y/c_0$. Taking $\delta t = \frac{1}{2}\delta y/c_0$ and $\delta y = \delta_\nu/5$, this results in a number of time steps $N_T$ per period of oscillation proportional to $\sqrt{L}$, $N_T = 1/(f\delta t) = 10\sqrt{\frac{2\pi c_0}{\nu}}\sqrt{L}$. Since transients of several hundreds of periods may be needed in order to reach stabilized steady streaming flow, simulations are very costly, and one must rely on numerical schemes that are sufficiently accurate in both space and time.

Finally, the mean flow is obtained from calculating a simple mean value for each physical quantity (velocity, pressure, temperature) over an acoustic period, which for example in the case of velocity can be written as

$$ \bar{\mathbf{v}}_{i,j}^E = \frac{1}{N_T} \sum_{n=1}^{N_T} \mathbf{v}_{i,j}^n. \qquad (7) $$

The mean velocity obtained in Eq. (7) can be understood as a “Eulerian” time-average velocity. Similarly, an “average mass transport” velocity, associated with the streaming mass transport, could also be computed, as the Eulerian time-average of mass flux vector divided by the Eulerian time-average of gas density:

$$ \bar{\mathbf{v}}_{i,j}^M = \frac{\sum_{n=1}^{N_T} \rho_{i,j}^n \mathbf{v}_{i,j}^n}{\sum_{n=1}^{N_T} \rho_{i,j}^n}. \qquad (8) $$

For weakly nonlinear regimes, the difference between these two velocities is significant only in the boundary layer, and in the case of wide channels they are almost the same. In order to observe the mechanism of cell distortion, either one of these velocities can be monitored.

### 3. Numerical results

We consider a channel initially filled with air at standard thermodynamic conditions, $p_0 = 101,325$ Pa, $\rho_0 = 1.2$ kg m$^{-3}$, $T_0 = 294.15$ K. The thermo-physical properties of air are $\mu = 1.795 \times 10^{-5}$ kg m$^{-1}$ s$^{-1}$ and $k = 0.025$ W m$^{-1}$ K$^{-1}$. Also for air, $\gamma = 1.4$ and $r = 287.06$ J kg$^{-1}$ K$^{-1}$. The Prandtl number $Pr$ is equal to 0.726. This results in an initial speed of sound $c_0 = 343.82$ m s$^{-1}$.

For a 1 m long channel, the limiting time step would correspond to $N_T \approx 25,000\sqrt{L}$, that is 25,000 iterations per period. Since transients of several hundreds of periods may be needed in order to reach stabilized streaming flow, several millions of iterations are necessary for each simulation. Considering these numerical constraints, a shorter channel is considered, with $L = 8.59$ mm. This corresponds to a high-frequency wave, with $f = 20,000$ Hz. The resulting boundary layer thickness is $\delta_\nu = 1.54 \times 10^{-5}$ m. The time step $\delta t = 8 \times 10^{-9}$ s is chosen in order to satisfy the numerical stability condition, corresponding to 6250 time iterations per period for a mesh involving 5 grid points across the boundary layer thickness. The axial particle velocity produced in the channel depends on the amplitude of the channel displacement and on the ratio $y_0/\delta_\nu$. It varies approximately linearly with the amplitude of the channel displacement, for a given ratio $y_0/\delta_\nu$. Table 1 summarizes the different parameter values corresponding to the simulations that are presented thereafter.

As mentioned earlier, the parameter identified as relevant in describing the regularity of streaming flow is the nonlinear Reynolds number $Re_{NL}$ introduced by Menguy and Gilbert [4]. In this paper we used a slightly different definition for $Re_{NL}$, because the definition of the viscous boundary layer thickness is different. Our Reynolds number corresponds to half of that of Menguy and Gilbert [4].

We first present results concerning the main acoustic field in the channel, for a small value of $Re_{NL}$ corresponding to slow streaming. In Fig. 1(a) is represented the velocity signal at the center of the channel, as a function of the number of periods elapsed. At this location, the axial particle velocity amplitude is maximum since it corresponds to the antinode. For this value of $Re_{NL}$, the problem is nearly linear and the final signal is purely sinusoidal, in agreement with the linear theory. The amplification of the initial perturbation until saturation can be observed. The periodic regime is established after about 20 periods. Fig. 1(b) shows the time evolution of the mean horizontal velocity (over an acoustic period) and of the mean temperature difference $\Delta T = T - T_0$ (also over an acoustic period) on the axis, at $x = \lambda/8$. At this location the streaming velocity is maximum. It can be noticed that the steady streaming field is established also after about 20 periods which is
<figure><img src="image_98.png" /><figcaption>**Fig. 1.** (a) Axial particle velocity at the channel's center, as a function of time counted by the number of periods elapsed. (b) Mean horizontal velocity and mean temperature variation, on the channel's axis at $x = \lambda/8$. Case $Re_{NL} = 0.041$ ($y_0/\delta_v = 10$, $M = 0.02$).</figcaption></figure>

<figure><img src="image_99.png" /><figcaption>**Fig. 2.** (a) $U_{st}$ as a function of $y/y_0$ at $x = \lambda/8$, numerical (present study) and analytical [20] results. (b) Dimensionless mean fluctuating pressure, $p - p_0$ along the channel's axis. Case $Re_{NL} = 0.041$ ($y_0/\delta_v = 10$, $M = 0.02$).</figcaption></figure>

of the same order of magnitude as the theoretical characteristic streaming timescale $\tau_c = (\frac{2y_0}{\pi})^{\frac{2}{\nu}}$ (see Amari, Gusev and Joly [19]) which in this case gives $n_c$ periods for reaching steady-state, with $n_c = 13$. In Fig. 2(a) is shown the variation of the axial dimensionless streaming velocity at $x = \lambda/8$ along the channel's width, compared with results computed using the analytical expressions of Hamilton, Ilinskii and Zabolotskaya [20]. In this figure, the reference velocity is the Rayleigh streaming reference velocity [2,5], $u_{Rayleigh} = \frac{3}{16}U_{max}^2/c_0$. The slight discrepancy between the numerical and the analytical profiles is probably due to the presence of the vertical end walls, which is not accounted for in the model of Hamilton, Ilinskii and Zabolotskaya [20]. In Fig. 2(b) is shown the stabilized mean pressure $p - p_0$ (over an acoustic period), scaled by $(\gamma/4)p_0 M^2$, along the channel's axis. It is the second order average pressure resulting from the streaming flow, which is clearly one-dimensional and has a cosine variation with respect to $x$, as expected in the linear regime of streaming. In the present case, there is an offset pressure $p_{off}$, corresponding to an increase of the mean pressure and temperature inside the channel, due to the harmonic forcing source term. When subtracting off this offset pressure, the theoretical result for the dimensionless hydrodynamic streaming pressure is obtained, $P_{2s} = \cos 4\pi \frac{x}{\lambda}$ (see Menguy and Gilbert [4]).

Simulations are then performed for several values of $Re_{NL}$ corresponding to configurations ranging from slow streaming flow ($Re_{NL} = 0.041$) to fast streaming flow ($Re_{NL} = 27.384$), for several values of the cavity width ($y_0/\delta_v = 10, 20, 40$), and for increasing acoustic velocities, with Mach numbers ranging from $M = 0.02$ to $M = 0.27$ so that shock waves can occur. This can be seen in Fig. 3(right) showing the axial velocity signal at channel's center as a function of time counted by the number of periods elapsed. In Fig. 3(a)(right) the signal contains only one frequency, but for all other cases, there are shock waves and the axial velocity signal is distorted in a “U” shape, because of the presence of odd harmonics (3, 5, etc.). Fig. 3(left) shows the streamlines of the streaming velocity field over the whole length and only over the top half width of the channel. As expected, in the case of small $Re_{NL}$ number values (Fig. 3(a)), four symmetric streaming cells develop over the length and the half width of the channel: two cells in the boundary layer, and two cells in the core of the channel, identified in the literature as Rayleigh streaming. These results are in agreement with the predictions of analytical models of streaming flows [5,6], and with experimental measurements [7]. The only noticeable difference is the slight asymmetry of cells with respect to the vertical lines $x = \lambda/8$ and $x = 3\lambda/8$, due to the presence of boundary layers along the vertical end walls. Indeed the effects of the vertical end walls are accounted for in the present simulations but are neglected in the analytical models and are very far from the measurement area in the experiments. Several simulations have shown that this asymmetry is independent of $Re_{NL}$ as long as the value of the latter remains small with respect to 1. One can note that 2-D steady (time-periodic) streamlines shown in Fig. 3(left) correspond to the Eulerian time-average velocity defined in Eq. (7). They are not closed curves, even though steady streaming motion is completely established. If the streamlines were calculated for the average mass transport velocity of Eq. (8), they would be closed curves, since they would show streamlines of mass transport.

For $Re_{NL} > 1$, the steady streaming flow is established after the same characteristic time as in the linear case. The re-circulation cells become very asymmetric as $Re_{NL}$ increases, and streaming flow becomes irregular (Fig. 3(b)–(e)(left)). This was also observed experimentally (with PIV measurements) by Nabavi, Siddiqui and Dargahi [8] in a rectangular enclosure. The centers of all streaming cells (boundary layer cells as well as central cells) are displaced toward the ends of the resonant
<figure><img src="image_100.png" /><figcaption>**Fig. 3.** Streamlines of mean flow on the top half of the channel (left) and axial particle velocity signal at channel's center as a function of time counted by the number of periods elapsed (right) (a) $Re_{NL} = 0.041$ ($y_0/\delta_v = 10, M = 0.02$). (b) $Re_{NL} = 4.257$ ($y_0/\delta_v = 10, M = 0.206$). (c) $Re_{NL} = 9.469$ ($y_0/\delta_v = 40, M = 0.077$). (d) $Re_{NL} = 12.636$ ($y_0/\delta_v = 20, M = 0.178$). (e) $Re_{NL} = 27.384$ ($y_0/\delta_v = 20, M = 0.262$). Lengths are normalized with L.</figcaption></figure>

channel, close to the boundary layers next to the vertical walls. PIV measurements by Nabavi, Siddiqui and Dargahi [8] show the same distortion of streamlines between an axial particle velocity node and an antinode. Fig. 4(a) shows the $x$ variation along the channel's central axis $y = 0$, of the axial dimensionless streaming velocity component, using as reference velocity the Rayleigh streaming reference velocity [2,5], $u_{Rayleigh} = \frac{3}{16}U_{max}^2/c_0$. There is a clear modification of the velocity profiles as $Re_{NL}$ increases: the sine function associated with slow streaming becomes steeper in the vicinity of both ends. The slope to the curve at the channel's center (axial particle velocity node) becomes smaller as $Re_{NL}$ increases, then becomes close to zero (curve parallel to the longitudinal axis) for a critical value between 13 and 27, and then changes sign, which indicates the emergence of new streaming cells (Fig. 4(a)). Another consequence of the distortion of streaming cells can be observed on the acoustic streaming axial velocity profiles along the width of the channel, shown in Fig. 4(b)–(d). The parabolic behavior in the center of the channel at $x = \lambda/8$ disappears as $Re_{NL}$ increases (see Fig. 4(b)), as a consequence of displacement of the center of each streaming cell toward the velocity node. Fig. 4(c), (d) also confirm the direction of the displacement of the streaming cells' centers. This distortion of streaming cells was already observed in experiments in rectangular or cylindrical geometries in wide channels [7–9]. Nabavi, Siddiqui and Dargahi [8] described it as irregular streaming and detected a critical nonlinear Reynolds number $Re_{NL} = 25$ that separates regular and irregular streaming, which is in agreement with our simulations. In the literature there is to our knowledge no other theoretical or numerical study confirming measurements in these streaming regimes. With the weakly nonlinear model of Menguy and Gilbert [4] the streaming can be calculated for a maximum value of $Re_{NL} = 2$ (in our definition), while the numerical simulations of Aktas and Farouk [13] show the existence of multiple streaming cells for a low value of $Re_{NL} = 1.4$, which is in contradiction with our results and with experiments. Moreover, these numerical simulations [13] do not analyze in detail the transition from two exterior streaming cells to more streaming cells.

According to Menguy and Gilbert [4], the fluid inertia causes distortion of streaming cells for large values of $Re_{NL}$. This was also verified through our simulations. For $Re_{NL} = O(1)$, the Mach number is still small (the wave is almost a monofrequency wave) and the mean temperature difference inside the channel is smaller than 0.1 K (the mean temperature gradient is negligible). The approximations of the model by Menguy and Gilbert [4] still apply here, so we can say that the distortion is caused only by inertial effects. In their experimental study, Thompson, Atchley and Maccarone [9] show the existence of some distortion of the streaming field that are not predicted by existing models of the literature in the nonlinear regime. They do not relate this distortion to fluid inertia but rather to the influence of the mean temperature
<figure><img src="image_101.png" /><figcaption>**Fig. 4.** Horizontal mean velocity component $U_{st}$, normalized with $\frac{3}{16}U_{\text{max}}^2/c_0$ for the 5 cases of Fig. 3. (a) $U_{st}$ along the channel's central axis. (b)–(d) $U_{st}$ as a function of $y/y_0$ for several sections x.</figcaption></figure>

<figure><img src="image_102.png" /><figcaption>**Fig. 5.** Mean temperature field on the top half of the channel (a) $Re_{NL} = 0.041$. (b) $Re_{NL} = 4.257$. (c) $Re_{NL} = 12.636$. (d) $Re_{NL} = 27.383$. The difference between minimum and maximum values of temperature is (respectively): (a) $\Delta T = 0.039$ K, (b) $\Delta T = 6.39$ K, (c) $\Delta T = 19.4$ K, (d) $\Delta T = 44$ K.</figcaption></figure>

field, and more specifically of the axial temperature gradient induced through a thermoacoustic effect along the horizontal walls of the resonating channel. In an experimental case with no shock waves Merkli and Thomann [17] showed that a mean temperature gradient is established inside the tube so that heat is removed close to the velocity antinodes, i.e. at the location of largest viscous dissipation, and heat is produced close to velocity nodes, along the lateral walls. Similar observations can
<figure><img src="image_103.png" /><figcaption>**Fig. 6.** Vector field for the second order hydrodynamic entropy flux, with only one out of twenty vectors shown in the axial direction, and one out of two vectors shown in the transverse direction, for $Re_{NL} = 0.041$.</figcaption></figure>

be made in our simulations. Fig. 5 shows the mean temperature field for increasing values of $Re_{NL}$. Also, for a small value of $Re_{NL}$, Fig. 6 shows the corresponding vector field for the second order hydrodynamic entropy flux (corresponding to the induced thermoacoustic heat transfer in the case of a standing wave [21]). It is defined as the discrete mean value over the acoustic period of $\rho_0 T_0 \Delta s v$, with $\Delta s$ the excess entropy defined as $\Delta s = s - s_0 = c_v \log[(p/p_0)(\rho_0/\rho)^\gamma]$. For small values of $Re_{NL}$, as can be seen in Figs. 5(a) and 6 the thermoacoustic heat transfer takes place at a distance of one thermal boundary layer thickness and then heat diffuses in the radial direction, yielding a temperature field almost one-dimensional in the central part of the tube in the steady-state. Note that the prescribed temperature condition at the vertical walls influences the mean temperature field only close to the walls. The associated thermal boundary layer width remains about the same in all simulation cases, suggesting that it is linked to the acoustic boundary layer width $\delta_v$.

When $Re_{NL}$ increases, periodic shocks appear and the mean temperature gradient becomes important, potentially modifying the streaming patterns. The mean temperature field clearly becomes two-dimensional, as a consequence of both convective heat transport by streaming flow and heat conduction in both directions (Fig. 5(b)). Within the considered range of values of the nonlinear Reynolds number, there is a change of regime for the temperature field before $Re_{NL} = 12$, corresponding to the confinement of outer streaming cell toward the axial particle velocity node. Consequently a zone of very small streaming velocities is generated in the middle of the cavity and that induces the accumulation of heat (Fig. 5(c)). The mean temperature gradient changes the orientation and can cause the splitting of the outer cell into several cells when further increasing $Re_{NL}$ (Fig. 5(d)).

Note that the streaming flow stabilizes in several stages in regimes with high values of the nonlinear Reynolds number. In a first and rapid stage (a few tens of periods), regular streaming flow appears. Then this regular streaming is destabilized along with increasing heterogeneity of the mean temperature field. The steady mean flow stabilizes much later, with timescales related to convection and heat conduction.

## 4. Conclusions

The numerical simulations performed demonstrate the transition from regular acoustic streaming flow toward irregular streaming, in agreement with existing experimental data. These are the first simulations, to our knowledge, in good alignment with experiments of nonlinear streaming regimes. Results show a sizable influence of vertical boundary layers for the chosen configuration. There is also intricate coupling between the mean temperature field and the streaming flow. This coupling effect will be the object of future work. Also, extension of current results in configurations with larger channels is currently in progress.

## Acknowledgments

The authors wish to acknowledge many fruitful discussions with Joël Gilbert and Hélène Bailliet.

## References

[1] Lord Rayleigh, On the circulation of air observed in Kundts tubes, and on some allied acoustical problems, Philos. Trans. R. Soc. Lond. 175 (1884) 1–21.

[2] W.L. Nyborg, Acoustic streaming, in: W.P. Mason (Ed.), Physical Acoustics, Vol. 2B, Academic Press, New York, 1965, pp. 265–331.

[3] H. Schlichting, Berechnung ebener periodischer Grenzschicht–Strommungen [calculation of plane periodic boundary layer streaming], Phys. Zcit. 33 (1932) 327–335.

[4] L. Menguy, J. Gilbert, Non-linear acoustic streaming accompanying a plane stationary wave in a guide, Acta Acust. 86 (2000) 249–259.

[5] M.F. Hamilton, Y.A. Ilinskii, E.A. Zabolotskaya, Acoustic streaming generated by standing waves in two-dimensional channels of arbitrary width, J. Acoust. Soc. Am. 113 (1) (2003) 153–160.

[6] H. Bailliet, V. Gusev, R. Raspet, R.A. Hiller, Acoustic streaming in closed thermoacoustic devices, J. Acoust. Soc. Am. 110 (2001) 1808–1821.

[7] S. Moreau, H. Bailliet, J.-C. Valière, Measurements of inner and outer streaming vortices in a standing waveguide using laser doppler velocimetry, J. Acoust. Soc. Am. 123 (2) (2008) 640–647.

[8] M. Nabavi, K. Siddiqui, J. Dargahi, Analysis of regular and irregular acoustic streaming patterns in a rectangular enclosure, Wave Motion 46 (2009) 312–322.

[9] M.W. Thompson, A.A. Atchley, M.J. Maccarone, Influences of a temperature gradient and fluid inertia on acoustic streaming in a standing wave, J. Acoust. Soc. Am. 117 (4) (2004) 1839–1849.

[10] D. Marx, P. Blanc-Benon, Computation of the mean velocity field above a stack plate in a thermoacoustic refrigerator, C. R. Mecanique 332 (2004) 867–874.
[11] A. Boufermel, N. Joly, P. Lotton, M. Amari, V. Gusev, Velocity of mass transport to model acoustic streaming: numerical application to annular resonators, Acta Acust. United Acust. 97 (2) (2011) 219–227.

[12] T. Yano, Turbulent acoustic streaming excited by resonant gas oscillation with periodic shock waves in a closed tube, J. Acoust. Soc. Am. 106 (1999) L7–L12.

[13] M.K. Aktas, B. Farouk, Numerical simulation of acoustic streaming generated by finite-amplitude resonant oscillations in an enclosure, J. Acoust. Soc. Am. 116 (5) (2004) 2822–2831.

[14] V. Daru, C. Tenaud, High order one-step monotonicity-preserving schemes for unsteady compressible flow calculations, J. Comput. Phys. 193 (2004) 563–594.

[15] V. Daru, X. Gloerfelt, Aeroacoustic computations using a high order shock-capturing scheme, AIAA J. 45 (10) (2007) 2474–2486.

[16] V. Daru, C. Tenaud, Numerical simulation of the viscous shock tube problem by using a high resolution monotonicity preserving scheme, Comput. Fluids 38 (3) (2009) 664–676.

[17] P. Merkli, H. Thomann, Thermoacoustic effects in a resonance tube, J. Fluid Mech. 70 (1) (1975) 161–177.

[18] P.L. Roe, Approximate Riemann solvers, parameter vectors and difference schemes, J. Comput. Phys. 43 (1981) 357–372.

[19] M. Amari, V. Gusev, N. Joly, Temporal dynamics of the sound wind in acoustitron, Acta Acust. United Acust. 89 (2003) 1008–1024.

[20] M.F. Hamilton, Y.A. Ilinskii, E.A. Zabolotskaya, Thermal effects on acoustic streaming in standing waves, J. Acoust. Soc. Am. 114 (6) (2003) 3092–3101.

[21] G.W. Swift, Thermoacoustics: A Unifying Perspective for Some Engines and Refrigerators, Acoustical Society of America, Melville, NY, 2002.
V. Daru, C. Weisman, D. Baltean-Carlès, I. Reyt, H. Bailliet

Inertial effects on acoustic Rayleigh streaming flow : Transient and established regimes
Wave Motion 74, 1-17, 2017
# Inertial effects on acoustic Rayleigh streaming flow: Transient and established regimes

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

CrossMark

Virginie Daru<sup>a,b</sup>, Catherine Weisman<sup>b,c</sup>, Diana Baltean-Carlès<sup>b,c,*</sup>, Ida Reyt<sup>a,d</sup>, Hélène Bailliet<sup>d</sup>

## HIGHLIGHTS

* Nonlinear inertial effects produce patterns different from DNS and experiments.

* Time evolutions of streaming cells near the wall and in the core are analyzed.

* The full streaming equations are solved numerically to obtain the established flow.

* Nonlinear inertial effects produce patterns different from DNS and experiments.

* Nonlinear inertial effects cannot explain the mutation of streaming at high levels.

## ARTICLE INFO

**Article history:**
Received 27 October 2016
Received in revised form 21 April 2017
Accepted 5 June 2017
Available online 19 June 2017

**Keywords:**
Rayleigh streaming
Standing wave
Inertial effects
Nonlinear streaming
Transient evolution

## ABSTRACT

The effect of inertia on Rayleigh streaming generated inside a cylindrical resonator where a mono-frequency standing wave is imposed, is investigated numerically and experimentally. To this effect, time evolutions of streaming cells in the near wall region and in the resonator core are analyzed. An analogy with the lid-driven cavity in a cylindrical geometry is presented in order to analyze the physical meanings of the characteristic times.

Inertial effects on the established streaming flow pattern are then investigated numerically using a code solving the time averaged Navier-Stokes compressible equations, where a mono-frequency acoustic flow field is used to compute the source terms.

It is shown that inertia of streaming cannot be considered as the leading phenomenon to explain the mutation of streaming at high acoustic levels.

© 2017 Elsevier B.V. All rights reserved.

## 1. Introduction

Acoustic streaming is often pointed out as being responsible for misunderstood thermal losses in thermoacoustic machines [1]. Among the different streaming flows found in such devices, Rayleigh streaming consists of a mean flow, of second order, which is produced by the interaction between the acoustic wave and a solid wall [2]. The pioneering work of Rayleigh was continued with the studies, among others, of Schlichting [3], Nyborg [4], Westervelt [5], Rudenko and Soluyan [6], yielding the steady streaming theory. Later, thermal effects were taken into account by Rott [7], the compressibility of the fluid by Qi [8] and more recently semi-analytical models adapted for thermoacoustic applications

<footnote><sup>a</sup> Arts et Métiers ParisTech, Lab. DynFluid, 151 Bd. de l'Hôpital, 75013 Paris, France</footnote>

<footnote><sup>b</sup> LIMSI, CNRS, Université Paris-Saclay, Bât. 508, Rue John Von Neumann, Campus Universitaire, F-91405 Orsay Cedex, France</footnote>

<footnote><sup>c</sup> Sorbonne Universités, UPMC Univ Paris 06, UFR d'Ingénierie, 4 Place Jussieu, 75252 Paris Cedex 05, France</footnote>

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

<footnote>* Corresponding author.
E-mail address: baltean@limsi.fr (D. Baltean-Carlès).</footnote>
<figure><img src="image_105.png" /><figcaption>**Fig. 1.** Streamlines of streaming velocity fields without (a) and with (b) the extra contra-rotating cell (results of numerical simulations $N_1$ and $N_2$).</figcaption></figure>

were developed (see for example [9,10]). These theoretical models, used to describe the streaming flow in the slow regime, are based on the method of successive approximations in powers of the Mach number, performed up to the second order and followed by time averaging over the acoustic period, yielding linear equations for streaming flow.

For high intensity acoustic waves, the streaming flow is referred to as nonlinear streaming or fast streaming. In this case, the effect of fluid inertia was first considered by [11–14] in unbounded or semi-infinite spaces, and more recently by Menguy and Gilbert [15] in resonators and thus for Rayleigh streaming. In the weakly nonlinear study proposed by [15], the method of successive approximations is performed up to the third order and separate solutions are obtained for the near-wall region (inner streaming) and in the resonator core (outer streaming). The global solution is obtained by applying matching conditions to the streaming velocity components. The inner streaming approximate equations were shown to remain linear in the weakly nonlinear case, while the outer streaming flow is affected by inertia [15].

Numerical simulations for acoustic Rayleigh streaming can be found in the literature starting in the 1990's. Two approaches are used: one is based on solving the full Navier-Stokes equations for low and high intensity waves and extract the mean flow from the instantaneous flow [16–18], and recently [19–22] and the other approach is based on solving streaming equations for the mean flow in the linear case [23,24] or in the nonlinear case [25]. Numerical simulations as well as experimental studies (for example [26,27]) in large resonators have shown the two regimes of streaming flows associated with standing waves: slow streaming similar to Rayleigh type streaming or fast streaming where the mean flow is strongly modified when the intensity of the acoustic wave is increased. The dimensionless number commonly used to differentiate between slow and fast streaming is the nonlinear Reynolds number $Re_{NL}$ [15], defined by $Re_{NL} = (M \times R/\delta_v)^2$, where M is the acoustic Mach number, $M = U_{max}/c_0$, with $U_{max}$ the acoustic velocity amplitude at the velocity antinode and $c_0$ the speed of sound, R the guide radius and $\delta_v = \sqrt{2\nu/\omega}$ the acoustic boundary layer thickness, $\omega$ being the pulsation and $\nu$ the kinematic viscosity of the gas. Fig. 1(a) shows the streaming flow pattern (i.e. the streamlines for the steady mean velocity field, with the average performed over an acoustic period) in a half-plane containing the axis of a cylindrical waveguide of axial extension $\lambda/2$ in the slow streaming regime (result from the numerical simulation $N_1$ further described in this paper). This pattern is symmetrical with respect to the guide center and is composed of four streaming cells: two inner toroidal cells close to the tube wall and two outer toroidal cells in the core. The established streaming flow behavior for high acoustic amplitudes (i.e. in the nonlinear regime of streaming) was previously analyzed numerically [19] (with a DNS code) and experimentally [28]: additional streaming cells were observed inside the resonator [29]. Fig. 1(b) shows the distorted streaming pattern (result from the numerical simulation $N_2$ further described in this paper) in the fast streaming regime, when two new streaming cells are observed in the center of the guide, near the acoustic velocity antinode ($z/L = 0$) [19,28,29]. Several sources could potentially be responsible for this mutation of streaming pattern: thermal effects were considered [26], as well as nonlinear propagation [29] and it was concluded that inertial effects on acoustic streaming could be the main source [29]. In the literature there is no study available that allows to conclude about the role of inertia in the behavior of streaming when the latter is “fast”. The current study aims at investigating the hypothesis that inertia is responsible for the mutation of streaming flow at high levels by using both experimental and numerical results.

In Section 2, numerical methods and experimental setup are presented. Section 3 describes the establishment of slow and fast streaming using numerical and experimental observations. Time evolutions of inner and outer streaming cells are discussed. In order to analyze the associated physics, Section 4 shows the streaming flow equations and a scale analysis used to point out the leading parameters. The characteristic time scales are estimated according to the configuration parameters and compared with numerical and experimental results.
<figure><img src="image_106.png" /><figcaption>Fig. 2. Sketch of experimental setup and wave pattern for the working frequency.</figcaption></figure>

In order to further isolate the driving effect of inner streaming on outer streaming, the problem of a flow inside a cylindrical tube with prescribed axial velocity on the walls is solved and analyzed in Section 5. Then the role of inertia on the established nonlinear streaming flow patterns is examined (Section 6) by solving the time-averaged isentropic Navier–Stokes compressible equations (AMS), where a mono-frequency acoustic flow field is used to compute the different terms that are sources of streaming flow. The effect of increasing the amplitude of the acoustic flow field on the streaming flow is analyzed. The results are compared with the predictions of the asymptotic model developed in [15] for the weakly nonlinear regime.

## 2. Methods and set-up

In this section, the experimental setup and the numerical configurations and methods used to investigate the effect of inertia during the establishment of the streaming flow and on the steady state flow are presented briefly (see [29] for more detail). The experimental and numerical configurations are very different in physical dimensions and in frequency but it was checked in [29] that they are asymptotically similar for acoustic and streaming flow. It was shown indeed that the same flow behavior is to be expected if $M \ll 1$, $Re_{ac} \gg 1$ and $Sh \ll 1$ (where $Re_{ac} = \frac{1}{2}(\frac{\lambda}{\pi\delta_v})^2$, $Sh = \frac{\delta_v}{R}$, $\lambda$ being the acoustic wavelength). Both experimental and numerical configurations match these conditions (see Table 1 in [29]).

### 2.1. Experimental set-up

The set-up used to study acoustic streaming consists of a cylindrical glass tube of radius $R_{exp} = 19.5$ mm filled with air at atmospheric pressure and at ambient temperature (Fig. 2). In order to reach high acoustic levels, the acoustic sound wave is generated by two loudspeakers (one at each end of the guide), tuned at the same frequency, that are driven in opposite phase. Those loudspeakers are connected to the main part of the guide via connecting tubes designed to avoid separation effects related to the singularities in change of section. The total length of the wave guide is $L_{exp} = 2.13$ m and the system is tuned at a frequency $f_{exp} = 240$ Hz so that to excite the third mode of the system. This choice of working frequency allows to perform measurements on a streaming cell away from singularities. The set-up geometry prevents the wave distortion even at high levels.

A wave generator provides the loudspeakers input signal, whose frequency and amplitude are controlled, as well as a trigger reference signal, used to synchronize the laser Doppler measurement system. Velocity measurements are started when the acoustic field is switched on in order to measure the transient evolution of the streaming velocity field. Those measurements are performed at different locations within the streaming cells. The laser system used and the processing of the laser signal are described in [28].

### 2.2. Numerical methods

For the numerical calculation, we consider a cylindrical tube of length $L$ and radius $R$, initially filled with air, considered as an ideal gas. An acoustic standing wave is initiated by shaking the tube in the axial direction ($z$) with given amplitude and frequency, so that a harmonic velocity law is imposed, $V(t) = V(t)\bar{e}_z$, with $V(t) = z_p\omega \cos(\omega t)$, $\omega$ being the angular frequency and $z_p$ the amplitude of the tube displacement.

The flow verifies the compressible Navier–Stokes equations expressed in the moving frame attached to the tube, with an additional forcing source term depending on $V$:

$$ \left\{ \begin{aligned} \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) &= 0 \\ \frac{\partial \rho \mathbf{v}}{\partial t} + \nabla \cdot (\rho \mathbf{v} \otimes \mathbf{v}) + \nabla p &= \nabla \cdot (\bar{\tau}) - \rho \frac{d\mathbf{V}}{dt} \end{aligned} \right. \quad (1) $$
<table><caption>**Table 1**
Parameters of different presented cases (N: numerical, E: experimental).</caption><thead><tr><th>Name</th><th>f (Hz)</th><th>R/δ<sub>v</sub></th><th>R/L</th><th>U<sub>max</sub> (m/s)(z = 0)</th><th>M</th><th>Re<sub>NL</sub></th></tr></thead><tbody><tr><td>N<sub>1</sub></td><td>20 000</td><td>50</td><td>0.09</td><td>0.5</td><td>1.45 × 10<sup>-3</sup></td><td>0.005</td></tr><tr><td>N<sub>2</sub></td><td>20 000</td><td>50</td><td>0.09</td><td>61.83</td><td>0.18</td><td>80.84</td></tr><tr><td>E<sub>2</sub></td><td>240</td><td>138</td><td>0.027</td><td>23</td><td>0.067</td><td>85.22</td></tr></tbody></table>

where *p* is the pressure, *ρ* is the density, **v** = *u*⃗*e*<sub>*z*</sub> + *v*⃗*e*<sub>*r*</sub> is the flow velocity, $\bar{\bar{\tau}} = -\frac{2}{3}\mu(\nabla \cdot \mathbf{v})\bar{\bar{I}} + 2\mu\bar{\bar{D}}$ the viscous stress tensor of a Newtonian fluid, $\bar{\bar{D}}$ the strain tensor, $\mu$ the dynamic viscosity which is supposed to be constant. In this formulation we have neglected the bulk viscosity since $\frac{M^2}{Sh^2Re_{ac}^2} \ll 1$ [30].

We consider an isentropic state law for the ideal gas, given by

$$p - p_0 = c_0^2(\rho - \rho_0), \qquad (2)$$

where $p_0$, $\rho_0$ and $c_0$ are respectively the constant initial pressure, density and sound velocity. Therefore the energy equation is decoupled from the other equations and we only solve the mass and momentum equations (Eqs. (1) and (2)), excluding thermal effects in our study.

No slip conditions are used on the solid boundary $r = R$, i.e. $\mathbf{v}(r = R, z) = 0$. On the axis $r = 0$, symmetry conditions are used, i.e. $\mathbf{v}(r = 0, z) = u(z)\vec{e}_z$ with $\frac{\partial u}{\partial r}(r = 0, z) = 0$. On the vertical boundaries $z = 0, L$, symmetry conditions are also used, i.e. $\mathbf{v}(r, z = 0, L) = v(r)\vec{e}_r$ with $\frac{\partial v}{\partial z}(r, z = 0, L) = 0$.

Due to the geometry, shock waves are generated in the simulations but it was checked that they are of small intensity and do not affect drastically the streaming flow [19,29]. The presence of harmonics will also be addressed in Section 6 of the present article.

The unsteady equations (Eq. (1)) are solved numerically using third order finite difference schemes [19,29]. The mean (streaming) flow is obtained from calculating a simple mean value for each physical quantity (velocity, pressure, density) over an acoustic period.

For simulations, the reference air pressure is $p_0 = 101 325$ Pa and the density is $\rho_0 = 1.2$ kg/m$^3$. The reference speed of sound is $c_0 = 343.82$ m/s. The air viscosity is $\mu = 1.795 \times 10^{-5}$ kg m$^{-1}$ s$^{-1}$. Resonant conditions are imposed, for which $L = \lambda/2$. Numerical values are $f_{num} = 20 000$ Hz, so that $\lambda/2 = 8.396$ mm and the acoustic period is $T_{num} = 50 \times 10^{-6}$ s. All results presented below are obtained using a regular mesh of rectangular cells composed of 500 points in the axial direction, and of $5 \times R/\delta_v$ points in the radial direction. The time step $\delta t$ is fixed equal to $8 \times 10^{-9}$ s.

### 3. Establishment of acoustic streaming

In order to study the establishment of the streaming flow generated by a sinusoidal standing wave independently from the transient of its acoustic source, the numerical code described in Section 2.2 is used in the following modified manner: for time stepping during the transient of acoustics, at the end of each acoustic period, the instantaneous fields are reinitialized by subtracting the mean flow over that period. Hence, the streaming flow is calculated here only once the acoustics is established, in order to dissociate the transient of streaming from the establishment of acoustics. Of course this cannot be done in experiments so we chose to discuss only experimental results for which the establishment of streaming has its own time scale, which is not the case of experimental results at low levels.

The different case parameters are summarized in Table 1. Subscript 1 is used to denote slow streaming and subscript 2 is used to denote fast streaming.

#### 3.1. Slow streaming

Figs. 3 and 4 address the case of slow streaming flow (for $N_1$). The established streaming flow pattern is very regular and corresponds to the representation in Fig. 1(a). Fig. 3 shows the time evolution of the streaming pattern: contours of axial streaming velocity are displayed (with color code) and streamlines are superimposed on the right half plane section. On Fig. 3(b) are also displayed two specific locations that will be used for further investigations: $I(r = R - 3\delta_v, z = L/4)$ representing the place of maximum inner streaming velocity and $O_c(r = 0, z = -L/4)$ representing the center of the plane section and the place of maximum outer streaming velocity in the slow regime. Note that the aspect ratio (height/length) of the displayed figures is equal to 0.25 which is much larger than the actual aspect ratio (height/length = 0.09) so that both inner and outer streaming are visible. One can observe that the establishment of the streaming flow occurs in several stages: The inner streaming cell is established first (Fig. 3(a)) and then does not evolve very much. The magnitude of the outer streaming velocity increases on a longer time scale (Fig. 3(b), (c) and (d)). This is in agreement with the results of previous numerical simulations [25].

The time evolution of the streaming velocity at some specific points can also be displayed. Fig. 4(a) shows the time evolution of the axial streaming velocity at $I$ (see Fig. 3(b)). Fig. 4(b) shows the time evolution of the axial streaming velocity
<figure><img src="image_107.png" /><figcaption>**Fig. 3.** Contours of the axial streaming velocity (with color codes in m/s) at different times during the transient evolution of $N_1$. Streamlines are superimposed on the right half plane section (real aspect ratio height/length=0.09). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

at the center of the plane section $O_c$ (see Fig. 3(b)). In agreement with Fig. 3 the inner streaming (Fig. 4(a)) is established much faster than outer streaming (Fig. 4(b)). As expected from Rayleigh calculations [2], the steady state values at $O_I$ and $O_c$ are roughly equal.

In the experiments, the establishment of acoustics at low amplitudes is very long (see Fig. 13 of [28]). In this case, the evolution of streaming is intrinsically linked to that of acoustics and the associated establishment time scales, determined by fitting an $\exp(-t/\tau)$ law to experimental data, are very close ($\tau_{ac} = 190$ s and $\tau_{st} = 233$ s). As will be shown in the next section, acoustics' transient is much shorter for high amplitudes ($\tau_{ac}$ is of the order of tenths of seconds) so that the establishment of streaming can be observed for different locations in the wave guide.
<figure><img src="image_108.png" /><figcaption>Fig. 4. Time evolution $N_1$. Axial streaming velocity at (a) $I$ and (b) $O_c$.</figcaption></figure>

## 3.2. Fast streaming

Figs. 5–7 address cases of fast streaming flow (for $N_2$ and $E_2$ respectively). The established streaming flow pattern is strongly distorted with respect to the slow streaming case, and two extra contra-rotating streaming cells are present near the acoustic velocity antinode (see Fig. 5(d)), in agreement with Fig. 1(b).

Fig. 5 shows the time evolution of the streaming patterns: contours of axial velocity are displayed on the left, and streamlines are displayed on the right. On Fig. 5(d) are also displayed two specific locations that will be used for further investigations: $O_e$ representing the place on the axis corresponding to the maximum velocity inside the extra contra-rotating cell and $O_{c-exp}$ representing the location on the axis, close to the position $z = L/4$, where measurements were performed to characterize the transient of the main outer streaming cell (both locations are on the right half side of the resonator). Here also the aspect ratio of the displayed figures is increased so that both inner and outer streaming are visible. One can observe that the establishment of the streaming flow occurs in several stages:

The inner streaming cell is established first, typically after $3 \times 10^{-3}$ s (see Fig. 5(a)) and then does not evolve very much. It is not as symmetrical as in the case of slow streaming (Fig. 3(a)), the cell centers being shifted towards the resonator closed ends (acoustic velocity node).

The magnitude of the outer streaming velocity increases on a longer time scale, and the main outer cell centers progressively shift towards the resonator closed ends (Fig. 5(b), (c) and (d)). Conversely to the slow streaming case, the maximum velocity of the outer streaming cell is smaller than that of inner streaming, as can be seen on the color code.

The emergence of the extra contra-rotating cell happens last (Fig. 5(d)).

On Fig. 6 are presented numerical time evolutions of the streaming velocity at $I$, $O_c$ (defined in the previous section, Fig. 3(b)) and at the place on the axis corresponding to the maximum velocity inside the extra contra-rotating cells $O_e$ (refer to Fig. 5(d)). First, this figure show that the time evolution of the inner streaming velocity at $I$ (Fig. 6(a)) is similar to the case of slow streaming (similar time scale).

Conversely in the core, the velocity increase for the location $O_c$ (Fig. 6(b)) depicts clearly an inflection point so that two characteristic time scales can be distinguished. Finally for the location $O_e$ (Fig. 6(c)) the velocity is first positive before becoming negative meaning that the contra-rotating cell appears after a certain time. A third time scale emerges from this time evolution. Note that the appearance of the extra cell does not modify the time evolution of streaming outside the extra cell (see Fig. 6(a) and (b)).

On Fig. 7 are presented experimental time evolutions for inner streaming at location $I$ (Fig. 7(a)), for outer streaming on the axis $r = 0$ at location $O_{c-exp}$ (Fig. 7(b)) and at the extra contra-rotating cell center $O_e$ (Fig. 7(c)) for the fast streaming case E2. Similar behavior to numerical results for equivalent positions can be observed on the experimental measurements and several distinct time scales are visible. The inner streaming is established first. Time evolution of outer streaming at $O_{c-exp}$ shows a distinct time scale. The magnitude of the contra-rotating cell maximum velocity is very small but the establishment of the extra streaming cell at $O_e$ is on a distinguished time scale from the last two.

In the experiment $E_2$ the established value in the inner cell ($U_s \approx 0.5$ m/s at $I$) is much greater than the established value in the core of the guide ($U_s \approx 0.01$ m/s at $O_{c-exp}$), as can be seen on Fig. 7. Similar observations were made on the simulation results in Fig. 5(b). Therefore, both experiments and numerical simulations show that nonlinear phenomena have very small influence on the streaming velocity in the boundary layer.

In summary, for slow streaming flow, only one time scale is necessary to describe the dynamics of streaming in the core of the waveguide. On the contrary for fast streaming flow, three different time scales were observed numerically and two experimentally in the core of the waveguide. At this stage the results seem to indicate that inner streaming flow is hardly modified by nonlinear effects (inertia included) at high levels. On the contrary the outer streaming is very strongly modified. This is in agreement with Menguy and Gilbert findings [15] which stated that inertia affects outer streaming but not inner streaming. However, it is necessary first to clearly identify the effect of inertia on the establishment of streaming. At this stage one cannot indeed assess whether inertial effects are truly responsible for the emergence of the new streaming cell.
<figure><img src="image_109.png" /><figcaption>Fig. 5. Contours of the axial streaming velocity (with color code) at different times during the transient evolution of N₂. Streamlines are superimposed on the right half plane section (real aspect ratio height/length = 0.09). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

In the next section, the equations governing the streaming flow are presented along with a scale analysis in order to exhibit governing dimensionless numbers and time scales for the establishment of inner and outer streaming flows. The expected time scales for the configuration under study are then calculated and compared with fits to the curves shown in the present section.

## 4. Scale analysis of slow and fast streaming flow

In order to analyze the dynamics of the streaming flow and particularly the competition between diffusion at the origin of streaming and convection associated to nonlinear inertial effects, a scale analysis based on the full equations governing
<figure><img src="image_110.png" /><figcaption>Fig. 6. Time evolution $N_2$. Axial streaming velocity at (a) $I$ (b) $O_c$ and (c) $O_e$.</figcaption></figure>

<figure><img src="image_111.png" /><figcaption>Fig. 7. Time evolution in the fast streaming case $E_2$. Axial streaming velocity at (a) $I$ (b) $O_{c-exp}$ and (c) $O_e$.</figcaption></figure>
the streaming flow is conducted below. Our approach is similar to that developed by Boufermel et al. [24] for computing the average mass transport velocity $\overline{\rho\dot{\nu}}$ in the slow streaming regime, as well as to that developed by Moudjed et al. [31] in the somehow different configuration of bulk streaming.

Focus is on outer streaming, whose transient and established regimes was shown to be drastically modified in the fast nonlinear regime.

## 4.1. Equations governing the streaming flow

Particle velocity $\mathbf{v}$, fluid density $\rho$ and pressure $p$ are decomposed into a fluctuating, periodic component, and a steady component (corresponding to the time average of the variable over the period, i.e. corresponding to the streaming flow)

$$ (\mathbf{v}, \rho, p) = (\mathbf{v}' + \bar{\mathbf{v}}, \rho' + \bar{\rho}, p' + \bar{p}), \quad (3) $$

the fluctuating component being denoted by a' superscript and the average operator by an overline. The time average of the fluctuating component is zero. Conversely to usual models of slow streaming, in the present approach no assumption is done on the magnitude of the fluctuation with respect to the steady component. The equations are therefore correct whether the streaming flow is slow or fast.

With this decomposition, the double and triple averaged products between any three variables $f, g, h$ become:

$$ \begin{cases} \overline{fg} = \bar{f}\bar{g} + \bar{f}'\bar{g}' \\ \overline{fgh} = \bar{f}\bar{g}\bar{h} + \bar{f}'\bar{g}'\bar{h} + \bar{f}'\bar{h}'\bar{g} + \overline{g'\bar{h}'\bar{f}} + \overline{f'\bar{g}'\bar{h}'} \end{cases} \quad (4) $$

The equations governing the streaming flow are obtained by averaging over the acoustic period the full compressible (isentropic) instantaneous Navier-Stokes equations (1) and the isentropic state law (2). Using the double and triple averaged products detailed in Eq. (4) then results in the following time averaged equation system:

$$ \left\{ \begin{aligned} \frac{\partial \bar{\rho}}{\partial t} + \nabla \cdot (\bar{\rho} \bar{\mathbf{v}}) &= -\nabla \cdot (\overline{\rho' \mathbf{v}'}) \\ \frac{\partial \bar{\rho}}{\partial t} \frac{\partial \bar{\mathbf{v}}}{\partial t} + \nabla \cdot (\bar{\rho} \bar{\mathbf{v}} \otimes \bar{\mathbf{v}}) + \nabla \bar{\rho} &= \nabla \cdot (\overline{\bar{\tau}}) + \bar{F}, \end{aligned} \right. \quad (5) $$

and $\bar{p} - p_0 = c_0^2(\bar{\rho} - \rho_0)$, where $\overline{\bar{\tau}}$ is the averaged viscous tensor and $\bar{F}$ is the excitation force given by:

$$ \bar{F} = -\nabla \cdot (\overline{\rho' \mathbf{v}' \otimes \mathbf{v}} + \bar{\mathbf{v}} \otimes \overline{\rho' \mathbf{v}'}) + \overline{\rho \overline{\mathbf{v}' \otimes \mathbf{v}'}} + \overline{\rho' \overline{\mathbf{v}' \otimes \mathbf{v}'}}. \quad (6) $$

Note that $\bar{F}$ is the sum of the source terms, representing the nonlinear contribution of the imposed acoustic field.

## 4.2. Scale analysis

The equation system (5) is written in dimensionless form, using standard fluid mechanics scalings, detailed in the following. As stated previously, focus is on the resonator core. In a cylindrical wave guide, streaming cells are toroidal of axial dimension $\lambda/4$ and radial dimension $R$ in the core of the resonator. Therefore the axial variable $z$ is scaled with $\lambda/4$ and the radial variable $r$ with $R$.

Mean velocity is scaled with a reference value of the axial streaming velocity $U_s$. Density is scaled with $\rho_0$. Pressure variations are scaled with $(\delta p)$, not specified at this stage of the scale analysis. The aim of the study is to determine the role of inertia on streaming in the fast regime. Therefore, time is scaled with the characteristic convection time scale in the axial direction $t_c = \lambda/(4U_s)$.

Scalings can then be summarized as follows:

$$ r^* = \frac{r}{R}, \quad z^* = \frac{z}{\lambda/4}, \quad \mathbf{v}^* = \frac{\bar{\mathbf{v}}}{U_s}, \quad t^* = \frac{t}{t_c}, \quad \rho^* = \frac{\bar{\rho}}{\rho_0}, \quad p^* = \frac{\bar{p} - p_0}{(\delta p)}. \quad (7) $$

Also, the fluctuation variables required to assemble the source terms are scaled with the acoustic amplitudes as follows:

$$ \mathbf{v}^{*\prime} = \frac{\mathbf{v}'}{U_{\max}}, \quad \rho^{*\prime} = \frac{\rho'}{\rho_{\max}}, \quad (8) $$

$\rho_{\max}$ being the amplitude of density fluctuations at acoustic antinode position in the standing wave. Whenever two different spatial derivatives are involved, the coefficient multiplying each dimensionless term is chosen with respect to the dominant spatial derivative under plane wave assumption: $R \ll \lambda$. In particular, scaling the viscous stress tensor yields a $\frac{\mu U_s}{R^2}$ factor, and scaling the divergence operator $\nabla \cdot (\mathbf{v}^*)$ yields a $\frac{U_s}{\lambda/4}$ coefficient. The following dimensionless equations are deduced from (5):

$$ \left\{ \begin{aligned} \frac{\rho_0}{t_c} \frac{\partial \rho^*}{\partial t^*} + \frac{\rho_0 U_s}{\lambda/4} \nabla \cdot (\rho^* \mathbf{v}^*) &= -\frac{\rho_{\max} U_{\max}}{\lambda/4} \nabla \cdot (\overline{\rho'^* \mathbf{v}'^*}) \\ \frac{\rho_0 U_s}{t_c} \frac{\partial \rho^* \mathbf{v}^*}{\partial t_c^*} + \frac{\rho_0 U_s^2}{\lambda/4} \nabla \cdot (\rho^* \mathbf{v}^* \otimes \mathbf{v}^*) + \frac{(\delta p)}{\lambda/4} \nabla p^* &= \frac{\mu U_s}{R^2} \nabla \cdot (\overline{\bar{\tau}^*}) + \bar{F}, \end{aligned} \right. \quad (9) $$
with

$$ \bar{F} = - \frac{\rho_{\max} U_{\max} U_s}{\lambda/4} \nabla \cdot (\overline{\rho'^* \mathbf{v}'*} \otimes \mathbf{v}^* + \mathbf{v}^* \otimes \overline{\rho'^* \mathbf{v}'*} + \frac{\rho_0 U_{\max}}{\rho_{\max} U_s} \overline{\rho^* \mathbf{v}'*} \otimes \mathbf{v}'^* + \frac{U_{\max}}{U_s} \overline{\rho'^* \mathbf{v}'*} \otimes \mathbf{v}^*). \quad (10) $$

Then, multiplying the mass equation with $\lambda/(4\rho_0U_s)$ and the momentum equation with $\lambda/(4\rho_0U_s^2)$ and using $t_c = \lambda/(4U_s)$ one can obtain:

$$ \left\{ \begin{aligned} \frac{\partial \rho^*}{\partial t^*} + \nabla \cdot (\rho^* \mathbf{v}^*) &= - \frac{\rho_{\max}}{\rho_0} \frac{U_{\max}}{U_s} \nabla \cdot (\overline{\rho'^* \mathbf{v}'*}) \\ \frac{\partial \rho^* \mathbf{v}^*}{\partial t^*} + \nabla \cdot (\rho^* \mathbf{v}^* \otimes \mathbf{v}^*) + \frac{(\delta p)}{\rho_0 U_s^2} \nabla p^* &= \frac{1}{\bar{R} e} \nabla \cdot (\bar{\bar{\tau}}^*) + \bar{F}^*, \end{aligned} \right. \quad (11) $$

where the dimensionless excitation force, gathering all the source terms is expressed as:

$$ \bar{F}^* = - \frac{\rho_{\max} U_{\max}}{\rho_0 U_s} \nabla \cdot (\overline{\rho'^* \mathbf{v}'*} \otimes \mathbf{v}^* + \mathbf{v}^* \otimes \overline{\rho'^* \mathbf{v}'*}) \\ - \frac{U_{\max}^2}{U_s^2} \nabla \cdot (\rho^* \overline{\mathbf{v}^* \otimes \mathbf{v}^*}) - \frac{\rho_{\max}}{\rho_0} \frac{U_{\max}^2}{U_s^2} \nabla \cdot (\overline{\rho'^* \mathbf{v}^* \otimes \mathbf{v}^*}). \quad (12) $$

The system equations (Eq. (11)) is the dimensionless form of the general governing streaming equations (Eq. (5)). The Reynolds number $Re = \frac{4U_s R^2}{\lambda \nu}$ that appears in the momentum equation is the Reynolds number characterizing the mean flow in a thin tube, defined as the ratio between the diffusion time scale in the radial direction $t_d = R^2/\nu$ and the convection time scale in the axial direction $t_c = \lambda/(4U_s)$, $Re = t_d/t_c$. Note that $t_d$ is equivalent (with a different length reference value) to the characteristic time for streaming development defined by Amari et al. in [32] and used to obtain the transient linear streaming equations. Under the assumptions that acoustic perturbations are small (i.e. of order $M = \frac{U_{\max}}{c_0} \ll 1$) and the streaming field is of order $M^2$ ($U_s = MU_{\max}$), this Reynolds number can be rewritten as $Re = \frac{4 M^2}{\pi \bar{R} \nu^2} = \frac{4}{\pi} Re_{NL}$, $Re_{NL}$ being the commonly introduced nonlinear Reynolds number.

Under these assumptions, the Reynolds number $Re$ and the nonlinear Reynolds number $Re_{NL}$ are thus of the same order of magnitude. However, since the Reynolds number $Re$ is defined directly with the value of the maximum streaming velocity, it describes more accurately the influence of inertia in the nonlinear regime in the core of the wave guide than $Re_{NL}$. Indeed, previous numerical and experimental studies [29] have shown that in the nonlinear regime, the maximum value of the axial streaming velocity on the tube axis is much smaller than slow streaming expectations $3U_{\max}^2/(8c)$. Therefore the nonlinear Reynolds number $Re_{NL}$ overestimates the influence of inertia on the streaming flow. For example, in the numerical case $N_2$, the Reynolds number based on the maximum value of $U_s$ on the axis is $Re = 17.1$ while the corresponding nonlinear Reynolds number is $Re_{NL} = 80.84$. Nonetheless $Re$ cannot be estimated before the simulations/experiments whereas $Re_{NL}$ can be calculated from the acoustics and thus remains useful.

Another Reynolds number used in the literature to distinguish linear from nonlinear regime was the streaming Reynolds number used by Stuart et al. [11] defined as $Re_s = \frac{U_s \lambda}{2 \pi \nu}$. It can be expressed as $Re_s = Re \frac{\lambda^2}{R^2}$ and also has the drawback to overestimate the nonlinear inertial effects on Rayleigh streaming.

In the case $Re \ll 1$ (which also means that convection time is much greater than diffusion time) and under the assumption of small perturbations ($M = \frac{U_{\max}}{c_0} = \frac{U_s}{U_{\max}} \ll 1$), the inertial term and the time depending term can both be neglected with respect to the viscous term. The magnitude of source terms is then estimated:

$$ \frac{\rho_{ac}}{\rho_0} \frac{U_{ac}}{U_s} \nabla \cdot (\overline{\rho'^* \mathbf{v}'*} \otimes \mathbf{v}^* + \mathbf{v}^* \otimes \overline{\rho'^* \mathbf{v}'*}) = O\left(M \cdot \frac{1}{M}\right) = O(1), \quad (13) $$

$$ \frac{U_{\max}^2}{U_s^2} \nabla \cdot (\rho^* \overline{\mathbf{v}^* \otimes \mathbf{v}^*}) = O\left(\frac{1}{M^2}\right), \quad (14) $$

$$ \frac{\rho_{\max}}{\rho_0} \frac{U_{\max}^2}{U_s^2} \nabla \cdot (\overline{\rho'^* \mathbf{v}'*} \otimes \mathbf{v}'^*) = O\left(M \cdot \frac{1}{M^2}\right) = O\left(\frac{1}{M}\right). \quad (15) $$

$\bar{F}^*$ is then of order $\frac{1}{M^2}$. The pressure gradient must have the same order of magnitude as $\bar{F}^*$ which yields also the pressure variations scale $(\delta p) = \rho_0 U_{\max}^2 = p_0 M^2$. The magnitude of the source term in the mass equation is estimated as previously: $\frac{\rho_{\max} U_{\max}}{\rho_0 U_s} \nabla \cdot (\overline{\rho'^* \mathbf{v}'*}) = O(1)$. Finally, the steady linear equations for the streaming flow (dimensional form) are then deduced from the first order approximation of the momentum equation and (5) gives:

$$ \left\{ \begin{aligned} \nabla \cdot (\overline{\rho \mathbf{v}}) &= -\nabla \cdot (\overline{\rho' \mathbf{v}'}) \\ \nabla \bar{p} &= \nabla \cdot (\bar{\bar{\tau}}) - \bar{F}, \end{aligned} \right. \quad (16) $$

where $\bar{F} \approx \nabla \cdot (\overline{\rho \mathbf{v}' \otimes \mathbf{v}'})$.

As the $Re$ number is increased i.e. $Re \gtrsim 1$, inertial and time depending terms cannot be neglected anymore and the complete unsteady form of streaming equations (5) with all the source terms should be used. In this case the diffusion time scale becomes larger than the convection time scale, thus the overall process is faster and the steady state is reached earlier.
<table><caption>**Table 2**
Reynolds numbers of different presented cases (N: numerical, E: experimental).</caption><thead><tr><th>Name</th><th>U<sub>s</sub> (at I) (m/s)</th><th>U<sub>s</sub> (at O<sub>c</sub>) (m/s)</th><th>Re<sub>NL</sub></th><th>Re(I)</th><th>Re(O<sub>c</sub>)</th></tr></thead><tbody><tr><td>N<sub>1</sub></td><td>2.42 × 10<sup>-4</sup></td><td>2.63 × 10<sup>-4</sup></td><td>0.005</td><td>2.24 × 10<sup>-3</sup></td><td>2.43 × 10<sup>-3</sup></td></tr><tr><td>N<sub>2</sub></td><td>3.18</td><td>1.84</td><td>80.84</td><td>29.44</td><td>17</td></tr><tr><td>E<sub>2</sub></td><td>0.5</td><td>0.01</td><td>85.22</td><td>35.26</td><td>0.7</td></tr></tbody></table>

### 4.3. Reynolds numbers and time scale values for slow and fast outer streaming

Table 2 displays the values of the axial streaming velocity at points I or O<sub>c</sub>, and of the corresponding Reynolds number of previously described cases N<sub>1</sub>, N<sub>2</sub> and E<sub>2</sub>.

The previous analysis can also be used to estimate quantitatively the expected time scales of the streaming time evolutions for outer streaming, in the slow and fast streaming regime, that was described qualitatively in Section 3.

In the slow streaming case, diffusion of momentum in the radial direction is expected to be dominant. The convection time variations take place over a very long scale and they are very small, therefore only time variations over a diffusion time scale are significant. The associated characteristic time scale for the establishment of outer streaming is $t_d = R^2/v$, which, using the parameter values in Table 2 for N<sub>1</sub>, amounts to $t_d = 0.026$ s. Taking into account the fact that the streaming velocity at O<sub>c</sub> in Fig. 4(b) is not established in 0.01 s, this is a coherent result.

Based on known analytical solutions of the unsteady pure diffusion equation, if diffusion were dominant the time evolution curves should be in the form $u(t) = u_∞(1 - e^{-t/τ})$, where the time constant $τ$ is the time needed for the value $u$ to reach 63% of the established value $u_∞$ and the time constant needed to achieve 99.999% of the established value is equal to $10τ$. The time evolution curve at O<sub>c</sub> in Fig. 4(b) is actually very well fitted with such a diffusion curve, and the quantitative fit gives $τ ≈ 0.0015$ s and therefore $10τ = 0.015$ s, which is the same order of magnitude as $t_d$.

In the fast streaming case N<sub>2</sub>, in the core, transverse diffusion is expected to compete with convection in the axial direction. The value of the theoretical diffusion time scale $t_d$ is still the same, since it does not depend on the excitation level ($t_d = 0.026$ s). The first velocity increase in the time evolution curve at O<sub>c</sub> (Fig. 6(b)) can be fitted by assuming dominant diffusion and the obtained characteristic time scale is the same as in the slow streaming case. Therefore transverse diffusion is dominant at first in the establishment of the outer streaming, just like in the slow streaming case. Assuming next that the second velocity increase in the time evolution curve at O<sub>c</sub> (Fig. 6(b)) is related to convection, the theoretical time scale associated with convection along the length $λ/4$ with velocity $U_s$ at the location I is $t_c = λ/(4U_s) = 1.6 × 10^{-3}$ s based on the parameters in Table 2. This increase is 20 times faster than that of the diffusion process. A linear fit of the second velocity increase at O<sub>c</sub> on Fig. 6(b) gives $t_{cfit} ≈ 10^{-3}$ s, which is of the same order as $t_c$. This shows that the second velocity increase is indeed related to convection in the axial direction. Note that convective transport is much faster than diffusive transport and that the steady state is reached earlier than in the slow streaming case N<sub>1</sub>, in agreement with the analysis given at the end of Section 4.2.

In the experimental case E<sub>2</sub>, Fig. 7 shows the evolution of streaming at the three different locations I, O<sub>c-exp</sub> and O<sub>e</sub>. Two separate slopes in the time evolution of streaming at O<sub>c-exp</sub> are also observed, but these are not due to the competition between diffusion and convection as was the case for numerical results. Actually our experimental set up does not allow to observe this competition: the theoretical value of the radial diffusion time scale is $t_d = 24$ s and the convection time scale in the axial direction calculated for E<sub>2</sub> based on the maximum value of U<sub>s</sub> at I location (read in Fig. 7(a)) is $t_{cE2} = 0.71$ s. It can be observed that acoustics requires about 3 min to be established. The transient of acoustics is therefore too long to allow to separately identify the establishment of streaming. Moreover one measurement point actually corresponds to an integration of results of particle velocity measurements over several seconds: indeed in order to obtain converged results of the particle velocity estimate, several thousands of seeding particles need to have crossed the measuring volume [27,33]. The associated measuring time depends on the data rate associated with the seeding density; in our set up this process can take as long as 10 s. Therefore it is impossible to distinguish time scales as short as $t_d$ or $t_c$. The evolution of streaming shown in Fig. 7 can be supposed to result from heat transfer effects. Indeed measurements give an estimate of the mean velocity, that is equal to streaming velocity only if no heat convection occurs. In the case of E<sub>2</sub> experiments, a 4 K temperature difference in the guide wall temperature was measured across a distance of $λ/4$. The temperature establishment was found to last for 20 min, probably due to the time needed to reach thermal equilibrium between thermoacoustic, conduction and convection heat fluxes. It can thus be assumed that the measured time evolution of mean velocity is not the only result of acoustic streaming. In the following, only the numerical results will be considered.

In the next section, the driving effect of inner streaming on the outer streaming will be analyzed. Indeed, an analogy with the lid-driven cavity in a cylindrical geometry is presented in order to identify the physical meanings of the characteristic times.

## 5. Lid-driven cavity analogy

Based on the observation that inner streaming drives outer streaming with competing diffusion and convection, and in order to mimic the establishment of outer streaming, we introduce an equivalent fluid-mechanics problem. We therefore consider the establishment of flow in a lid-driven cavity, but with a cylindrical geometry, and assume that the “lid” is the tube envelope where the inner streaming velocity is imposed.
<figure><img src="image_112.png" /><figcaption>Fig. 8. Schematics of the lid-driven cylindrical flow problem: geometry and boundary conditions.</figcaption></figure>

<figure><img src="image_113.png" /><figcaption>Fig. 9. (a) Time evolution (—) of $U(z = L/4, r = 0)$ in lid-driven cavity for $U_w$ from $N_1$, and (—) of $U_s$ at $O_c$ from Fig. 4. (b) Same with both curves normalized to reach the same steady-state value. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

## 5.1. Description of the problem

The problem being axisymmetric, the geometry is defined by $z \in [0, L/2]$ and $r \in [0, R - 3\delta_v]$. The viscous flow is governed by the following equations, which are the same as in Eq. (1) without the acoustic source terms:

$$ \left\{ \begin{aligned} \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) &= 0 \\ \frac{\partial \rho \mathbf{v}}{\partial t} + \nabla \cdot (\rho \mathbf{v} \otimes \mathbf{v}) + \nabla p &= \nabla \cdot (\bar{\tau}) \end{aligned} \right. \qquad (17) $$

This system is completed with the isentropic state law Eq. (2). Symmetry conditions are imposed at the vertical walls $z = 0, L/2$ and on the horizontal axis $r = 0$. The boundary condition at the wall $r = R - 3\delta_v$ (the “lid”) is that the velocity is only axial. Its value $U_w(z)$ is equal to the streaming axial velocity profile $U_s(z)$ along the tube at $r = R - 3\delta_v$ given by direct numerical calculations described in the previous sections (see Fig. 8). There is initially no flow in the cavity. The numerical procedure is the same as that described in Section 2.2. The simulations show that the flow is almost incompressible (density variations are less than 0.006%).

## 5.2. Results

A first simulation (low level) is run with $U_w$ corresponding to the established $U_s$ at $r = R - 3\delta_v$ for $N_1$ case (slow streaming). Fig. 9(a) shows (in red) the time evolution of the axial velocity $U$ at $r = 0$ and $z = L/4$ (center of the mimicked streaming cell). The time evolution of the streaming velocity at $O_c$ (same location) in Fig. 4 is reproduced in blue dots. Fig. 9(b) shows both curves normalized so that they reach the same steady state value. Both curves are very similar with similar time scales, which confirms that transverse diffusion of momentum is the dominant physical phenomenon for the establishment of the slow streaming pattern.

A simulation at higher level is run with $U_w$ corresponding to the established $U_s$ at $r = R - 3\delta_v$ for $N_2$ case (fast streaming). Fig. 10(a) shows (solid red line) the time evolution in the lid-driven cavity of the axial velocity $U$ at $r = 0$ and $z = L/4$. The time evolution of the streaming velocity at $O_c$ (same location) given by Fig. 6(b) is reproduced in blue dots. Both curves are very similar in shape, although the maximum value obtained for the streaming curve is about half of that in the lid-driven cavity. Fig. 10(b) shows the same curves normalized to reach the same steady-state. Both curves are almost superimposed, which shows that the time scales are the same and that the change of slope in the curve is indeed related to the inertial effect, when convection becomes strong enough to compete efficiently with diffusion and to modify the time evolution. However, as stated above the established streaming velocity is much smaller than that obtained in the lid-driven cavity. This suggests that the inertial effects of the lid-driven cavity model are not sufficient to correctly model the change in the streaming flow field.

Fig. 11 shows several snapshots of contours of the axial velocity in the lid-driven cavity during the transient evolution. Streamlines are superimposed on the right half plane section (the simulations are performed on a length $L/2$ and the solution
<figure><img src="image_114.png" /><figcaption>Fig. 10. (a) Time evolution (—) of $U(z = L/4, r = 0)$ in lid-driven cavity for $U_w$ from $N_2$, and (—) of $U_s$ at $O_c$ from Fig. 6. (b) Same with both curves normalized to reach the same steady-state value. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

in the right half plane section is constructed by symmetry). Comparison with Fig. 5 shows that the time evolution of patterns are very similar which confirms that the transient evolution of the streaming patterns comes from competing diffusion and convection in a manner analogous to the lid-driven cavity. However, an important difference exists: there is no extra streaming cell in the lid-driven cavity.

Fig. 12 shows axial profiles of $U$ on the axis $r = 0$ for corresponding times of Fig. 11. The axial profile imposed at the wall $r = R - 3\delta_v$ is also shown. At early times the axial profile of $U$ is similar to the profile of $U_w$, which confirms that momentum diffusion in the radial direction is dominant. The distortion observed in the last profile is due to convection in the axial direction.

A similar simulation can be conducted using the experimental configuration of $E_2$. Fig. 13 shows the time evolution of the axial velocity $U$ at $r = 0$ and $z = L/4$ in the lid-driven cavity (the experimental inner streaming axial profile is imposed at the lid). The associated time scales are of the same order as $t_d$ and $t_c$ calculated in Section 4.3, and the transient behavior of the flow is analogous to that observed in Fig. 10. As mentioned earlier, this evolution cannot be highlighted experimentally.

This analysis shows that although diffusion and convection do compete in the establishment of the streaming flow in a similar manner as in the lid-driven cavity problem, this is not enough to explain the emergence of the new streaming cell. Therefore, at this stage, we can conclude that the lid-cavity model is not complete enough to describe the mutation of streaming at high levels.

In the next section, the complete streaming equations (5) are solved numerically, with a mono-frequency acoustic source used to compute the different source terms. The effect of increasing the amplitude of the acoustic flow field is analyzed.

## 6. Averaged model simulation for the streaming flow

The high-level mono-frequency acoustic field source used in this section is obtained in two steps. Firstly a DNS simulation at very low-level amplitude is run until a steady streaming flow is obtained. The resulting acoustic field, calculated by subtracting the streaming flow from the total velocity field ($v' = v - \bar{v}$, $\rho' = \rho - \bar{\rho}$), contains insignificant level of harmonics. Then this acoustic field ($v'$, $\rho'$) is multiplied by a coefficient $\alpha$ chosen to obtain high values of the nonlinear Reynolds number $Re_{NL}$. This acoustic field is used to compute the source terms so that the following averaged Navier-Stokes equations are solved for the unknown streaming field ($v_s$, $\rho_s$, $p_s$):

$$ \left\{ \begin{aligned} \frac{\partial \rho_s}{\partial t} + \nabla \cdot (\rho_s \mathbf{v}_s) &= -\alpha^2 \nabla \cdot (\bar{\rho}' \mathbf{v}') \\ \frac{\partial \rho_s \mathbf{v}_s}{\partial t} + \nabla \cdot (\rho_s \mathbf{v}_s \otimes \mathbf{v}_s) + \nabla p_s &= \nabla \cdot (\bar{\tau}_s) + \bar{F}_1, \end{aligned} \right. \qquad (18) $$

with $\bar{F}_1$ defined by

$$ \bar{F}_1 = -\nabla \cdot (\alpha^2 \overline{\rho' \mathbf{v}'} \otimes \mathbf{v}_s + \alpha^2 \mathbf{v}_s \otimes \overline{\rho' \mathbf{v}'} + \alpha^2 \rho_s \overline{\mathbf{v}'} \otimes \mathbf{v} + \alpha^3 \overline{\rho' \mathbf{v}'} \otimes \overline{\mathbf{v}'}) \qquad (19) $$

and $p_s - p_0 = c_0^2(\rho_s - \rho_0)$. The boundary conditions are no slip on the outside wall parallel to the horizontal symmetry axis, and symmetry conditions on the vertical boundaries and horizontal axis.

These equations, very similar to Eqs. (5), are solved with the same numerical method as presented in Section 2.2.

This AMS approach can be used to study the influence of inertial effects on the streaming flow for a monofrequency high amplitude acoustic forcing field. This is equivalent to the strategy followed in [15]. In this way we are able to exclude the influence on the streaming flow of higher harmonics in the acoustic signal. This numerical approach is used here to investigate nonlinear Reynolds numbers $Re_{NL}$ ranging from 0.002 up to 94 (corresponding to $Re$ ranging from 0.0025 to 30.81).

Fig. 14 displays the axial streaming velocity component $U_s$ obtained with the AMS simulation scaled by the Rayleigh streaming maximum velocity $3 \times U_{\text{max}}^2/(8c_0)$ along the scaled central axis $z/L$ for $R/\delta_v = 40$ and $0.002 \le Re_{NL} \le 94$. For
<figure><img src="image_115.png" /><figcaption>**Fig. 11.** Contours of the axial velocity (with color code) at different times during the transient evolution in the lid-driven cavity for $U_w$ from $N_2$. Streamlines are superimposed on the right half plane section. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

increasing $Re_{NL}$ numbers, the sine function associated to slow streaming becomes steeper next to the resonator center, at the acoustic velocity antinode location, and not next to the resonator ends (acoustic velocity nodes locations). This behavior is in agreement with the predictions of the nonlinear model developed by Menguy and Gilbert [15] at moderate nonlinear Reynolds numbers $Re_{NL} = O(1)$. The direction of distortion is opposite to that observed in the full DNS simulations and experiments [28,29].

In Fig. 15, the dimensionless streaming velocity is plotted along the axis for DNS, AMS and Lid-driven cavity (LDC) simulations, for slow-streaming ($N_1$) and fast-streaming ($N_2$). For slow-streaming, all three simulations give analogous
<figure><img src="image_116.png" /><figcaption>**Fig. 12.** Axial profiles of $U$ at $r = 0$ (in blue) for corresponding times expressed in $10^{-3}$ s and at $r = R - 3\delta_v$ (in red) in the lid-driven cavity for $U_w$ from $N_2$. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_117.png" /><figcaption>**Fig. 13.** Time evolution of $U(z = L/4, r = 0)$ in lid-driven cavity for $U_w$ from $E_2$.</figcaption></figure>

<figure><img src="image_118.png" /><figcaption>**Fig. 14.** AMS results for streaming velocity scaled by $3 \times U_{\text{max}}^2/(8c_0)$ along the axis, $R/\delta_v = 40$.</figcaption></figure>

profiles even though the LDC streaming velocity is slightly smaller. This is because the lid-cavity model does not include the acoustic source terms.

For fast-streaming, the AMS result shows the same behavior as the lid-driven cavity. Once again the LDC streaming velocity is slightly smaller, because that model does not include the acoustic source terms. For the AMS simulations, the associated axial velocity contours and streamlines are displayed in Fig. 16. They are similar to those obtained for the lid driven cavity flow, in the steady state (see Fig. 11(d)). On the contrary, the DNS streaming velocity on the axis is globally much smaller, and shows the presence of new streaming cells.

Both the lid-driven cavity and the AMS approaches show that when the intensity of a monofrequency wave is increased, the streaming flow becomes nonlinear because of inertial effects, but the distortion of the Rayleigh streaming profile is different than in DNS simulations or experiments.

This shows that in order to explain the mutation of the streaming flow at high levels, other effects must be analyzed. The only effects that have not been considered yet are the consequences of nonlinear interaction between the streaming flow and acoustics, i.e. the modification of the acoustic flow by the streaming flow. This is currently under study.

## 7. Conclusion

In this paper the effect of inertia on Rayleigh streaming was investigated numerically and experimentally.

Time evolution of streaming patterns were analyzed. Based on an analysis of the equations driving the streaming flow, two time scales were identified: the diffusion and the convection time scales. The ratio of these time scales is a Reynolds
<figure><img src="image_119.png" /><figcaption>**Fig. 15.** Dimensionless streaming velocity along the axis. Comparison of DNS (solid lines), AMS (dashed lines) and Lid-Driven Cavity, LDC (dashed-dotted lines) results, $R/\delta_v = 50$, for slow-streaming (in blue) and fast-streaming (in red). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_120.png" /><figcaption>**Fig. 16.** AMS results, $R/\delta_v = 50$. Isolines of the axial streaming velocity colored by its value and the streamlines of the streaming flow. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

number $Re = \frac{4U_s R^2}{\lambda \nu}$ that characterizes the streaming flow and separates the linear from the nonlinear regime. This number is proportional to the nonlinear Reynolds number $Re_{NL}$ only when the streaming flow is of order $M^2c_0$.

The time evolution of outer streaming was mimicked by analyzing a lid-driven flow. It was confirmed that in a first stage, the outer fast streaming cells are created due to momentum diffusion in the radial direction and that in a second stage, the streaming cells are distorted due to convection in the axial direction. But no extra cell is created in the lid-driven cavity.

Nonlinear inertial effects on the established streaming flow were further investigated using a numerical approach that solves the full streaming flow equations with an imposed monofrequency acoustic wave source computed from the DNS. These simulations as well as the lid-driven cavity simulations show similar results: The streaming patterns in the core of the waveguide are different than those observed in DNS simulations and experiments.

Therefore inertial effects cannot be considered as the leading phenomenon to explain the mutation of streaming at high acoustic levels. Nonlinear interactions between the streaming flow and acoustics have to be considered, which will be the scope of further studies.

## Acknowledgment

The authors wish to thank J. Gilbert for many fruitful discussions on this work.

## References

[1] G.W. Swift, D.L. Gardner, S. Backhaus, Acoustic recovery of lost power in pulse tube refrigerators, J. Acoust. Soc. Am. 105 (1999) 711–724.

[2] L. Rayleigh, On the circulation of air observed in Kundt's tubes, and on some allied acoustical problems, Philos. Trans. R. Soc. 175 (1884) 1–21.

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

[4] W.L. Nyborg, Acoustic streaming, in: W.P. Mason (Ed.), Physical Acoustics, Vol. 2B, Academic Press, New York, 1965, pp. 265–331.

[5] P.J. Westervelt, The theory of steady rotational flow generated by sound fields, J. Acoust. Soc. Am. 25 (1953) 60–67.

[6] O.V. Rudenko, S.I. Soluyan, Theoretical Fundations of Non- Linear Acoustics, English Translation, Plenum, New York, 1977, p. 187–211.

[7] N. Rott, The influence of heat conduction on acoustic streaming, Z. Angew. Math. Phys. 25 (1974) 417–421.

[8] Q. Qi, The effect of compressibility on acoustic streaming near a rigid boundary for a plane traveling wave, J. Acoust. Soc. Am. 94 (1993) 1090–1098.

[9] H. Bailliet, V. Gusev, R. Raspet, R.A. Hiller, Acoustic streaming in closed thermoacoustic devices, J. Acoust. Soc. Am. 110 (2001) 1808–1821.

[10] M.F. Hamilton, Y.A. Ilinskii, E.A. Zabolotskaya, J. Acoust. Soc. Am. 114 (6) (2003) 3092–3101.

[11] J.T. Stuart, Double boundary layer in oscillatory viscous flow, J. Fluid Mech. 24 (1966) 673–987.

[12] J. Lighthill, Acoustic streaming, J. Sound Vib. 61 (1978) 391–418.

[13] V.E. Gusev, O.V. Rudenko, Nonsteady quasi-one-dimensional acoustic streaming in unbounded volumes with hydrodynamic nonlinearity, Sov. Phys. - Acoust. 25 (1979) 493–497.
[14] N. Riley, Streaming from a cylinder due to an acoustic source, J. Fluid Mech. 180 (1987) 319–326.

[15] L. Menguy, J. Gilbert, Non-linear acoustic streaming accompanying a plane stationary wave in a guide, Acta Acust. 86 (2000) 249–259.

[16] A. Gopinath, A.F. Mills, Convective heat transfer due to acoustic streaming across the ends of a Kundt tube, J. Heat Transfer 116 (1994) 47–53.

[17] M. Kawahashi, M. Arakawa, Nonlinear phenomena induced by finite-amplitude oscillation of air column in closed duct, JSME Int. J. 39 (1996) 280–286.

[18] M.K. Aktas, B. Farouk, Numerical simulation of acoustic streaming generated by finite-amplitude resonant oscillations in an enclosure, J. Acoust. Soc. Am. 116 (5) (2004) 2822–2831.

[19] V. Daru, D. Baltean-Carlès, C. Weisman, P. Debesse, G. Gandikota, Two-dimensional numerical simulations of nonlinear acoustic streaming in standing waves, Wave Motion 50 (2013) 955–963.

[20] A.A. Gubaidullin, A.V. Yakovenko, Effects of heat exchange and nonlinearity on acoustic streaming in a vibrating cylindrical cavity, J. Acoust. Soc. Am. 137 (2015) 3281–3287.

[21] A.A. Gubaidullin, A.V. Pyatkova, Acoustic streaming with allowance for heat transfer, Acoust. Phys. 62 (3) (2016) 300–305.

[22] C. Scalo, S.K. Lele, L. Hesselink, Linear and nonlinear modelling of a theoretical travelling-wave thermoacoustic heat engine, J. Fluid Mech. 766 (2015) 368–404.

[23] M. Cervenka, M. Bednarik, Variety of acoustic streaming in 2D resonant channels, Wave Motion 66 (2016) 21–30.

[24] A. Boufermel, N. Joly, P. Lotton, M. Amari, V. Gusev, Velocity of mass transport to model acoustic streaming: Numerical application to annular resonators, Acta Acust. United Acust. 97 (2) (2011) 219–227.

[25] S. Boluriaan, P.J. Morris, Numerical simulation of acoustic streaming in high amplitude standing wave, in: AAIA Paper 2003-3152 (9th AIAA/CEAS Aeroacoustics Conference and Exhibit, 12-14 May 2003, Hilton Head, South Carolina).

[26] M. Thompson, A. Atchley, M. Maccarone, Influences of a temperature gradient and fluid inertia on acoustic streaming in a standing wave, J. Acoust. Soc. Am. 117 (2005) 1839–1849.

[27] S. Moreau, H. Bailliet, J.-C. Valière, Measurements of inner and outer streaming vortices in a standing waveguide using laser doppler velocimetry, J. Acoust. Soc. Am. 123 (2008) 640–664.

[28] I. Reyt, H. Bailliet, J.C. Valière, Experimental investigation of acoustic streaming in a cylindrical wave guide up to high streaming Reynolds numbers, J. Acoust. Soc. Am. 135 (2014) 27–37.

[29] I. Reyt, V. Daru, H. Bailliet, S. Moreau, J.C. Valière, D. Baltean-Carlès, C. Weisman, Fast acoustic streaming in standing waves: Generation of an additional outer streaming cell, J. Acoust. Soc. Am. 134 (3) (2013) 1791–1801.

[30] C. Olivier, G. Penelet, G. Poignand, J. Gilbert, P. Lotton, Weakly nonlinear propagation in thermoacoustic engines: A numerical study of higher harmonics generation up to the appearance of shock waves, Acta Acust. United Acust. 101 (5) (2015) 941–949.

[31] B. Moudjed, V. Botton, D. Henry, H. Ben Hadid, J.-P. Garandet, Scaling and dimensional analysis of acoustic streaming jets, Phys. Fluids 26 (2014) 093602.

[32] M. Amari, V. Gusev, N. Joly, Temporal dynamics of the sound wind in acoustitron, Acta Acust. United Acust. 89 (2003) 1008–1024.

[33] J.-C. Valière, S. Moreau, H. Bailliet, Development of laser techniques for acoustic boundary layer measurements. Part I: LDV signal processing for high acoustic displacements, Acta Acust. United Acust. 95 (2009) 585–594.
D. Baltean-Carlès, V. Daru, C. Weisman, S. Tabakova, H. Bailliet

An unexpected balance between outer Rayleigh streaming sources
J. Fluid Mech. 867, 985-1011, 2019.
# An unexpected balance between outer Rayleigh streaming sources

D. Baltean-Carlès<fnref n="1" />,<fnref n="3" />, V. Daru<fnref n="2" />,<fnref n="3" />, C. Weisman<fnref n="1" />,<fnref n="3" />,†, S. Tabakova<fnref n="4" /> and H. Bailliet<fnref n="5" />

(Received 25 July 2018; revised 20 December 2018; accepted 4 February 2019)

Acoustic streaming generated by a plane standing wave between two infinite plates or inside a cylindrical tube is considered, under the isentropic flow assumption. A two-dimensional analysis is performed in the linear case of slow streaming motion, based on analytical formal solutions of separate problems, each associated with a specific source term (Reynolds stress term). In order to obtain these analytical solutions, a necessary geometrical hypothesis is that $(R/L)^2 \ll 1$, where R and L are the guide half-width (or radius) and length. The effect of the two source terms classically taken into account is quantified in order to derive the dependence of the maximum axial streaming velocity on the axis as a function of the ratio $R/\delta_\nu$, where $\delta_\nu$ is the acoustic boundary layer thickness. The effect of two other source terms that are usually neglected, is then analysed. It is found that one of these terms can generate a counter-rotating streaming flow. While negligible for very narrow guides, this term can become important for some values of the aspect ratio $L/R$.

**Key words:** acoustics, boundary layers, low-Reynolds-number flows

## 1. Introduction

Rayleigh streaming generated by a plane standing wave between two infinite plates or inside a cylindrical tube is considered in the present paper. This acoustic wind is a mean flow generated by Reynolds stresses in the Stokes boundary layer formed along the solid walls of an acoustic guide (with no other mean flow and a no-slip boundary condition on the guide walls). It is usually referred to as Rayleigh streaming because Lord Rayleigh (1884) was the first to develop an analytical solution that describes the steady vortices generated in the core of the fluid in a wide channel formed by two

<footnote n="1">¹Sorbonne-Université, Faculté des Sciences et Ingénierie, UFR d'Ingénierie, 4 Place Jussieu, 75005 Paris, France</footnote>

<footnote n="2">²DynFluid Lab., ENSAM, 151 boulevard de l'hôpital, 75013, Paris, France</footnote>

<footnote n="3">³LIMSI, CNRS, Université Paris-Saclay, Bât. 508, Rue John Von Neumann, Campus Universitaire, F-91405 Orsay CEDEX, France</footnote>

<footnote n="4">⁴Institute of Mechanics, BAS, 4 acad. G. Bontchev, 1113 Sofia, Bulgaria</footnote>

<footnote n="5">⁵Institut Pprime, CNRS - Université de Poitiers - ISAE-ENSMA, ENSIP, 6 Rue Marcel Doré, Bâtiment B17 - BP 633, 86022 Poitiers CEDEX, France</footnote>

<footnote>† Email address for correspondence: catherine.weisman@limsi.fr</footnote>
infinite parallel adiabatic plates. These vortices, also called outer cells, are the main object of the present study.

Rayleigh-Nyborg-Westervelt (RNW) classical approach. Consider a standing plane acoustic wave of wavelength $\lambda$ and angular frequency $\omega$ propagating in a wide channel formed by two infinite parallel plates. Assuming that the viscous penetration depth is small compared to the channel width and that they are both small compared to the wavelength, Rayleigh (1884) established the equations that govern the second-order velocity (responsible for the streaming motion), based on successive approximations. These second-order equations were obtained in the case of slow streaming motion (associated with low acoustic amplitudes) and therefore they are linear. The source terms that create the streaming motion are provided by the solution of the first-order problem (linear acoustics). In Rayleigh's solution the first-order motion is assumed to be divergence free in the viscous boundary layer. The resulting reference solution for the second-order velocity, far away from the guide wall is

$$ \bar{u}_{Rayleigh} = -\frac{3U_{ac}^2}{16c_0} \sin(2kx) \left[ 1 - 3 \left( \frac{y}{R} \right)^2 \right], \quad (1.1) $$

$$ \bar{v}_{Rayleigh} = -\frac{3U_{ac}^2}{8c_0} kR \cos(2kx) \left[ \frac{y}{R} - \left( \frac{y}{R} \right)^3 \right], \quad (1.2) $$

with $\bar{u}_{Rayleigh}$ (respectively $\bar{v}_{Rayleigh}$) the axial (respectively vertical) component of the streaming velocity i.e. the time average over an acoustic period of the velocity component, $U_{ac}$ the acoustic velocity amplitude at its antinode, $c_0$ the speed of sound, $R$ the channel half-width, $x$ the axial coordinate, $k$ the complex wavenumber and $y$ the distance from the guide axis. Rayleigh's analysis was applied to the case of a large cylindrical guide by Schuster & Matz (1940) yielding the following second-order velocity components:

$$ \bar{u}_{Rayleigh} = -\frac{3U_{ac}^2}{8c_0} \sin(2kx) \left[ 1 - 2 \left( \frac{y}{R} \right)^2 \right], \quad (1.3) $$

$$ \bar{v}_{Rayleigh} = -\frac{3U_{ac}^2}{8c_0} kR \cos(2kx) \left[ \frac{y}{R} - \left( \frac{y}{R} \right)^3 \right], \quad (1.4) $$

where here $R$ is the cylinder radius. In order to obtain these solutions, terms in $e^{-(R-y)/\delta_v}$ were neglected in the complete Rayleigh solution, with $\delta_v = \sqrt{2\nu/\omega}$ the acoustic viscous boundary layer thickness. Here $\nu$ is the gas kinematic viscosity, related to the gas dynamic viscosity $\mu$ by the classical relation $\nu = \mu/\rho_0$, where $\rho_0$ is the reference gas density. This is valid far enough from the wall. These expressions show the main characteristics of Rayleigh streaming: velocities are of second order in powers of acoustic Mach number $M = U_{ac}/c_0$, with counter-rotating vortices that span $\lambda/4$ in the axial direction. The transverse streaming velocity is much smaller and has a $\lambda/4$ spatial phase shift with the axial streaming velocity.

Later Westervelt (1953) and Nyborg (1953) improved this analytical model by proposing a complete solution to the problem considered by Rayleigh. This solution is based on a vorticity equation and there is no assumption on the divergence of the acoustic velocity. The theoretical description issued from the works of Rayleigh, Westervelt and Nyborg is usually associated with the so-called RNW streaming theory Lighthill (1978).
*Sources of streaming.* The analytical studies mentioned above (and most of the subsequent ones) are based on the same process of derivation of the streaming velocity: write fundamental laws of fluid motion, use a perturbation method with asymptotic expansions and then time average. As discussed by Lighthill (1978), this process leads to an equation of the form

$$ \mu \nabla^2 \bar{u} = -\bar{F} + \nabla \bar{p}, \qquad (1.5) $$

where $\bar{p}$ is the second-order pressure and $\bar{F}$ is a product of first-order quantities and their derivatives. The source terms composing $\bar{F}$ are based on acoustic quantities. They are also referred to as Reynolds stresses.

In all studies mentioned above, the source terms that appear in the transverse momentum equation are neglected, as discussed by Lighthill (1978). The averaged acoustic velocity source term $\bar{F}$ (Reynolds stresses) is approximated as the sum of the two largest source terms, that are associated with the axial momentum equation. The variation of these two source terms with the radius was analysed analytically in annular resonators, in the case of a pure travelling wave by Amari, Gusev & Joly (2003). An investigation of the streaming flow generation through the analysis of the contribution of each source term was conducted recently by Paridaens, Kouidri & Jebali Jerbi (2013), in the case of a thermoacoustic system. Again, only the two dominant Reynolds stress terms were considered.

One of the goals of the present study is to show that one of the terms that compose $\bar{F}$ and that is classically neglected in the analysis of Rayleigh streaming can be responsible for drastic change of acoustic streaming patterns within a given range of guide geometry. It will be shown that it can be the case even in the linear regime. A modified analytical expression for the streaming velocity is proposed, that takes into account the other two acoustic source terms associated with the transverse component of the momentum equation. This represents an improved approximation for the Rayleigh problem. This enriched solution exhibits an unexpected balance between outer Rayleigh streaming sources.

*Beyond RNW streaming.* As stated above, RNW analysis concerns streaming in wide channels outside the viscous boundary layers that results in large counter-rotating vortices usually referred to as outer vortices. These vortices are actually generated by a driving mean flow taking place inside the viscous boundary layers, associated with slender counter-rotating vortices usually referred to as inner streaming vortices. A sketch of inner and outer streaming cells, as well as transverse and axial profiles of the axial streaming velocity are shown in figure 1, in the case of a $\lambda/2$ resonator. Schlichting (1932) first described the inner streaming vortices, under the hypothesis of incompressible flow. Reviews of theoretical investigations for inner and outer streaming were proposed by Nyborg (1965) and Zaremba (1971). As shown by figure 1, the limit between inner and outer streaming cells is generally considered to be equal to $3\delta_v$ from the wall. In the present study the evolution of this limit with the dimensions of the channel will be briefly addressed.

Rott (1974) improved the description of outer streaming by taking into account thermal effects and Qi (1993) included the fluid compressibility and temperature effects to improve the description of inner streaming. Investigations of streaming in channels of arbitrary width with a mean temperature gradient were performed by Waxler (2001) and Bailliet et al. (2001) in the context of thermoacoustic applications where heat transport associated with the streaming flow is an important mechanism that may limit the efficiency of the systems. Sugimoto (2016) obtained a nonlinear
<figure><img src="image_121.png" /><figcaption>FIGURE 1. (Colour online) RNW theory. (II) sketch of inner and outer streaming cells in the case of two infinite parallel plates, showing lines (c) and (d). (I) axial acoustic velocity $u'$ profile along line (c). (III) axial streaming velocity $\bar{u}$ profile along line (d). (IV) axial streaming velocity $\bar{u}$ profile along line (c).</figcaption></figure>

thermoacoustic wave equation inside narrow channels and the associated streaming flow results can be applied to pore modelling in thermoacoustic engines.

Also in thermoacoustics the overall mass transport velocity is usually considered instead of the time-averaged velocity:

$$ \bar{u}^M = \frac{\bar{\rho} \bar{u}}{\bar{\rho}}. \qquad (1.6) $$

Here the previously considered $\bar{u}$ velocity is called Eulerian streaming velocity. Hamilton, Ilinskii & Zabolotskaya (2003a), Olson & Swift (1997) adopted a Lagrangian approach, considering the density weighted average velocity $\bar{u}^M$, also referred to as the Favre average (Favre 1965), because it is responsible for effective mass and heat transport. Hamilton et al. (2003a), Hamilton, Ilinskii & Zabolotskaya (2003b) proposed an analytical solution for the average mass transport velocity with or without temperature effects. They studied the influence of the channel width on the relative size of the inner and outer streaming vortices. Most of the present study is conducted following an Eulerian approach but the consequences of the results presented thereafter will also be briefly discussed following a Lagrangian approach.

Thermoacoustic applications which require large amplitude acoustic waves also led the authors to reconsider the hypothesis of slow streaming that is one of the bases of RNW streaming, as discussed in the extensive review proposed by Boluriaan & Morris (2003). In the last decades experimental studies (Thompson, Atchley & Maccarone
2004; Moreau, Bailliet & Valière 2008) and numerical studies (Menguy & Gilbert 2000a; Daru et al. 2013; Reyt et al. 2013) on acoustic streaming from slow to fast regimes were performed (the fast regime is associated with high amplitude acoustics). These studies exhibited streaming behaviour departing from the RNW theory: in the fast regime the dependence of the streaming velocity on the acoustic velocity is no longer parabolic and new streaming cells can emerge near the guide axis. In the slow regime, it was found that the axial streaming velocity along the axis, when normalized by Rayleigh's solution, is a linear function of the ratio $\delta_v/R$ (Daru et al. 2017a). Another goal of the present study is to quantify further this linear function, and also to show a dependency on the aspect ratio (length over width or length over radius) of the wave guide.

The hypotheses are listed in § 2, along with a description of the geometries under study. Section 3 then focuses on the rectangular case. The classical Rayleigh solution of the (first order) acoustic problem is described (§ 3.1) and the source terms for the streaming flow are deduced (§ 3.2). An analysis of the influence of each of the source terms is conducted and their effect on streaming flow patterns is highlighted (§§ 3.2.1–3.2.6). The effect of both Reynolds stress source terms, classically considered, is shown first, and then the effects of the third and fourth source terms, classically neglected, are also discussed.

The influence of the aspect ratio of the channel geometry on the streaming velocity is quantified. The results show that for specific values of the aspect ratio new streaming counter-rotating cells appear next to the guide axis. Then other characteristics of the streaming flow are examined in §§ 3.2.7–3.2.9: slip velocity, mass transport velocity, transverse component of streaming velocity. The axisymmetric case is briefly addressed in § 4.

## 2. Position of the problem

In the present approach, the discussed geometry is valid for both a cylindrical acoustic guide and acoustic flow between two infinite parallel plates. We consider a rectangular (or cylindrical) channel, of length $L$, and half-width (or radius) $R$, filled with a perfect gas of density $\rho_0$ and pressure $p_0$. The sound velocity in the ideal gas is $c_0 = \sqrt{\gamma p_0 / \rho_0}$, $\gamma$ being the specific heat ratio ($\gamma = 1.4$ for air). We use a coordinate system $(x, y)$ (or $(x, r)$) with the origin located at the centre of the channel, that is $-L/2 \le x \le L/2$, $-R \le y \le R$ (or $0 \le r \le R$).

By exploiting the symmetry of the problem, only the upper half-channel $0 \le y \le R$ is considered in the plane case. The symbolic computational software Mathematica (Wolfram 2018) is used to solve the equations throughout this study.

Throughout this work, the angular frequency is always considered to be equal to $\omega_0$, the angular frequency corresponding to the acoustic standing wave resonating at its first mode along $x$ (so-called $\lambda/2$ mode), in the case of a non-viscous gas ($\omega_0 = \pi c_0 / L$).

It is assumed that an acoustic wave propagates laminarly along the $x$ axis of the guide. There is no mean flow apart from acoustic streaming. A main hypothesis is that the wave is plane. More precisely the condition that only plane modes propagate is that the working frequency is lower than the first mode cutoff frequency. In cylindrical guides the cutoff frequency is $f_{c(1,0)} = 1.84c_0/2\pi R$ (e.g. Pierce 1989) and the condition of a plane wave is equivalent to $\lambda > 2\pi/1.84R$. In the case of a half-wavelength guide, the plane wave condition becomes $L/R > 1.7$. For rectangular channels the cutoff frequency is $c_0/2R$ or equivalently $L/R > 1$.
Another important hypothesis is that the effect of viscosity is assumed to be concentrated in the near wall region, that is, volume viscosity is neglected because the real part of the wavenumber associated with volume viscosity is negligible compared to the one associated with the viscous boundary layer. This hypothesis is valid if $1/(ReSh) \ll 1$, where $Re = \rho_0 c_0^2/(\mu\omega)$ is the acoustic Reynolds number and $Sh = \delta_v/R$ is the Shear number (Menguy & Gilbert 2000b).

We also make the assumption of an isentropic flow, in order to exclude thermal effects on streaming that are not within the scope of this work. Under these hypotheses, acoustic pressure is independent of $y$, and, as a consequence of isentropy, so is acoustic density.

Finally, several hypotheses on streaming flow are imposed. Only channels with moderate to large width values ($R/\delta_v \gtrsim 6$) are considered here in which outer streaming vortices do exist (Bailliet et al. 2001; Hamilton et al. 2003a). Another hypothesis is that $(R/L)^2 \ll 1$ which means that momentum diffusion in the axial direction is negligible with respect to momentum diffusion in the transverse direction. Moreover the present study concerns the slow streaming regime, in which the streaming flow equations are linear. This hypothesis is valid if inertial effects are negligible, which is usually related to small values of the nonlinear Reynolds number $Re_{NL} = (M/Sh)^2 \ll 1$ (Menguy & Gilbert 2000b).

### 3. Flow between two parallel plates

#### 3.1. Acoustics

Let us note $u'$ and $v'$ the axial and vertical components of the acoustic velocity, and $\rho' = \rho - \rho_0$, $p' = p - p_0$ the associated density and pressure fluctuations. Using the hypothesis $\partial^2 u'/\partial x^2 \ll \partial^2 u'/\partial y^2$, the linearized equations governing acoustics are

$$ \frac{\partial \rho'}{\partial t} + \rho_0 \left( \frac{\partial u'}{\partial x} + \frac{\partial v'}{\partial y} \right) = 0 \qquad (3.1a) $$

$$ \rho_0 \frac{\partial u'}{\partial t} + \frac{\partial p'}{\partial x} = \mu \frac{\partial^2 u'}{\partial y^2} \qquad (3.1b) $$

$$ \frac{\partial p'}{\partial y} = 0. \qquad (3.1c) $$

Using the assumption of isentropic flow, we have $p' = c_0^2 \rho'$.

The following change of variables is introduced

$$ \hat{y} = y/\delta_v, \quad \hat{x} = x/L, \quad \hat{t} = t/T, \quad \hat{R} = R/\delta_v, \quad \hat{L} = L/\delta_v, \qquad (3.2a-e) $$

$T$ being the acoustic time period so that $L/T = (1/2)c_0$. The half-channel is then described by $-1/2 \le \hat{x} \le 1/2$, $0 \le \hat{y} \le \hat{R}$ and the dependency on $\hat{x}$ is of $e^{i\pi\hat{x}}$ form. Let us note that this is only a convenient change of variables that allows us to work with non-dimensional parameters $\hat{R}$ and $\hat{L}$ instead of dimensional values of $R$ and $L$. Velocity, pressure and density are left in dimensional units.

The boundary and symmetry conditions for $u'$ and $v'$ with respect to $\hat{y}$ are

$$ u' = 0 \quad \text{at} \quad \hat{y} = \hat{R}; \quad \frac{\partial u'}{\partial \hat{y}} = 0 \quad \text{at} \quad \hat{y} = 0, \qquad (3.3a,b) $$

$$ v' = 0 \quad \text{at} \quad \hat{y} = 0 \quad \text{and} \quad \hat{y} = \hat{R}. \qquad (3.3c) $$
Following Rayleigh (1884), the axial acoustic velocity is expressed as

$$u'(\hat{x}, \hat{y}, \hat{t}) = U_{ac} \cos(\pi \hat{x}) [(-1 + e^{\hat{y}-\hat{R}} \cos(\hat{R} - \hat{y})) \cos(2\pi \hat{t}) \\ + e^{\hat{y}-\hat{R}} \sin(\hat{R} - \hat{y}) \sin(2\pi \hat{t})], \quad (3.4)$$

where $U_{ac}$ is the acoustic amplitude at the velocity antinode.

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

$$v'(\hat{x}, \hat{y}, \hat{t}) = -\frac{\pi U_{ac}}{2\hat{R}} \frac{1}{\hat{L}} \sin(\pi \hat{x}) [\hat{y}(\cos(2\pi \hat{t}) + \sin(2\pi \hat{t})) + \hat{R}e^{\hat{y}-\hat{R}}(\sin(\hat{R} - \hat{y} - 2\pi \hat{t}) \\ - \cos(\hat{R} - \hat{y} - 2\pi \hat{t})) + e^{-\hat{R}}(\hat{R} - \hat{y})(\cos(\hat{R} - 2\pi \hat{t}) - \sin(\hat{R} - 2\pi \hat{t}))]. \quad (3.5)$$

Finally the density fluctuation $\rho'$ is obtained by integrating (3.1a) with respect to time

$$\rho'(\hat{x}, \hat{t}) = \frac{\rho_0 U_{ac}}{2c_0 \hat{R}} \sin(\pi \hat{x}) [(1 - 2\hat{R}) \sin(2\pi \hat{t}) - \cos(2\pi \hat{t}) \\ + e^{-\hat{R}}(\cos(\hat{R} - 2\pi \hat{t}) + \sin(\hat{R} - 2\pi \hat{t}))]. \quad (3.6)$$

### 3.2. Streaming flow

The averaged equations are derived from the Navier–Stokes equations with each variable $f$ divided into a fluctuating, periodic, component $f'$, and a steady component $\bar{f}$ (corresponding to the streaming flow) according to

$$f = \bar{f} + f'. \quad (3.7)$$

where the overline denotes the average operator in time over one period. For any two variables $f$ and $g$

$$\overline{fg} = \overline{f}\overline{g} + \overline{f'}\overline{g'}. \quad (3.8)$$

Using the hypotheses that $\partial^2/\partial x^2 \ll \partial^2/\partial y^2$ and $Re_{NL} = (M/Sh)^2 \ll 1$, the Navier–Stokes equations are averaged in time over one period, linearized and simplified as

$$\left\{
\begin{aligned}
\rho_0 \left( \frac{\partial \bar{u}}{\partial x} + \frac{\partial \bar{v}}{\partial y} \right) &= -\frac{\partial}{\partial x} (\bar{\rho}' u') - \frac{\partial}{\partial y} (\bar{\rho}' v') \\
\mu \frac{\partial^2 \bar{u}}{\partial y^2} &= \frac{\partial \bar{p}}{\partial x} + \rho_0 \frac{\partial}{\partial x} (\overline{u' u'}) + \rho_0 \frac{\partial}{\partial y} (\overline{u' v'}) \\
\mu \frac{\partial^2 \bar{v}}{\partial y^2} &= \frac{\partial \bar{p}}{\partial y} + \rho_0 \frac{\partial}{\partial x} (\overline{u' v'}) + \rho_0 \frac{\partial}{\partial y} (\overline{v' v'}).
\end{aligned}
\right.
\quad (3.9)$$

where $\bar{u}$, $\bar{v}$ and $\bar{p}$ are the averaged velocity components and pressure.

Making the same change of variables as in § 3.1, equations (3.9) above are rewritten as

$$\left\{
\begin{aligned}
\rho_0 \left( \frac{\partial \bar{u}}{\partial \hat{x}} + \hat{L} \frac{\partial \bar{v}}{\partial \hat{y}} \right) &= -\frac{\partial}{\partial \hat{x}} (\bar{\rho}' \hat{u}') - \hat{L} \frac{\partial}{\partial \hat{y}} (\bar{\rho}' \hat{v}') \\
\rho_0 \frac{\pi}{2} c_0 \frac{\partial^2 \bar{u}}{\partial \hat{y}^2} &= \frac{\partial \bar{p}}{\partial \hat{x}} + \rho_0 \frac{\partial}{\partial \hat{x}} (\overline{\hat{u}' \hat{u}'}) + \rho_0 \hat{L} \frac{\partial}{\partial \hat{y}} (\overline{\hat{u}' \hat{v}'}) \\
\rho_0 \frac{\pi}{2} c_0 \frac{\partial^2 \bar{v}}{\partial \hat{y}^2} &= \hat{L} \frac{\partial \bar{p}}{\partial \hat{y}} + \rho_0 \frac{\partial}{\partial \hat{x}} (\overline{\hat{u}' \hat{v}'}) + \rho_0 \hat{L} \frac{\partial}{\partial \hat{y}} (\overline{\hat{v}' \hat{v}'})
\end{aligned}
\right.
\quad (3.10)$$
The boundary conditions for $\bar{u}$ and $\bar{v}$ are

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

$$ \bar{v} = 0 \quad \text{at} \quad \hat{y} = 0 \quad \text{and} \quad \hat{y} = \hat{R}, \qquad (3.11c) $$

$$ \int_{0}^{\hat{R}} \bar{u} \, d\hat{y} = 0. \qquad (3.11d) $$

The third condition in (3.11) is obtained by integrating the continuity equation in (3.10) over the half-width of the channel.

The expression of the five averaged products (Reynolds stresses) on the right-hand side of (3.10) are obtained from (3.4), (3.5) and (3.6):

$$ \overline{\rho' u'} = \frac{\rho_0 U_{ac}^2}{8c_0 \hat{R}} \sin(2\pi\hat{x}) \{1 - e^{\hat{y}-\hat{R}}(\cos(\hat{R}-\hat{y}) + (2\hat{R}-1)\sin(\hat{R}-\hat{y})) \\ + e^{\hat{y}-2\hat{R}}(\cos\hat{y} + \sin(\hat{y})) - e^{-\hat{R}}(\cos(\hat{R}) + \sin(\hat{R}))\}, \qquad (3.12) $$

$$ \overline{\rho' v'} = \frac{\pi \rho_0 U_{ac}^2}{4c_0 \hat{R}} \frac{1}{\hat{L}} \sin^2(\pi \hat{x}) \{\hat{y} + e^{\hat{y}-\hat{R}}(-\hat{R} \cos(\hat{R}-\hat{y}) + \sin(\hat{R}-\hat{y})(1-\hat{R})) \\ + e^{\hat{y}-2\hat{R}} \sin \hat{y} + e^{-\hat{R}}((\hat{R}-\hat{y}) \cos \hat{R} + (\hat{R}-\hat{y}-1) \sin \hat{R})\}, \qquad (3.13) $$

$$ \overline{u'u'} = \frac{1}{2} U_{ac}^2 \cos^2(\pi \hat{x}) (1 + e^{2(\hat{y}-\hat{R})} - 2e^{\hat{y}-\hat{R}} \cos(\hat{R}-\hat{y})), \qquad (3.14) $$

$$ \overline{u'v'} = \frac{\pi U_{ac}^2}{8 \hat{R} \hat{L}} \sin(2\pi \hat{x}) \{\hat{y} + e^{\hat{y}-\hat{R}}(-\hat{R} + \hat{y}) \cos(\hat{R}-\hat{y}) + (\hat{R}-\hat{y}) \sin(\hat{R}-\hat{y}) \\ + \hat{R} e^{2(\hat{y}-\hat{R})} + e^{\hat{y}-2\hat{R}}((\hat{y}-\hat{R}) \cos \hat{y} + (\hat{R}-\hat{y}) \sin \hat{y}) \\ + e^{-\hat{R}}((\hat{R}-\hat{y}) \cos \hat{R} - (\hat{R}-\hat{y}) \sin \hat{R})\}, \qquad (3.15) $$

$$ \overline{v'v'} = \frac{\pi U_{ac}^2}{\hat{L}^2} \frac{1}{4\hat{R}^2} \sin^2(\pi \hat{x}) \{\hat{y}^2 - 2e^{\hat{y}-\hat{R}}\hat{y}\hat{R} \cos(\hat{R}-\hat{y}) + e^{2(\hat{y}-\hat{R})}\hat{R}^2 \\ + 2\hat{R}e^{\hat{y}-2\hat{R}}(\hat{y}-\hat{R})\cos\hat{y} + 2e^{-\hat{R}}\hat{y}(\hat{y}+\hat{R})\cos\hat{R} + e^{-2\hat{R}}(\hat{R}-\hat{y})^2\}. \qquad (3.16) $$

Under the hypothesis $\hat{R} \gtrsim 6$, terms of the same magnitude as $e^{-\hat{R}}$ or smaller can be neglected in the previous expressions of the averaged products. Nevertheless, all symbolic calculations were performed using the full expressions.

Following Westervelt (1953), we take the curl of the momentum equations, which eliminates pressure, and use the assumption $\partial \bar{v}/L \partial \hat{x} \ll \partial \bar{u}/\delta_v \partial \hat{y}$, to obtain

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

Classical studies (e.g. Rayleigh 1884, Westervelt 1953, Bailliet et al. 2001, Hamilton et al. 2003a) always consider an added approximation in which the last two source terms on the right-hand side of (3.17) are negligible. Neglecting these two terms was justified by an analysis of the orders of magnitude, based on the variations in the
viscous boundary layer. The usual assumptions are $\partial/\partial x \propto 1/\lambda \propto 1/L$ and $\partial/\partial y \propto 1/\delta_v$ so that $\partial/\partial \hat{x} \propto 1$ and $\partial/\partial \hat{y} \propto 1$. Therefore equations (3.14), (3.15), (3.16) yield

$$ \hat{L} \frac{\partial^2(\overline{u'u'})}{\partial \hat{x} \partial \hat{y}} \propto \hat{L} U_{ac}^2, \quad \hat{L}^2 \frac{\partial^2(\overline{u'u'})}{\partial \hat{y}^2} \propto \hat{L} U_{ac}^2, \quad \frac{\partial^2(\overline{u'u'})}{\partial \hat{x}^2} \propto \frac{U_{ac}^2}{\hat{L}}, \quad \hat{L} \frac{\partial^2(\overline{v'v'})}{\partial \hat{x} \partial \hat{y}} \propto \frac{U_{ac}^2}{\hat{L}}. \quad (3.18a-d) $$

It follows that the third and fourth terms are respectively $1/\hat{L}^2$ smaller than the first two terms.

In his analysis, Lighthill (1978) states that making a decision to keep or neglect a term should be based on an order of magnitude analysis and not a perturbation expansion alone. In the present work it will be shown that even an order of magnitude analysis on the acoustic source terms is not sufficient to neglect a term. This is because the different terms in (3.17) do not have the same $\hat{y}$ dependence; some are proportional to $\hat{y}$, others to $\hat{y}^2$ and others do not have any $\hat{y}$ dependence so that the successive integrations over $\hat{y}$ needed to obtain $\bar{u}$ from (3.17) together with the boundary conditions will change the balance between the effects of these different source terms. Even though they do not have the same order of magnitude, the source terms can have effects of the same order on the streaming velocity. Therefore, they will all be considered. (The change in balance between the source terms can take place independently of the fact that $\partial \bar{v}/L \partial \hat{x} \ll \partial \bar{u}/\delta_v \partial \hat{y}$, which is given by the fundamental equations for streaming (3.9). Thus the process of keeping all the terms on the right-hand side of (3.17) is consistent with neglecting the $\bar{v}$ term on its left-hand side.)

In order to analyse in detail the effect of each of the source terms, we will consider separately the four problems linked to each source term. We will normalize the solutions with Rayleigh's solution on the axis $\bar{u}_R = -(3/16)U_{ac}^2/c_0 \sin(2\pi\hat{x})$. The solution of problem i, associated with the ith source term, will be denoted by $\bar{u}_i$, $i = 1, 2, 3, 4$. In the following, for the sake of simplicity, the expressions for each solution $\bar{u}_i$ will be presented in an approximated form, in which terms of the same magnitude as $e^{-\hat{R}}$ or smaller will be neglected with respect to terms of order $1/\hat{R}^n$, for $n \in \mathbb{N}$. However the symbolic computations are still performed using the full expressions.

### 3.2.1. Problem 1: first source term

When reducing equation (3.17) by taking into account the first source term only, the equation to be solved is

$$ \frac{\pi}{2} c_0 \frac{\partial^3 \bar{u}_1}{\partial \hat{y}^3} = \frac{\partial^2(\overline{u'u'})}{\partial \hat{x} \partial \hat{y}}, \quad (3.19) $$

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

$$ \bar{u}_1(\hat{x}, \hat{y}) = -\frac{U_{ac}^2}{16c_0} \sin(2\pi\hat{x}) \left\{ 15 \frac{\hat{y}^2}{\hat{R}^3} - \frac{15}{\hat{R}} - 6 \frac{\hat{y}^2}{\hat{R}^2} + 2(1 + 2e^{2(\hat{y}-\hat{R})} + 8e^{\hat{y}-\hat{R}} \sin(\hat{R}-\hat{y})) \right\}. \quad (3.20) $$
Along the axis, by normalizing with Rayleigh's solution and neglecting smaller terms, we obtain at $O(1/\hat{R})$:

$$ \frac{\bar{u}_1(\hat{x}, 0)}{\bar{u}_R} = \frac{2}{3} - \frac{5}{\hat{R}}. \qquad (3.21) $$

The second-order pressure gradient can also be obtained, and as expected is independent of $\hat{y}$. It is given by

$$ \frac{d\bar{p}_1}{d\hat{x}} = \frac{\rho_0 \pi U_{ac}^2}{16} \sin(2\pi\hat{x}) \left(8 + \frac{6}{\hat{R}^2} - \frac{15}{\hat{R}^3}\right). \qquad (3.22) $$

### 3.2.2. Problem 2: second source term

The equation to be solved for the acoustic streaming issued from the second source term in equation (3.17) is

$$ \frac{\pi}{2} c_0 \frac{\partial^3 \bar{u}_2}{\partial \hat{y}^3} = \hat{L} \frac{\partial^2 (\bar{u}'\hat{v}')}{\partial \hat{y}^2}, \qquad (3.23) $$

with boundary conditions (3.11a,b)–(3.11d) for $\bar{u}_2$. Integrating (3.23) gives the following simplified expression for $\bar{u}_2$

$$ \bar{u}_2(\hat{x}, \hat{y}) = -\frac{U_{ac}^2}{32c_0} \sin(2\pi\hat{x}) \left\{ \frac{1}{\hat{R}} + 3\frac{\hat{y}^2}{\hat{R}^3}(1-2\hat{R}) + 2(1-2e^{2(\hat{y}-\hat{R})}) \right. \\ \left. + \frac{4}{\hat{R}}e^{\hat{y}-\hat{R}}[-(-1+2\hat{y})\cos(\hat{R}-\hat{y}) + (1-2\hat{R})\sin(\hat{R}-\hat{y})] \right\}. \qquad (3.24) $$

Along the axis, we obtain at the order $O(1/\hat{R})$:

$$ \frac{\bar{u}_2(\hat{x}, 0)}{\bar{u}_R} = \frac{1}{3} + \frac{1}{6\hat{R}}, \qquad (3.25) $$

and the corresponding pressure gradient is

$$ \frac{d\bar{p}_2}{d\hat{x}} = -\frac{\rho_0 \pi U_{ac}^2}{32} \sin(2\pi\hat{x}) \left(\frac{4}{\hat{R}} - \frac{6}{\hat{R}^2} + \frac{3}{\hat{R}^3}\right). \qquad (3.26) $$

### 3.2.3. Analysis of the effect of the first two source terms

It is interesting to analyse the right-hand side of the axial component of the momentum equation in (3.10). For the first problem, the right-hand side is $\text{RHS}_1 = \partial\bar{p}_1/\partial\hat{x} + \rho_0(\partial/\partial\hat{x})(\bar{u}'\bar{u}')$. For the second problem, the right hand side is $\text{RHS}_2 = \partial\bar{p}_2/\partial\hat{x} + \rho_0\hat{L}(\partial/\partial\hat{y})(\bar{u}'\bar{v}')$.

The right-hand sides are therefore composed of two terms, a pressure gradient term and a Reynolds stress term. Figure 2 shows, for $\hat{R} = 20$, at $\hat{x} = -1/4$, the transverse profile of both terms composing $\text{RHS}_1$ figures 2(a) and $\text{RHS}_2$ 2(b), normalized by $\rho_0 U_{ac}^2$. Comparing figure 2(a,b) shows that the terms in $\text{RHS}_2$ are much smaller than in $\text{RHS}_1$, and also that the pressure gradient and Reynolds stress terms almost compensate each other in the core of the guide, for both problems 1 and 2. The pressure gradients of problems 1 and 2 have opposite signs, and also $d\bar{p}_2/d\hat{x}$ is much
<figure><img src="image_122.png" /></figure>

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

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

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

FIGURE 2. (Colour online) Variation with $\hat{y}$, for $\hat{x} = -1/4$ and $\hat{R} = 20$ of (a) Reynolds stresses term (red solid line) and normalized pressure gradient term (red dashed line) in RHS$_1$, (b) Reynolds stresses term (blue solid line) and normalized pressure gradient term (blue dashed line) in RHS$_2$, (c) total normalized RHS$_1$ (red solid line) and RHS$_2$ (blue dashed line), (d) normalized axial streaming velocity $\bar{u}_1/\bar{u}_R$ (red solid line) and $\bar{u}_2/\bar{u}_R$ (blue dashed line).

smaller than $d\bar{p}_1/d\hat{x}$, with $d\bar{p}_2/d\hat{x} \to 0$ when $\hat{R}$ becomes large (refer to equations (3.22) and (3.26)).

Figure 2(c) shows the transverse profile of the complete right-hand sides RHS$_1$ and RHS$_2$, normalized by $\rho_0 U_{ac}^2$, again for $\hat{R}=20$, at $\hat{x}=-1/4$. This figure shows that the two complete right-hand sides are of the same order, with significant values only near the wall. This stresses the fact that streaming 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. As stated by Menguy & Gilbert (2000a) (p. 253): pressure is created directly by the acoustic wave, and not by the mean flow. Outer streaming flow can thus be considered as a driven flow, similar to a Couette flow (see the driven cavity analogy developed in Daru *et al.* (2017b)).

The resulting normalized axial streaming velocity transverse profiles are plotted on figure 2(d). The following observation can be made: the first problem is responsible for the inner streaming (Schlichting streaming), while the second problem only contributes to the outer streaming (Rayleigh streaming), as also noticed in Paridaens *et al.* (2013). This is related to the fact that RHS$_2$ vanishes for $\hat{y} = \hat{R}$, while RHS$_1$ does not, as can be seen in figure 2(c).

Let us now consider the superposition of problems 1 and 2 considered above, which is usually assumed in the literature to give the total axial streaming velocity $\bar{u}_{1+2} = \bar{u}_1 + \bar{u}_2$. Figure 3 shows the normalized axial streaming velocity transverse profiles $\bar{u}_1$, $\bar{u}_2$ and $\bar{u}_{1+2}$. This figure also illustrates that the first problem is responsible for the inner streaming, while the second problem only contributes to the outer streaming: only $\bar{u}_1$ takes a positive value near the wall, while $\bar{u}_2$ remains negative.
<figure><img src="image_126.png" /><figcaption>FIGURE 3. (Colour online) Axial streaming velocities $\bar{u}_1(-1/4, \hat{y})/\bar{u}_R$ (red solid line), $\bar{u}_2(-1/4, \hat{y})/\bar{u}_R$ (blue dashed line) and $\bar{u}_{1+2}(-1/4, \hat{y})/\bar{u}_R$ (black dot-dashed line). $\hat{R}=50$.</figcaption></figure>

The expression for $\bar{u}_{1+2}$ obtained on the axis is approximated as:

$$ \frac{\bar{u}_{1+2}(\hat{x}, 0)}{\bar{u}_R} = 1 - \frac{29}{6} \frac{1}{\hat{R}}. \qquad (3.27) $$

For large guides ($\hat{R} \gg 1$, but still $(\hat{R}/\hat{L})^2 \ll 1$), we thus recover Rayleigh's solution. As noticed by Lighthill (1978) and as shown by equations (3.21) and (3.25), the contributions of the first and second source terms in Rayleigh's solution are respectively 2/3 and 1/3. Equation (3.27) also shows a linear variation with $1/\hat{R}$, which is not negligible for moderate values of $\hat{R}$. The first source term is the main contributor to this dependency: -5 in (3.21) versus 1/6 in (3.25).

The total axial pressure gradient is

$$ \frac{d\bar{p}_{1+2}}{d\hat{x}} = \frac{\rho_0 \pi U_{ac}^2}{2} \sin(2\pi\hat{x}) \left( 1 - \frac{1}{8\hat{R}} + \frac{3}{16\hat{R}^2} - \frac{27}{32\hat{R}^3} \right). \qquad (3.28) $$

It follows that, for large values of $\hat{R}$

$$ \frac{d\bar{p}_{1+2}}{d\hat{x}} \approx \frac{d\bar{p}_1}{d\hat{x}} = \frac{\rho_0 \pi U_{ac}^2}{2} \sin(2\pi\hat{x}). \qquad (3.29) $$

Integrating equation (3.29) with respect to $\hat{x}$ yields

$$ \bar{p}_{1+2} \approx \bar{p}_1 = p_0 - \frac{\rho_0 U_{ac}^2}{4} \cos(2\pi\hat{x}). \qquad (3.30) $$

When scaling $\bar{p}_{1+2} - p_0$ by $(\gamma/4)p_0(U_{ac}/c_0)^2$, we recover the theoretical result for the dimensionless hydrodynamic streaming pressure known in the literature $P_{2s} = -\cos(2\pi\hat{x})$ (e.g. Menguy & Gilbert 2000a).