# Phenomenology of laminar acoustic streaming jets

Bjarne Vincent, Daniel Henry, Abhishek Kumar, Valéry Botton, Alban

Pothérat, Sophie Miralles

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

► To cite this version:

Bjarne Vincent, Daniel Henry, Abhishek Kumar, Valéry Botton, Alban Pothérat, et al.. Phenomenology of laminar acoustic streaming jets. 2025. hal-04674347v2

HAL Id: hal-04674347

https://hal.science/hal-04674347v2

Preprint submitted on 18 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
Phenomenology of laminar acoustic streaming jets

Bjarne Vincent,<sup>1,2</sup> Daniel Henry,<sup>1</sup> Abhishek Kumar,<sup>2</sup> Valéry Botton,<sup>1,*</sup> Alban Pothérat,<sup>2</sup> and Sophie Miralles<sup>1, †</sup>

<sup>1</sup>INSA Lyon, CNRS, Ecole Centrale de Lyon, Université Claude Bernard Lyon 1,
Laboratoire de Mécanique des Fluides et d'Acoustique, UMR5509, 69621, Villeurbanne France

<sup>2</sup>Fluid and Complex Systems Research Centre, Coventry University, Coventry CV15FB, United Kingdom

This work identifies the physical mechanisms at play in the different flow regions along an Eckart acoustic streaming jet by means of numerical simulation based on a novel modeling of the driving acoustic force including attenuation effects. The flow is forced by an axisymmetric beam of progressive sound waves attenuating over a significant part of a closed cylindrical vessel where the jet is confined. We focus on the steady, axisymmetric and laminar regime. The jet typically displays a strong acceleration close to the source before reaching a peak velocity. At further distances from the transducer, the on-axis jet velocity smoothly decays before reaching the opposite wall. For each of these flow regions along the jet, we derive scaling laws for the on-axis velocity with the magnitude of the acoustic force and the diffraction of the driving acoustic beam. These laws highlight the different flow regimes along the jet and establish a clear picture of its spatial structure, able to inform the design of experimental or industrial setups involving Eckart streaming jets.

# I. INTRODUCTION

This study aims to characterize the jet flow produced by an axisymmetric beam of progressive, attenuated acoustic waves. We focus here on determining accurate scaling laws for the velocity along the jet in a fluid domain significantly larger than the attenuation length of the sound waves, a geometry that is representative of a number of applications of acoustically-driven flows.

Driving flows by means of acoustic waves, or acoustic streaming, relies on a nonlinear phenomenon where momentum is generated by the attenuation of sounds waves [1]. In this work, we focus on the Eckart-type acoustic streaming, that is a flow forced by traveling waves being attenuated in the bulk of a fluid [2, 3] (from now, any mention to acoustic streaming shall refer to Eckart streaming, unless specified otherwise). Because of its contactless nature, acoustic streaming has aroused considerable interest in applications in which the flow actuation by mechanical means is rendered difficult by tight geometrical constraints. This is commonly the case of microfluidics, in which the recent developments make an extensive use of sound waves to drive fluid flows and to possibly promote mixing in microcavities [4–7], in droplets [8, 9] and around oscillating bubbles [10–12]. Acoustic waves are also often used to trap and manipulate small objects by means of “acoustic tweezers” [13–15].

At a larger scale, acoustic waves constitute an appropriate way of remotely interacting with delicate fluids, for which the contamination by intrusive bodies must be avoided. This is for instance the case in metallurgy, where the access to the melts, in solidification processes, is further hindered by the important temperatures encountered and by tight thermal insulation constraints. In this context, electromagnetic [16, 17] and mechanical [18] stirring of the melt have proven to be beneficial for the quality of the obtained solid phase. Eckart streaming can thus be considered as an alternative to the more intrusive actuations. It may be used to promote mixing in the liquid phase [19, 20] or to directly act on the solid-liquid interface itself [21–23]. In such industrial applications, large-scale Eckart streaming flows (i.e., flows for which the characteristic length scales are significantly larger than the size of the source) are commonly obtained with high power and low frequency ultrasound waves (typically hundreds of watts and kilohertz). Due to the high acoustic wave amplitudes involved, a cloud of cavitation bubbles is produced close to the transducer and acts as the main sound attenuation mechanism driving the flow [24–26]. While this approach leads to strong jets, its main drawback is the loss of coherence of the sound beam over distances comparable to the size of the acoustic source. It thus restricts the controlled actuation of the fluid to regions close to the actuator.

We are instead interested in higher frequencies and lower powers, for which the beam remains coherent over larger distances. This makes possible to act on regions of fluid located far from the acoustic transducer through the generation of elongated jets [27–31]. These jets can reach distances even greater than the acoustic beam attenuation length $\alpha^{-1}$, with $\alpha$ the acoustic pressure wave attenuation coefficient depending on the wave frequency and the material properties of the fluid. An important parameter for such flows is thus the quantity $N = \alpha L_s$, which compares the distance $L_s$

<footnote>* valery.botton@insa-lyon.fr</footnote>

<footnote>† sophie.miralles@insa-lyon.fr</footnote>
over which the Eckart streaming flow is observed to $\alpha^{-1}$. To ensure that most of the available energy radiated by the acoustic source is used to create flow kinetic energy, $N > 1$ shall be aimed for.

Still, as most of the experimental studies focused on streaming flows generated on small distances compared to the attenuation length ($N < 1$), the current knowledge of large-scale Eckart streaming remains too scarce to efficiently implement it in an engineering context where large actuation distances are sought. From the water experiments of Kamakura et al. [27] ($N = 0.17$) and Moudjed et al. [31] ($N = 0.05$), the long streaming jet strongly accelerates over distances greater than the transducer diameter but smaller than $\alpha^{-1}$. In this region, the on-axis velocity scales with the square root of the acoustic power, as a result of a balance between inertia and the acoustic force [31]. At further distances from the source, Mitome [29] ($N = 0.4$ and $N = 0.45$) and Dentry et al. [30] ($0.63 \le N \le 1400$) observed that the jet velocity reaches a maximum before decreasing further away from the source; however, the authors did not explain these jet velocity features. A global picture of the evolution of the streaming velocity with the distance from the source is thus still lacking. Aside from the cited experimental works, a coherent description of a jet obtained numerically, with justified assumptions for the modeling of the acoustic force, is scarce [14].

These Eckart streaming experiments call for numerical simulations to understand the scaling of streaming jets with the acoustic forcing. Direct Numerical Simulations (DNS) of the compressible Navier-Stokes equations, i.e., solving the wave propagation problem and waiting for a mean flow to build up, may appear as a natural approach to compute streaming jet flows [32, 33]. However, as they require to solve acoustic time and length scales that are usually several orders of magnitudes smaller than those of the flow [31], DNS are only tractable for streaming jet lengths that do not exceed a few wavelengths (e.g., up to 100 wavelengths in Ref. [33]). DNS are thus unsuitable for the long-range Eckart streaming we consider. Another approach relies on invoking the time scale discrepancy between the flow and acoustic problems to derive and solve a reduced set of equations describing the streaming flow [31, 34, 35]. With a method similar to RANS (Reynolds-Averaged Navier-Stokes) modeling in turbulence, Moudjed et al. [31] showed that the streaming flow can be described by the incompressible Navier-Stokes together with an acoustic body force driving the flow. The acoustic problem then decouples from the streaming one, thus greatly alleviating the range limit encountered with DNS. Nevertheless, this approach requires the acoustic field to be accurately modeled. In particular for $N > 1$, alteration of the acoustic field structure caused by sound attenuation is expected to fundamentally change the structure of the jet. In the present work, we adapt the usual equations governing the unattenuated sound field radiated by a baffled source [36, 37] to include attenuation effects while keeping a linear acoustic propagation workframe.

We present numerical simulations of Eckart streaming jets flowing in a closed cylindrical cavity significantly longer than the acoustic attenuation length, i.e., $N > 1$. This allows us to explore flow regions in which the acoustic force is either driving the flow (at distances to the transducer that are smaller than $\alpha^{-1}$) or nearly ineffective (at distances greater than $\alpha^{-1}$, where the acoustic beam is almost entirely attenuated). In the present study, we derive scaling laws for the velocity along the entire jet as a function of the forcing magnitude and the distance from the source. In particular, our analysis characterizes the velocity maximum and the downstream region of decaying velocity where the acoustic force becomes negligible.

This paper is organized as follows: first, we outline our modeling strategy in Sec. II. In particular, we define the considered problem, and provide the governing equations along with our numerical methodology. We then stress the central role of sound attenuation in Sec. III. Section IV is devoted to the analysis of a long-range streaming flow. Finally, concluding remarks are given in Sec. V.

## II. MODELING OF ACOUSTIC STREAMING JETS

### A. Definition of the streaming problem

We consider a steady and axisymmetric Eckart streaming flow of dimensional velocity field $u_s$ in a closed cylindrical cavity of diameter $D_c$ and length $L_c$ (Fig. 1). The fluid is Newtonian of density $\rho$ and kinematic viscosity $\nu$. A plane circular transducer of diameter $D_s$ is placed at $x = 0$ and radiates a beam-shaped sound field centered on the cylinder axis. The sound waves, of power $P_{ac}$ and of wavelength $\lambda$, propagate through the medium at a phase speed $c$. In addition to the decay caused by diffraction, the amplitude of the acoustic pressure waves drops exponentially along the axis due to viscous and thermal effects occurring within the bulk of the fluid. This decay rate is given by the sound attenuation coefficient [31, 36]:

$$ \alpha = \frac{2\pi^2\nu}{c\lambda^2} \left( \frac{4}{3} + \frac{\eta}{\mu} + \frac{\gamma - 1}{Pr} \right), \quad (1) $$

where $\eta$ and $\mu = \rho\nu$ are the bulk and dynamic viscosities, respectively, $\gamma$ is the adiabatic index and $Pr = \mu c_p / \kappa$ is the Prandtl number with $c_p$ the specific heat capacity at constant pressure and $\kappa$ the thermal conductivity. As formerly
<figure><img src="image_3.png" /><figcaption>FIG. 1. Sketch of the computational domain in the $(x_d, r_d)$ plane, where $x_d$ and $r_d$ refer to the dimensional axial and radial coordinates, respectively. The domain is a closed cylindrical cavity filled with a Newtonian fluid of density $\rho$ and kinematic viscosity $\nu$. It is fitted with a plane circular transducer (gray) at $x = 0$ that radiates an axisymmetric beam-shaped acoustic field (green). The sound waves of wavelength $\lambda$ are emitted with an initial power $P_{ac}$ and travel at the phase speed $c$. The pressure amplitude decays at a rate $\alpha$ along the wave propagation path. The attenuation of the sound beam creates a body force that drives a flow of velocity field $\mathbf{u}_s$. The symmetry axis of the cylindrical container and the flow field is the $r = 0$ line, and all boundaries are impermeable walls that completely absorb the sound waves.</figcaption></figure>

derived by Lighthill [1], the theoretical axial momentum flux of an Eckart streaming jet flowing in a semi-infinite domain increases with the distance $x_d$ from the transducer, and asymptotically reaches its maximum $M_{th,max} = P_{ac}/c$ at infinite $x_d$. In the case where the effective beam length matches $L_c$, only a fraction $(1 - e^{-2\alpha L_c})$ of $M_{th,max}$ is actually injected into the fluid volume. In the present work, we focus on a $L_c > (1/\alpha)$ case to maximize the jet momentum flux. Using $(\alpha, \lambda, M_{th,max}, \rho, \nu, D_s, D_c, L_c, \mathbf{u}_s)$ as dimensional design parameters, and $D_s, \rho D_s^3, D_s^2/\nu$ as length, mass and time scales, respectively, we describe the acoustic propagation and hydrodynamics problems using six dimensionless groups [31]:

$$N = \alpha L_c, \quad S = 1.22 \frac{\lambda}{D_s}, \quad M = \frac{P_{ac}}{\rho c \nu^2}, \quad L = \frac{L_c}{D_s}, \quad D = \frac{D_c}{D_s}, \quad \mathbf{u} = \frac{\mathbf{u}_s D_s}{\nu}. \qquad (2)$$

In Eq. (2), $S$ compares the wavelength to $D_s$. As the distance from the transducer increases, the beam radius increases due to diffraction. The acoustic field thus takes the shape of a cone of half-angle $\theta = \sin^{-1} S$ (Fig. 1). In other words, the parameter $S$ expresses how diffraction acts on the shape of the acoustic force field. The Eckart streaming problem is described by the dimensionless incompressible continuity and Navier-Stokes equations:

$$\nabla \cdot \mathbf{u} = 0, \qquad (3a)$$

$$\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla) \mathbf{u} = -\nabla p + \nabla^2 \mathbf{u} + Gr_{ac} \tilde{\mathbf{I}}. \qquad (3b)$$

The last term in Eq. (3b) is the acoustic force driving the flow and is valid far from solid boundaries and for any acoustic wave shape [9]. The spatial structure of the forcing is set by the normalized acoustic intensity $\tilde{\mathbf{I}}$ dictated by the acoustic problem (see Sec. II B). $\tilde{\mathbf{I}}$ is normalized by the maximum unattenuated plane wave intensity $I_{max}$ on the transducer axis (Eq. (A14)), giving $\|\tilde{\mathbf{I}}\| \le 1$ everywhere in the fluid domain. The magnitude of the forcing, relative to the magnitude of the viscous forces, is then given by the acoustic Grashof number $Gr_{ac}$ [23]:

$$Gr_{ac} = \frac{32\alpha P_{ac} D_s}{\pi \rho c \nu^2} = \frac{32 N M}{\pi L} = \frac{16 M}{\pi L_{\alpha}}, \qquad (4)$$

where $L_{\alpha} = L/(2N) = 1/(2\alpha D_s)$ is the dimensionless attenuation length defining the exponential decay of the force along the acoustic axis. A small value of $L_{\alpha}$ corresponds to significant sound attenuation, while $L_{\alpha} \to \infty$ denotes a vanishingly small attenuation.
<table><caption>TABLE I. Values of the dimensionless parameters considered in the present study, as defined in Eqs. (2) and (4). The values of the parameters *S* and *N/L* correspond to the water experiments of Moudjed *et al.* [31], for which *f* = 2 × 10⁶ Hz, *α* ≈ 0.1 m⁻¹, *D*<sub>s</sub> = 2.85 × 10⁻² m and *c* = 1480 m.s⁻¹.</caption><thead><tr><th>S</th><th>N</th><th>L</th><th>D</th><th>L<sub>α</sub> = L/(2N)</th><th>Gr<sub>ac</sub></th></tr></thead><tbody><tr><td>0.03</td><td>4.5</td><td>1500.7</td><td>548.3</td><td>166.7</td><td>Varied between 10<sup>3</sup> and 5 × 10<sup>4</sup></td></tr></tbody></table>

The cavity walls are modeled as impermeable boundaries, imposing **u** = 0 at *x* ∈ {0, *L*} and *r* = *D*/2. Besides, the flow is axisymmetric and swirl-free, so that for all *x* and *r*:

$$u_{\theta} = \frac{\partial u_x}{\partial \theta} = \frac{\partial u_r}{\partial \theta} = \frac{\partial p}{\partial \theta} = 0,$$

where *θ* is the azimuthal coordinate, and *u*<sub>*x*</sub>, *u*<sub>*r*</sub> and *u*<sub>*θ*</sub> refer to the axial, radial and azimuthal velocity components, respectively. Axisymmetry of both **u** and *p* further imposes

$$u_r = \frac{\partial u_x}{\partial r} = \frac{\partial p}{\partial r} = 0$$

at *r* = 0.

The values of the dimensionless parameters (2) are listed in Table I. We choose the same diffraction parameter *S* = 0.03 and attenuation coefficient *N/L* = 3 × 10⁻³ as in the water experiments of Moudjed *et al.* [31]. Similar parameters values are found in the experiments of Refs. [27, 29, 38], in which *S* ∈ [0.014, 0.036] and *N/L* ∈ [5.9 × 10⁻³, 18.8 × 10⁻³]. The acoustic forcing structure in our simulations is thus representative of those commonly found in streaming jets experiments.

An important aspect of our work is the moderately large *N*. In contrast to Refs. [27, 29, 31, 38], we choose *N* > 1. The beam is thus entirely attenuated before reaching the end of the fluid domain, yielding forcing-free regions at far distances from the source. Typically, at *x* = *L*, the remaining acoustic power is about 0.01 % of the power radiated by the transducer. To minimize the effects of the lateral wall on the jet, we choose *D* ≈ *L*/3.

Finally, we take 10³ ≤ Gr<sub>ac</sub> ≤ 5 × 10⁴, which fits in the ranges explored by Kamakura *et al.* [27], Mitome [29], Moudjed *et al.* [31], Frenkel *et al.* [38]. Some of these authors also investigated larger Gr<sub>ac</sub>. We limited ourselves to Gr<sub>ac</sub> ≤ 5 × 10⁴, as computing steady flows at greater Gr<sub>ac</sub> in such a large cavity becomes prohibitively expensive.

## B. Equations governing the spatial distribution of the attenuated acoustic force field

We shall now explicate the acoustic intensity field $\tilde{I}$, which depends only on the dimensionless group *S* and *N/L*. As *N* > 1, attenuation is expected to significantly affect $\tilde{I}$, so that the usual equations [36, 37] for the unattenuated sound field radiated by a plane circular transducer are no longer valid and must be corrected.

Although the streaming flow is sufficiently slow to be seen as incompressible, acoustic waves are a compressible process [39], and shall thus be modeled accordingly. We further assume that:

(i) Sound waves are irrotational, i.e., $\nabla \times u_{ac} = 0$. This is a common assumption in acoustics [36, 37, 40]; it holds in the bulk of the fluid but breaks down in a thin layer near solid boundaries, including at the source. There, a viscous layer of dimensionless thickness $\delta = \sqrt{\nu\lambda/(\pi c D_s^2)}$ forces the acoustic velocity fluctuation $u_{ac}$ to match the boundary velocity. This is typically the case near the vibrating membrane of the transducer radiating the acoustic waves. For the experiment in Ref. [31], on which we base our forcing field structure, $\delta \approx 1.3 \times 10^{-5}$. As $\delta$ is several orders of magnitudes smaller than the cavity and flow length scales, the supplementary attenuation and possible streaming flow within that layer can be safely neglected [41].

(ii) We assume that $L_α$ is smaller than the acoustic shock formation distance: the sound waves are thus linear and the superposition principle applies. This modelling choice removes the need for defining a frequency [31] and acoustic nonlinearity parameters [40], thus making our scaling study as generic as possible.

(iii) The waves are time-harmonic monochromatic signals with an angular frequency $\omega$.

(iv) The length and time scales of the acoustic propagation problem are significantly smaller than those of the streaming flow, so that the flow does not affect sound propagation [42]. These scale discrepancies are indeed verified in long-range laminar Eckart streaming [27, 29, 31].

(v) Sound attenuation is due to the viscosity of the medium only. Asymptotic analysis [40] shows that thermally-induced attenuation is proportional to ($\gamma - 1$). For common liquids such as water, $\gamma \approx 1$, hence thermally-induced attenuation can be dropped [9] (this approximation however does not hold for liquid metals). We also
discard any relaxation-based attenuation, e.g., through the excitation of the vibrational and rotational modes of the fluid molecules [36].

(vi) The sound attenuation length is large compared to the wavelength (i.e., $NS/(1.22L) \ll 1$).

We model the transducer as a disk vibrating at a uniform velocity $u_0e^{j\omega t}$ (with $j^2 = -1$). That disk is modeled as a theoretically infinite number of acoustic point sources, each emitting attenuated spherical waves contributing to the acoustic pressure $\tilde{p}_{ac}$, velocity $\tilde{u}_{ac}$ and $\tilde{I}$ at $(x, r)$ through (see appendix A for the complete derivation):

$$ \tilde{I} = \Re\{\tilde{p}_{ac} \tilde{u}_{ac}^*\}, \qquad (5a) $$

$$ \tilde{p}_{ac} = j \frac{k^2 S}{8 \times 1.22 \pi^2} \int_0^{1/2} \int_0^{2\pi} \frac{e^{-jk\|d\|}}{\|d\|} r' dr' d\theta', \qquad (5b) $$

$$ \tilde{u}_{ac} = \frac{1}{4\pi} \int_{0}^{1/2} \int_{0}^{2\pi} (1 + jk \|d\|) \frac{e^{-jk \|d\|}}{\|d\|^2} \frac{d}{\|d\|} r' dr' d\theta', \qquad (5c) $$

$$ k = \frac{2.44\pi}{S} - j \frac{N}{L}, \qquad (5d) $$

where * refers to the complex conjugate and $\Re$ to the real part. $d = x e_x + [r - r' \cos(\theta')] e_r - r' \sin(\theta') e_\theta$ is the distance between a point source located at the polar coordinates $(r', \theta')$ from the transducer center and a point at $(x, r)$. In Eqs. (5), $\tilde{p}_{ac}$, $\tilde{u}_{ac}$ and $\tilde{I}$ are normalized by their maximum unattenuated magnitudes, which are $2\rho c u_0$, $2u_0$ and $2\rho c u_0^2$, respectively. The effects of attenuation are accounted for through the complex wavenumber $k$ (Eq. (5d)), whose imaginary part $N/L = 1/(2L_\alpha) = \alpha D_s$ is the dimensionless attenuation coefficient. Finally, since the geometry is axisymmetric, $\tilde{u}_{ac} \cdot e_\theta$ vanishes when integrating Eq. (5c). Therefore, the acoustic fields are indeed axisymmetric and free of any azimuthal component.

The resulting acoustic intensity field exhibits strong gradients for $x \le L_F$ (see Fig. 11 (b) and Moudjed et al. [43]), with $L_F = 1.22/(4S) \approx 10$ the Fresnel length locating the transition between the near ($x < L_F$) and far ($x > L_F$) acoustic fields. Since (i) the beam and source radii are equal to 0.5 at $x = 0$ and (ii) the beam expands radially at a rate $S$ due to diffraction [31], an approximate description of the beam radius variation along $x$ is:

$$ R_{beam} = 0.5 + Sx. \qquad (6) $$

The acoustic intensity $\tilde{I}_x(x, r = 0)$ along the beam axis is decaying mainly due to this radial expansion effect but also from attenuation. The drop of $P_{ac}$ along the beam due to attenuation yields an exponential contribution to the decay of $\tilde{I}_x(x, r = 0)$. From Eq. (A16), the length scale for this exponential decay is thus defined by $L_\alpha = L/(2N) \approx 166.7$. The ratio $L_F/L_\alpha \approx 0.06$ is thus marginal, meaning that the near acoustic field structure is barely affected by attenuation. In our work, we choose to focus on the scaling of the jet velocity with $Gr_{ac}$ and thus keep $L_F/L_\alpha$ constant. For more details on the effect of $L_F/L_\alpha$ on the jet velocity near the source, the interested reader shall refer to Daru et al. [33].

### C. Computational methodology

The problem (Eq. (3)) is solved numerically using the open source spectral-element code Semtex [44] on a 2D mesh [45]. The exponential convergence of the solution is then obtained by increasing the polynomial degree $N_p$ of the expansion basis represented at the Gauss-Lobatto-Legendre points.

The computational grid used throughout this study is shown in Fig. 2 and comprises 880 quadrilateral elements. The axial direction is discretized by 44 elements, with localized refinements at both ends of the domain [46] to capture the flow gradients at the jet origin and the jet impingement. The grid stretching in the $x$-direction is chosen so that the maximum ratio between the lengths of two consecutive elements does not exceed 1.35. The radial direction is discretized using 20 elements, with one element for $r \in [0, 0.5]$. The element height is then continuously increased by a factor 1.3 in the ascending $r$-coordinates.

The acoustic force $Gr_{ac}\tilde{I}$ is computed at the grid points prior to every flow calculations. The integrals in Eqs. (5) are evaluated numerically by discretizing the transducer with a uniform grid made of $N_s$ points in both azimuthal and radial directions. The total number of point sources on the transducer is thus equal to $N_s^2$.

Eqs. (3b) and (3a) are initiated with $u = 0$, and are marched forward in time using the splitting scheme of Karniadakis et al. [47] that integrates the inertial terms explicitly and the viscous terms implicitly. As only the steady state is of interest, the numerical scheme is used in its first-order formulation. To ensure numerical stability, the time step $\Delta t$ is adjusted for each $Gr_{ac}$ so that the maximum CFL number in the domain remains below unity. For instance, for the $Gr_{ac} = 1000$ case, $\Delta t = 10^{-3}$, yielding a maximum CFL number of 0.68. A steady state is then considered to
<figure><img src="image_4.png" /><figcaption>FIG. 2. Discretization of the fluid domain with elements of polynomial degree $N_p = 1$. The main lobe of the unattenuated acoustic beam is represented by the green-shaded area, of approximate radius given from Eq. (6). The red dots located at $(x, r) = (400, 15)$, $(1000, 70)$ and $(1400, 20)$, are the points where the time series of $u$ and $p$ are recorded to monitor the convergence towards a steady state.</figcaption></figure>

be reached when the time series of the flow variables, recorded at different points of the domain (located by the red dots in Fig. 2), do not show variations with time up to a precision of five significant figures.

Grid sensitivity studies were made for both the discretization of the transducer and the discretization of the fluid domain, and are given in Appendix B. All the results presented in this document were obtained with $N_s = 300$ and $N_p = 8$. This choice of $N_s$ ensures a discrepancy of less than 1 % between the numerical evaluation of Eqs. (5) for $r = 0$ and the analytical on-axis intensity profile given by Eq. (A15). Regarding the discretization of the fluid domain, $N_p = 8$ yields a marginal error of only about 0.01 % on both the location and the amplitude of the maximum on-axis jet velocity.

### III. THE ROLE OF THE ACOUSTIC FORCE ATTENUATION

In the literature, the acoustic attenuation has frequently been discarded in the model of the acoustic force to derive analytical 1D-1C (one-dimensional, one-component) Eckart streaming velocity profiles within a weak attenuation assumption [2, 3, 21, 35, 41, 48–50]. We want to emphasize the importance of the contribution of the attenuation in the spatial structure of the acoustic force, even for cases where the jet observation length is small compared to $L_α$. For this purpose we will focus on numerical calculations of the on-axis velocity profile, the radius of the jet and the momentum flux of the jet, by computing Eqs. (5) considering two cases: $N/L = 0$ in Eq. (5d) (without attenuation, as it is often considered in the literature) and $N/L \neq 0$ (with attenuation, as it will be considered from Sec. IV onwards). We shall then compare the numerical results to the theoretical predictions of Lighthill [1] regarding the jet momentum flux.

We first compare on Fig. 3 (a), the on-axis jet velocity at $Gr_{ac} = 10^3$ for the model including attenuation (solid line) and discarding attenuation (dashed line). For $x < 8$ (i.e., $x/L_α < 5 \times 10^{-2}$), no significant difference between the two velocity profiles is observed. However, the discrepancy builds up further downstream. In the absence of attenuation in $\tilde{I}$, $u_x$ reaches a maximum at a greater distance from the transducer, and with a larger magnitude. The velocity drop downstream the peak is also greatly affected, as $u_x$ decays at a significantly slower rate with $x$. In the vicinity of the downstream wall, the observed velocity magnitudes typically differ by a factor 10 between the attenuated and non-attenuated cases.

The large discrepancy observed between the on-axis velocity profiles calculated with and without attenuation of $\tilde{I}$ is accompanied by a different radial spreading rate of the jet (Fig. 3 (b)). In the near-acoustic field region ($x \le L_F$), the jet radius $R_{jet}$ (defined as the radial coordinate where $u_x$ is 50 % of the on-axis velocity) remains nearly constant for both the attenuated and unattenuated cases. Beyond $L_F$, $R_{jet}$ increases similarly to $R_{beam}$ and up to $x \approx 2 \times 10^2$ for an unattenuated $\tilde{I}$, whereas $R_{jet}$ experiences a larger increase in the attenuated case. This difference is related to the jet velocity: when high (as for the unattenuated $\tilde{I}$), viscous diffusion weakly affects the radial spread of the jet. Thus, the shape of the jet tends to follow more closely the shape of the beam, which is, in turn, dictated by acoustic diffraction. By contrast, low velocities (as for an attenuated $\tilde{I}$), increase the role of viscous diffusion. This leads to a greater radial enlargement of the jet as $x$ increases, as observed in Fig. 3 (b).

From these numerical results, we can now study the evolution along $x$ of the axial jet momentum flux defined as $\int_0^{2\pi} \int_0^{D/2} u_x^2 r dr d\theta$. As illustrated in Fig. 4 for $Gr_{ac} = 10^3$, when $N/L \neq 0$, the axial jet momentum flux (filled symbols) saturates for $x \ge 500$ up to $x = 1300$. On the contrary, setting $N/L = 0$ leads to a linear increase with $x$ of the axial momentum flux (empty symbols). This is in agreement with the theoretical prediction by Lighthill [1] for an axisymmetric Eckart streaming jet driven by an attenuated sound beam in a semi-infinite medium: the evolution with $x$ of the axial momentum flux, then defined as $\int_0^{2\pi} \int_0^{\infty} u_x^2 r dr d\theta$, is equal to the integral of the acoustic force
<figure><img src="image_5.png" /><figcaption>FIG. 3. Effects of acoustic force attenuation in the expression of $\tilde{I}$ (Eqs. (5)) on steady flows computed for $Gr_{ac} = 10^3$. The cases of attenuated ($N/L = 0.003$) and unattenuated ($N/L = 0$) $\tilde{I}$ are plotted with solid and dashed lines, respectively. (a) Longitudinal profiles of the on-axis velocity $u_x$, with a detailed view close to the transducer shown in the inset. (b) Longitudinal profiles of the jet radius $R_{jet}$ defined as $u_x(x, R_{jet}) = 0.5 u_x(x, 0)$. The approximate beam radius (Eq. (6)) is shown as a black dotted curve. The purple and red vertical lines locate the Fresnel distance $L_F$ and the force attenuation distance $L_\alpha$, respectively. In (b), the vertical blue line locates the on-axis velocity peak position $x_p$ in the attenuated $\tilde{I}$ case.</figcaption></figure>

$$ \bullet \int_{0}^{2\pi} \int_{0}^{D/2} u_x^2 r dr d\theta \text{ with } N/L = 0.003 \qquad \circ \int_{0}^{2\pi} \int_{0}^{D/2} u_x^2 r dr d\theta \text{ with } N/L = 0 \\ -M_{th} = \frac{\pi L}{32N} Gr_{ac} (1 - e^{-2N \frac{x}{L}}) \qquad \lim_{\frac{N}{L} \to 0} M_{th} = \frac{\pi}{16} Gr_{ac} x $$

<figure><img src="image_6.png" /><figcaption>FIG. 4. Evolution of the momentum flux accumulated by the jet with the distance $x$ from the transducer for $Gr_{ac} = 10^3$. The numerical results obtained when attenuation is accounted for ($N/L = 0.003$, Table I) or not ($N/L = 0$) in $\tilde{I}$ (Eqs. (5)) are represented by the filled and empty symbols, respectively. The numerical results are overlayed to the profiles of the injected momentum flux $M_{th}$, given by Eq. (7) and evaluated for $N/L = 0.003$ (solid blue line) and for $N/L \to 0$ (dashed blue line). The red vertical line locates the force attenuation length $L_\alpha = L/(2N)$. The horizontal black dashed line indicates the maximum momentum flux $M$ (from Lighthill [1], Eq. (2)). A detailed view of the momentum profiles, closer to the transducer, is shown in the inset.</figcaption></figure>

along the beam if the pressure gradient is neglected. This momentum injected by the acoustic forcing in the fluid will be denoted $M_{th}$, and obeys to the following expression:

$$ M_{th} = M \left(1 - e^{-2N \frac{x}{L}}\right) = M \left(1 - e^{-\frac{x}{L_\alpha}}\right). \qquad (7) $$

For $x/L_\alpha \gg 1$, the injected momentum flux saturates towards the value $M$. The corresponding maximum dimensional momentum flux to be acquired by the jet is then $P_{ac}/c$. On the other hand, in the limit of vanishing attenuation,
for $x/L_\alpha \ll 1$, Eq. (7) reduces to $\lim_{\frac{N}{L}\to 0} M_{th} \approx \frac{Mx}{L_\alpha} = \pi Gr_{ac}x/16$. In other words, beam attenuation directly dictates a length scale $L/(2N) = L_\alpha$ over which the jet momentum flux is created, and ignoring attenuation in the acoustic force leads to an unimpeded build-up of the momentum flux along the jet. These theoretical predictions are well recovered in the flow simulations (Fig. 4): for $x < 1100$, the theoretical and numerical profiles differ by less than 3%. For $x > 1300$, the larger discrepancies for both the attenuated and unattenuated cases are due to the recirculation imposed by the downstream wall, which is absent in the free-jet setting of Lighthill [1]. Finally, for very small values of $x$ (inset of Fig. 4), the discrepancy between the two numerical momentum flux profiles is weak. Nevertheless, the relative error $\epsilon$ on the momentum flux increases rapidly with $x$: for instance, $\epsilon \approx 0.05$ for $x = 0.1L_\alpha$ and $\epsilon \approx 0.6$ for $x = L_\alpha$. Hence, discarding the effects of attenuation in $\tilde{I}$ greatly affects the flow, even at $x \ll L_\alpha$.

To conclude, the maximum portion of acoustic power used to create momentum flux is limited to $P_{ac}/c$ which does not depend on $\alpha$. Attenuation of the beam thus only sets the length scale over which the momentum flux is created, not its maximum value; and only a fraction of $P_{ac}$ generates momentum. Neglecting the attenuation of $\tilde{I}$, even when weak, results in large errors on the velocity of the jet and on its radial enlargement, even at short range ($x \ll L_\alpha$). It is then of primary importance to include the attenuation while modeling the acoustic force even for small $x/L_\alpha$. From now, all the results are obtained for an attenuated expression of $\tilde{I}$.

## IV. CHARACTERIZATION OF LONG-RANGE ACOUSTIC STREAMING JETS

We shall now thoroughly characterize the streaming jet. A typical illustration of the studied velocity fields is shown in Fig. 5. The flow is steady, axisymmetric and laminar. The evolution of the jet velocity along its axis is characterized

<figure><img src="image_7.png" /><figcaption>FIG. 5. Typical steady velocity field obtained for $Gr_{ac} = 5 \times 10^3$. Top: Map of velocity magnitude together with streamlines. The approximate beam radius, defined by Eq. (6), is displayed in blue, and the inset on the right focuses on the flow region close to the transducer. Bottom: longitudinal profile of $u_x$ on the jet axis on which the Fresnel distance $L_F = 1.22/(4S)$, the on-axis velocity peak position $x_p$ and the acoustic force attenuation distance $L_\alpha = L/(2N)$ are reported. The right plot focuses on the region framed in black in the left plot.</figcaption></figure>

by a strong acceleration near the transducer (yellow and red regions in the bottom plot of Fig. 5). The irregularities witnessed in the on-axis velocity profile for $1 \le x \le 10$ are caused by the large longitudinal gradients of the force in the acoustic near field (Fig. 11 (b)) and were discussed in previous works [27, 43]. On the axis, the jet velocity then reaches a peak value, whose position $x_p$ does not coincide with either of the acoustic quantities $L_F$ or $L_\alpha$. At $x_p$, the Reynolds number is maximum and is equal to 150.9 for $Gr_{ac} = 5 \times 10^3$. Downstream $x_p$, the on-axis jet velocity decays smoothly as the distance from the transducer increases (blue region in the bottom plot of Fig. 5). Despite flowing over large distances compared to $L_\alpha$ (i.e., up to distances over which the sound beam has been nearly entirely attenuated), the jet impinges the downstream wall with a velocity of nearly 10% of its peak value. The jet
impingement occurring at $x = L$ then gives rise to a large vortex of weak velocity amplitude. In the following, we shall only consider the jet itself and each of the flow features identified on the axial jet velocity profile, from the fluid region close to the transducer up to the downstream wall. We shall rely on an order-of-magnitude approach [51] (i) to determine the approximate magnitude of the on-axis velocity $u_x$ in the jet regions identified in the bottom plots of Fig. 5, (ii) to explain the evolution of $u_x$ in these regions based on a force balance, and (iii) to determine how $u_x$ scales with the parameters of the problem.

## A. Inertia-dominated region of accelerating fluid

We shall now investigate the changes of jet velocity when $Gr_{ac}$ is varied from $10^3$ to $5 \times 10^4$ while keeping all other parameters constant. We shall first focus on the $0.5 \le x \le x_p$ region over which the jet accelerates. We will not consider the region $x \le 0.5$ (shown in yellow in the inset of the bottom plot in Fig. 5): it was already studied in Ref. [43] and is extremely small compared the length of the fluid domain.

For $0.5 \le x \le x_p$, former studies have shown that the jet accelerates due to a balance between the acoustic and inertia forces [31, 52]. In these previous works, the authors considered that since the length of the acceleration region is small compared to $L_\alpha$, it was legitimate to neglect both the attenuation and the beam diffraction [33]. The order of magnitude $U_{x,i}$ of the inviscid on-axis velocity was then given by [31]:

$$U_{x,i} \sim \frac{1}{2} \sqrt{Gr_{ac} x}. \qquad (8)$$

Equation (8) fairly describes the increase of $u_x$ on $0.5 \le x \le x_p$ as can be seen on Fig. 6 (a), in which an exact description of $u_x(x, r = 0)$ would give a straight horizontal line at $u_x/U_{x,i} = 1$. However, it either overestimates $u_x$ by nearly 30% for $Gr_{ac} = 10^3$, or instead underestimates it by 55% for $Gr_{ac} = 5 \times 10^4$. This discrepancy on the dependency on $Gr_{ac}$ between Eq. (8) and the numerical profiles may be attributed to viscosity effects, neglected in Ref. [31]. An order of magnitude accounting for viscous forces can be obtained from the first component of the

<figure><img src="image_8.png" /><figcaption>FIG. 6. Longitudinal profiles of the on-axis jet velocity normalized by different scaling laws, from the source ($x = 0$) to the velocity peak position $x_p$ (black points). The profiles are shown for different acoustic Grashof numbers $Gr_{ac}$. (a) On-axis velocity normalized by the inviscid scaling law $U_{x,i}$ (Eq. (8), [31]). (b) On-axis velocity normalized by the improved scaling law accounting for viscosity effects (Eq. (11)). As the $u_x$ is normalized in both cases by the derived scaling law, an exact description of the profile would appear as a straight horizontal line at a unit ordinate.</figcaption></figure>

Navier-Stokes equation, Eq. (3b). On $0.5 \le x \le x_p$, the jet is almost purely axial and the order of magnitude of the inertia term in Eq. (3b) is given by $u_x \partial u_x / \partial x$ [31]. Furthermore, the pressure forces are negligible due to the large cavity diameter $D$. From an order of magnitude standpoint, the Navier-Stokes equations (3b) thus give:

$$\frac{1}{2} \frac{\partial (U_{x,v}^2)}{\partial x} \approx -\frac{U_{x,v}}{R_{jet}^2} + Gr_{ac} \tilde{I}_x, \qquad (9)$$

where $U_{x,v}$ is an order of magnitude of $u_x(x, r = 0)$ accounting for viscous forces, and $\tilde{I}_x = \tilde{I} \cdot e_x$. From the acoustic power definition (Eq. (A14)) and its attenuation along a beam of varying cross-section:

$$\tilde{I}_x \sim \frac{e^{-\frac{x}{L_\alpha}}}{16R_{beam}^2}. \qquad (10)$$
Since both diffraction and attenuation are weak over the jet region of interest (i.e., $S \ll 1$ and $x \ll L_\alpha$ respectively), the spatial variations of $\tilde{I}_x$ caused by attenuation and diffraction can be neglected in Eq. (10). Besides, Fig. 3 (b) shows that $R_{jet} \sim R_{beam} \sim 1/2$ near the transducer. Further taking $\frac{\partial(U_{x,v}^2)}{\partial x} \sim U_{x,v}^2/x$, Eq. (9) then reduces to a quadratic polynomial of positive root

$$U_{x,v} \approx -4x + \frac{\sqrt{2}}{2} \sqrt{32x^2 + Gr_{ac}x}. \quad (11)$$

Equation (11) improves on Eq. (8) as it incorporates the effect of viscous forces while still neglecting acoustic attenuation and diffraction. That new scaling for $U_{x,v}$ is compared to the numerical data in Fig. 6 (b) in which, as in Fig. 6 (a), an exact description of $u_x(x, r = 0)$ would give a straight horizontal line at $u_x/U_{x,v} = 1$. A glance at this plot is enough to conclude that Eq. (11) significantly improves the prediction of the on-axis velocity profile $u_x(x, r = 0)$ upstream the velocity peak: the discrepancy between $U_{x,v}$ and $u_x(x, r = 0)$ is reduced to at most 25 %, confirming that the viscous forces play a non-negligible role in the $0.5 \le x \le x_p$ jet region. We consider this a fair improvement in the frame of the order-of-magnitude approach adopted here.

## B. Scaling of the on-axis velocity peak position and amplitude

The region of jet acceleration is followed by a velocity peak $u_{x,p}$ at $x_p$. Both $x_p$ and $u_{x,p}$ display a weak dependency on $Gr_{ac}$ that we shall now investigate. At $x = x_p$ and $r = 0$, both $\partial u_x/\partial x$ and $\partial u_x/\partial r$ are zero, causing the left-hand side of the steady form of the Navier-Stokes equation (3b) to vanish. For any $Gr_{ac}$, $u_{x,p}$ is thus controlled by a local balance between the acoustic and viscous forces:

$$2 \frac{U_{x,p}}{R_{jet,p}^2} \sim Gr_{ac} \tilde{I}_x(x, r = 0), \quad (12)$$

where $U_{x,p}$ and $R_{jet,p}$ are orders of magnitude of the velocity amplitude and jet radius at $x_p$, respectively. Based on the definition of the axial momentum flux and Eqs. (7), $R_{jet,p}$ can be related to $U_{x,p}$ through:

$$U_{x,p}^2 R_{jet,p}^2 \sim \frac{L_\alpha}{32} Gr_{ac} \left(1 - e^{-\frac{x_p}{L_\alpha}}\right). \quad (13)$$

We thus assume that $D$ is sufficiently large so that the sidewall does not affect $x_p$. In the present work, the on-axis velocity peak occurs in the acoustic far field where the acoustic intensity asymptotically evolves as $\tilde{I}_x(x, r = 0) \sim (\frac{1.22\pi}{8S})^2 \frac{e^{-\frac{x}{L_\alpha}}}{x^2}$ (see Eqs. (A16)). Assuming $x_p \ll L_\alpha$, both that far-field expression of $\tilde{I}_x$ and Eq. (13) can be simplified using Taylor's expansions so that the viscous vs acoustic force balance (12) becomes:

$$U_{x,p} \sim \left( \frac{1.22\pi Gr_{ac}}{64S} \right)^{\frac{2}{3}} x_p^{-\frac{1}{3}}. \quad (14)$$

Finally, substituting the purely inertial velocity estimate $U_{x,i}$ (Eq. (8)) for $U_{x,p}$ in Eq. (14) yields:

$$x_{p_i} \sim \left[ \left( \frac{1.22\pi}{8S} \right)^4 \frac{Gr_{ac}}{64} \right]^{\frac{1}{5}}, \quad (15a)$$

$$U_{x,p_i} \sim \frac{1}{2} \left[ \frac{1}{64} \left( \frac{1.22\pi}{8S} \right)^4 \right]^{\frac{1}{10}} Gr_{ac}^{\frac{3}{5}}. \quad (15b)$$

Equation (15a) is plotted as a blue line in Fig. 7 (a). Neglecting the viscosity effects in the fluid region upstream the velocity peak is seen to lead to overestimate the increase of $x_p$ with $Gr_{ac}$ observed in the simulations (represented as filled symbols). This leads to an error on $x_p$ ranging from 1 % to 46 % over the $Gr_{ac}$ range.

As seen in section IV A, for higher values of $Gr_{ac}$, Eq. (11) is expected to give a better estimate of $U_{x,p}$ than Eq. (8), then an estimate of the peak position $x_{p_v}$ is obtained numerically from:

$$4x_{p_v}^{\frac{4}{3}} - \frac{\sqrt{2}}{2} \left( 32x_{p_v}^{\frac{8}{3}} + Gr_{ac}x_{p_v}^{\frac{5}{3}} \right)^{\frac{1}{2}} + \left( \frac{1.22\pi Gr_{ac}}{64S} \right)^{\frac{2}{3}} = 0. \quad (16)$$
<figure><img src="image_9.png" /><figcaption>FIG. 7. Evolution of the on-axis velocity peak position $x_p$ and amplitude $u_{x,p}$ with $Gr_{ac}$. (a) Comparison between the peak positions observed in the numerical simulations (black filled symbols), and the predictions when neglecting viscous forces ($x_{p,i}$, blue line, Eq. (15a)) and considering them ($x_{p,v}$, red line, Eq. (16)) in the region upstream $x_p$. As the velocity peak lies in the far-acoustic field, both $x_{p,i}$ and $x_{p,v}$ are determined using the far-field approximation (A16) for the on-axis intensity $\tilde{I}_x$. The quantity represented on the x-axis is chosen to assess the validity of Eq. (15a). (b) Comparison between the peak amplitudes observed in the numerical simulations (black filled symbols), the experiments from Dentry et al. [30] (empty symbols) and the theoretical predictions (solid lines) when viscous forces upstream $x_p$ are neglected ($U_{x,p_i}$, blue, Eq. (15b)) and taken into account ($U_{x,p_v}$, red, Eqs. (16) and (14)). The quantity plotted on the x-axis is chosen to assess the validity of Eq. (15b). In (b), the data of Dentry et al. [30] correspond to a jet of rectangular cross section. The beam driving the jet is radiated by a rectangular Surface Acoustic Waves (SAW) device of height $H_s$ that varies with frequency. For the data points of Dentry et al. [30], $Gr_{ac}$ and $S$ are thus based on $H_s$ (instead of $D_s$ in our numerical simulations).</figcaption></figure>

The corresponding velocity magnitude $U_{x,p_v}$ at $x_{p_v}$ is then evaluated using either Eq. (14) or (11) together with the obtained value of $x_{p_v}$.

The obtained prediction $x_{p_v}$, plotted as a red line in Fig. 7 (a), appears as a better estimate of the peak position for larger values of the scaling parameter. At lower $Gr_{ac}$, that discrepancy is however larger. The reason is that $R_{jet}$ varies slightly with $Gr_{ac}$. For decreasing values of $Gr_{ac}$, the jet velocity also decreases. As a consequence, the effect of viscous diffusion on the radial enlargement of the jet becomes more important: at $L_F$, $R_{jet}$ measured for $Gr_{ac} = 10^3$ is 64 % larger than for $Gr_{ac} = 5 \times 10^4$ (not shown). Thus, considering the same $R_{jet}$ for all $Gr_{ac}$ in Eq. (9) leads to an overestimation of the viscous forces at low $Gr_{ac}$.

The same loss of accuracy at low $Gr_{ac}$ is found for $U_{x,p_v}$, albeit to a lesser extent (Fig. 7 (b)). Although $U_{x,p_i}$ provides an acceptable discrepancy of 8.7 % with the velocity peak magnitudes observed in the numerical simulations, its variation with $Gr_{ac}$ is further improved when viscous forces are taken into account upstream the peak. The difference between $U_{x,p_v}$ and $u_{x,p}$ indeed falls below 3 % for $Gr_{ac}^{3/5}/S^{2/5} \ge 4 \times 10^2$. For comparison, we added the experimental data points of Dentry et al. [30] (empty symbols) to Fig. 7 (b). In their experiments, the authors used rectangular Rayleigh Surface Acoustic Waves (SAW) devices to drive streaming jets. Contrary to our numerical setup, the jets in Ref. [30] were forced by a beam of rectangular cross section and were therefore non-axisymmetric. These differences in the forcing and flow geometries inevitably affect the magnitudes of the forces to which the jet is subjected, causing the values of $u_{x,p}$ to differ. Nevertheless, the data points of Ref. [30] scale similarly with $Gr_{ac}$ and $S$: this proves that $u_{x,p}$ indeed satisfies a local balance between the acoustic and viscous forces, as assumed by the scaling laws (15b) and (14). Improving the match between these scaling laws and the $u_{x,p}$ of Ref. [30] would require developments specific to their peculiar jet geometry that are out of the scope of the present work.

### C. Viscosity-dominated flow region downstream the velocity peak

We shall finally consider the fluid region located downstream $x_p$, where the on-axis jet velocity decays smoothly before reaching the wall at $x = L$ (blue region in the bottom plot of Fig. 5). Since inertial effects become negligible near $x_p$, the on-axis velocity $U_{x,d}$ in the downstream region results from the balance between the acoustic and the viscous forces, so using the definition of the acoustic power (Eq. (10)):

$$ \frac{U_{x,d}}{R_{jet}^2} \sim \frac{Gr_{ac}}{16R_{beam}^2} e^{-\frac{x}{L_{\alpha}}} \qquad (17) $$
In Ref. [31] a similar longitudinal rate of change is assumed for $R_{jet}$ and $R_{beam}$, so the variations of $U_{x,d}$ could only arise from the attenuation of the beam. This is not observed here (Fig. 3 (b)), and beam attenuation is too weak to be the sole reason for the decay of $u_x$ downstream $x_p$ (Fig. 5). It is thus of primary importance to rely on a good estimate of the axial jet profile to assess viscosity effects. Contrary to Ref. [31], we do not assume $R_{jet} = R_{beam}$, but instead base $R_{jet}$ on equation (13) describing the build up of the jet momentum flux along $x$. Equation (17) then gives:

$$U_{x,d} \sim Gr_{ac}^{\frac{2}{3}}X^{\frac{1}{3}}, \quad \text{with} \quad X = \frac{L_{\alpha}}{512R_{beam}^2}e^{-\frac{x}{L_{\alpha}}} \left(1 - e^{-\frac{x}{L_{\alpha}}}\right). \qquad (18)$$

As the considered flow region is located downstream $L_F$, the enlargement of the beam due to diffraction is significant and needs to be taken into account in the definition of $X$. This means that $R_{beam}$ can no longer be considered as constant, but shall instead follow the linear expansion given by Eq. (6). Thus, from Eq. (18), the longitudinal variations of $u_x$ are caused by both diffraction and attenuation.

As shown in Fig. 8, $U_{x,d}$ (plotted as a thick straight black line) correctly captures the slope of the on-axis $u_x$ (solid colored curves) on $x_p \le x \le L_\alpha$ and for all $Gr_{ac}$. Besides, $U_{x,d}$ differs from $u_x$ by a factor at most 3, meaning that Eq. (18) also correctly estimates the order of magnitude of $u_x$ on that $x$ interval.

### D. Negligible acoustic force region far downstream the velocity peak

Beyond $L_\alpha$, the magnitude of the acoustic force relative to inertia decreases due to attenuation, causing $u_x$ to deviate from $U_{x,d}$. Equation (18) is thus unable to describe the on-axis $u_x$ over the entire $x \ge x_p$ jet region, and another approach is needed to recover the on-axis profile of $u_x$ at $x \gg L_\alpha$.

At $x \gg L_\alpha$, most of the energy of the acoustic beam has been attenuated, meaning that the forcing is negligible whatever the $Gr_{ac}$. As the pressure forces also are negligible, the jet necessarily satisfies a balance between inertial and viscous forces at $x \gg L_\alpha$ according to the Navier-Stokes equations (3b). Furthermore, as shown in Sec. III, the beam attenuation causes the axial momentum of the acoustic streaming jet to asymptotically reach the constant value $M = \pi LGr_{ac}/(32N) = \pi L_\alpha Gr_{ac}/16$. The situation thus becomes similar to the case of a straight jet of constant axial momentum flux $M$ emerging from an orifice and described by the boundary layer equation:

$$u_x \frac{\partial u_x}{\partial x} + u_r \frac{\partial u_x}{\partial r} = \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial u_x}{\partial r} \right), \qquad (19)$$

<figure><img src="image_10.png" /><figcaption>FIG. 8. On-axis velocity profiles downstream the velocity peaks $x_p$, which are represented by the black points. The profiles are shown for different $Gr_{ac}$. The quantity plotted on the $x$-axis includes the longitudinal variations of the acoustic force caused by both diffraction and attenuation of the beam (see the right-hand side of Eq. (18)). The red vertical line locates the force attenuation distance $L_\alpha = L/(2N)$. The scaling law (18) for the on-axis velocity $U_{x,d}$ at $x_p \le x \le L_\alpha$ is represented by the thick straight black line. The colored solid lines correspond to the numerical data. The dashed curves represent the self-similar jet solution $U_{x,S}$ of Schlichting [53] given by Eq. (20) and based on the maximum axial momentum flux $M = \pi LGr_{ac}/(32N) = \pi L_\alpha Gr_{ac}/16$ asymptotically acquired by the jet.</figcaption></figure>
for which a divergence-free self-similar solution is [53]:

$$u_x = \frac{3 M}{8 \pi x} \frac{1}{(1 + \zeta^2)^2},$$

where

$$\zeta = \frac{1}{8} \sqrt{\frac{3M}{\pi}} \frac{r}{x}$$

is the self-similar variable. The on-axis velocity of the streaming jet is then:

$$U_{x,S} = \frac{3L_{\alpha}Gr_{ac}}{128x}, \qquad (20)$$

which is thus expected to accurately describe the on-axis jet velocity as $x \to +\infty$.

The scaling law (20) for $U_{x,S}$ is plotted as dashed curves for all $Gr_{ac}$ in Fig. 8 alongside the numerical velocity profiles (solid curves). As expected, the decrease of the acoustic forcing relative to inertia as $x$ increases (or equivalently, as $X$ reduces) causes the numerical profiles to come closer to $U_{x,S}$. This strong agreement confirms the suitability of $U_{x,S}$ for describing the on-axis jet velocity as $x \to +\infty$.

## E. Global picture of the streaming velocity profile on the jet axis

All the scaling laws derived throughout this work for a steady and axisymmetric Eckart streaming jet are summarized in Fig. 9 in the $Gr_{ac} = 5 \times 10^3$ case. The scaling law given in Ref. [43] and describing the jet formation at $x \le 0.5$ has also been plotted. Downstream this small region dominated by 3D effects, our quasi-one dimensional approach

<figure><img src="image_11.png" /><figcaption>FIG. 9. Longitudinal profile of the on-axis axial velocity $u_x$ computed for $Gr_{ac} = 5 \times 10^3$ (solid black line), together with the scaling laws (dashed lines) and the force balance in each jet region. The scaling law in the yellow region is taken from Moudjed et al. [43]. The predicted peak position $x_{p,v}$ (Eq. (16)) and amplitude $U_{x,p_v}$ (from Eq. (14) with $x_{p_v}$) are represented by the dash-dotted blue vertical line and the black star, respectively. The solid blue and red vertical lines respectively locate the actual on-axis velocity peak $x_p$ and the force attenuation distance $L_\alpha = L/(2N)$.</figcaption></figure>

sheds lights on the following successive regions:

(i) First, the acoustic force strongly accelerates the fluid near the source (red region in Fig. 9). In this region of the jet, the acoustic forcing is essentially balanced by inertial forces [28, 31, 52]. However, for the investigated values of $Gr_{ac}$, viscous forces are not entirely negligible, causing a discrepancy of nearly 50 % between the simulated profiles and the scaling law of Moudjed et al. [31]. Accounting for the viscous friction forces leads to Eq. (11), as proposed in Sec. IV A. This new scaling law clearly improves the velocity scaling on $Gr_{ac}$ and reduces to 25 % the error on the predicted velocity. Finally, the oscillations of $u_x$ observed in Fig. 11(b) are directly correlated to the strong longitudinal gradients of the acoustic force in the near-acoustic field [43].
(ii) The acceleration region is followed by a velocity peak on the jet axis. Beam diffraction causes the acoustic force on the axis to significantly decrease locally, so that the acceleration becomes weaker, up to the point where inertia completely vanishes at the velocity peak position $x_p$. By balancing the acoustic and viscous forces, we obtained a scaling law for $x_p$ and the peak velocity magnitude $u_{x,p}$. The corresponding star on Fig. 9 has been plotted using the derived laws given by Eqs. (16) and (14). These two relations clearly allow to tailor streaming jets to e.g. maximize the velocity at a given point in the fluid domain by tuning $S$ or $Gr_{ac}$ for instance.

(iii) Further downstream, as the beam keeps opening up due to diffraction, the acoustic force becomes locally weaker and is balanced by viscous forces. As described by Eq. (18), the beam diffraction and attenuation jointly drive the smooth velocity decay along the jet axis for $x_p \le x \lesssim L_\alpha$ (deep blue region of Fig. 9).

(iv) For $x \gg L_\alpha$, the energy of the acoustic beam is almost entirely attenuated. Thus, the momentum balance reduces to an equilibrium between viscous and inertial forces. Besides, the axial momentum of the jet saturates to the constant value $M = \pi LGr_{ac}/(32N)$ (from Eq. (7)). As a consequence, the scaling law of the jet recovers that of a self-similar axisymmetric jet of constant momentum (Eq.(20)) adapted from [53] and represented by the dashed black line in Fig. 9.

V. CONCLUSION

We presented a numerical study of an Eckart-type acoustic streaming jet representative of former experiments [27, 29, 31, 38]. The flow develops in a wide and elongated cylindrical cavity entirely filled with a viscous liquid. The acoustic source is a plane circular transducer, of axis aligned with the cylinder axis, placed at one end of the cavity in direct contact with the fluid. It radiates an axisymmetric beam-shaped acoustic field yielding a body force that drives the flow. An important aspect of this work is that the characteristic attenuation length $L_\alpha$ of the acoustic force is smaller but comparable to the cavity length, so that the jet features both a region close to the source with prominent acoustic force, and regions further away where it is not. The diameter of the cavity is sufficiently large to avoid the influence of the lateral wall on the jet. The range of acoustic force magnitudes, measured by the acoustic Grashof number $Gr_{ac}$, is limited to moderate values, so that the resulting streaming jet remains axisymmetric, steady and laminar. Hence, this configuration provides a generic representation of applications where a transducer is directed into a large volume so as to act over a distance greater than the attenuation length but smaller than the size of the cavity (this would prove particularly useful in metallurgical and crystal growth applications [19, 54]). In this type of configuration, acoustic attenuation plays a central role for the modeling of the acoustic field even at small distance from the source. Indeed, the error incurred on the injected momentum flux when neglecting attenuation is already 10 % at a distance $0.2L_\alpha$ from the source, and further increases with the on-axis distance $x$ from the transducer. Despite being frequently used in the literature, the unattenuated beam approximation leads to significant errors: it severely overestimates the jet velocity along its axis, which in turn affects the flow structure by reducing the impact of viscous diffusion on the radial spreading of the jet.

To characterize the jet, several flow regions and features were identified based on the mechanisms driving the flow. For each of these regions, we derived scaling laws describing the on-axis velocity profile. The obtained relations are all summarized in Fig. 9 and section IV E. We thus give a complete picture of the different regimes encountered along an acoustic streaming jet that is several times longer than the attenuation length. These scaling laws offer an advantageous alternative to simulations for obtaining an accurate representation of the on-axis jet velocity profile. This paves the way for tailored applications informed by accurate predictions on acoustically-driven jets. For instance, we may use these scaling laws to optimize configurations in which several successive beam reflections are used to create a fully 3D forcing [20]. Matching the velocity peak position $x_p$ with a beam reflection on a wall would then maximize the velocity at the corresponding jet impingement, hence maximizing the 3D flow intensity and mixing efficiency [20]. A similar choice may also be made to maximize the shear stresses at the impingement to e.g. improve the mass transfer at a wall [23], or for surface cleaning [55]. The jet can also be optimized to favor transition to a chaotic or turbulent state for the purpose of increasing mixing by contactless means [56].

Finally, it would be interesting to see whether the different approaches used in the present work could be applied to the study of other flows forced by external beam-shaped fields. Despite obvious differences related to the type of waves driving the flow, Eckart streaming indeed shares common aspects with flows driven by the attenuation of laser light [57–59] or internal gravity wave beams [60–62]. These similarities may potentially make it possible to transpose some of the derivations from this paper to such flows.
# ACKNOWLEDGMENTS

This work was carried out as part of the BRASSOA project supported by the Institut Carnot Ingénierie@Lyon and by a Royal Society International Exchange grant (Ref. IES\R2\202212). The support from the PMCS2I of Ecole Centrale de Lyon for the numerical calculations is gratefully acknowledged. We would particularly like to thank both Laurent Pouilloux (Ecole Centrale de Lyon) and Alex Pedcenko (Coventry University) for their availability and for providing help at any stage of this project. We also would like to thank Hugh M. Blackburn for helping with the use of the spectral element code Semtex.

For the purpose of Open Access, a CC-BY public copyright licence has been applied by the authors to the present document and will be applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission.

## Appendix A: Derivation of the equations for the acoustic field radiated by a plane circular transducer

We shall derive the equations governing the attenuated acoustic field radiated by a plane circular transducer. For the sake of clarity, the derivation is carried out using dimensional equations. We discuss at the end of the appendix the nondimensionalization and normalization of the forcing field equations using the same length, mass and time scales as those of the hydrodynamics problem. We restrict ourselves to the case of linear acoustic propagation but include dissipation.

### 1. Governing equations of acoustics

Although acoustic streaming refers to an incompressible flow, it relies on the attenuation of sound waves, which are a compressible process [39]. Based on the assumptions listed in Sec. II B, the propagation of linear sound waves far from solid boundaries can be described by the irrotational compressible Navier-Stokes and continuity equations linearized around a resting fluid of density $\rho$, together with an isentropic equation of state [36, 37]:

$$ \rho \frac{\partial \mathbf{u}_{ac}}{\partial t} = -\nabla p_{ac} + \rho c^2 \tau \nabla (\nabla \cdot \mathbf{u}_{ac}), \qquad (A1a) $$

$$ \frac{\partial \rho_{ac}}{\partial t} + \rho \nabla \cdot \mathbf{u}_{ac} = 0, \qquad (A1b) $$

$$ p_{ac} = c^2 \rho_{ac}, \qquad (A1c) $$

where $p_{ac}$ and $\rho_{ac}$ are the acoustic fluctuations of pressure and density, respectively, $c$ is the speed of sound and $\tau = \frac{1}{\rho c^2} \left( \frac{4}{3}\mu + \eta \right)$. The propagation of linear attenuated sound waves is then governed by a wave equation for $\mathbf{u}_{ac}$:

$$ \frac{1}{c^2} \frac{\partial^2 \mathbf{u}_{ac}}{\partial t^2} - \nabla^2 \mathbf{u}_{ac} = \tau \frac{\partial}{\partial t} (\nabla^2 \mathbf{u}_{ac}), $$

which is obtained by (i) taking the time-derivative of Eq. (A1a), (ii) taking the gradient of Eq. (A1b), (iii) combining the obtained equations and (iv) making an extensive use of $\nabla \times \mathbf{u}_{ac} = \mathbf{0}$. For a monochromatic wave $\mathbf{u}_{ac} = \hat{\mathbf{u}}_{ac}e^{j\omega t}$, the above wave equation then reduces to a complex Helmholtz equation:

$$ \nabla^2 \hat{\mathbf{u}}_{ac} + k^2 \hat{\mathbf{u}}_{ac} = \mathbf{0}. \qquad (A2) $$

In Eq. (A2), the effects of attenuation appear through the imaginary part of the complex wavenumber $k$ defined as:

$$ k = k_a - j\alpha = \frac{\omega}{c\sqrt{1 + j\omega\tau}}, \qquad (A3) $$

where $j^2 = -1$ and $\omega\tau$ is a dimensionless group comparing the wavelength to the sound attenuation length (i.e., $\omega\tau = NS/(1.22L)$). The real and imaginary parts of $k$ are given by [36]:

$$ k_a = \frac{\omega}{c\sqrt{2}} \left[ \frac{\sqrt{1+(\omega\tau)^2+1}}{1+(\omega\tau)^2} \right]^{\frac{1}{2}} \underset{\omega\tau \ll 1}{\sim} \frac{\omega}{c} \qquad (A4) $$

$$ \alpha = \frac{\omega}{c\sqrt{2}} \left[ \frac{\sqrt{1+(\omega\tau)^2-1}}{1+(\omega\tau)^2} \right]^{\frac{1}{2}} \underset{\omega\tau \ll 1}{\sim} \frac{\omega^2\tau}{2c} = \frac{\omega^2}{2\rho c^3} \left( \frac{4}{3}\mu + \eta \right). \qquad (A5) $$
<figure><img src="image_12.png" /><figcaption>FIG. 10. Sketch of the circular transducer of diameter $D_s$ (gray), modeled as a plane disc made of an infinite number of point sources. Each source is centered on an elementary surface $dS$ and radiates spherical acoustic waves contributing to the net acoustic pressure and velocity (hence intensity) levels at $P$ in the fluid domain.</figcaption></figure>

Finally, under the assumptions (ii) and (iii) of Sec. II B, the momentum equation (A1a) can be rewritten to give:

$$j \frac{\rho c^2 k^2}{\omega} \hat{u}_{ac} = -\nabla \hat{p}_{ac}, \qquad (A6)$$

which can be seen as a generalized specific acoustic impedance. Equation (A6) proves particularly useful when computing the acoustic intensity $I$:

$$I = \frac{1}{2} \Re \{ p_{ac} \mathbf{u}_{ac}^* \}, \qquad (A7)$$

where * is the complex conjugate.

## 2. Integral relations for the acoustic intensity field radiated by a plane circular transducer

To obtain the integral expression for $p_{ac}$ and $\mathbf{u}_{ac}$ (and hence to infer $I$), we shall closely follow the derivation of Blackstock [37] and adapt it to account for sound attenuation. The idea is to model the vibrating surface $S_t$ of the transducer as a plane disc made of an infinite number of point sources. Each of these sources is centered on an elementary surface $dS$ (Fig. 10) and drives a net flow rate $Q_0$ through that surface. We shall thus first focus on the spherical pressure wave emitted by a single point source, and determine how the wave amplitude relates to $Q_0$. The acoustic pressure and velocity levels in the fluid domain are then obtained by summing up the spherical waves radiated by all point sources on $S_t$. We checked that the same pressure Rayleigh integral as Eq. (A10) can be obtained using Green's functions (see Ref. [63] for the derivation in the context of unattenuated acoustics).

A point source may be seen as a special case of a pulsating sphere of vanishing mean radius $R_S$. As its surface vibrates at a velocity $U_{R_S} e^{j\omega t} e_{d'}$, where $e_{d'}$ is the unit radial vector, the sphere produces spherical acoustic pressure waves defined as [37]:

$$p_{ac} = \frac{A}{d'} e^{j[\omega t - k(d' - R_S)]}, \qquad (A8)$$

with $d'$ the distance from the center of the sphere, and $A$ is a constant defining the initial wave amplitude. The pulsating sphere drives a volumetric flow rate of amplitude $Q_{R_S} = 4\pi R_S^2 U_{R_S}$ through its surface. From Eqs. (A8) and (A6), one obtains:

$$Q_{R_S} = -j \frac{4\pi A \omega (1 + jkR_S)}{\rho c^2 k^2}.$$

The strength $Q_0$ of a point source is then determined by taking the limit of $Q_{R_S}$ as $R_S \to 0$, allowing to relate $Q_0$ to $A$. Equation (A8) can thus be rewritten as:
$$p_{ac} = j \frac{\rho c^2 k^2 Q_0}{4\pi\omega d'} e^{j(\omega t - kd')} . \qquad (A9)$$

On $S_t$, $Q_0 = u_0 \cdot n_t dS$, with $u_0$ the vibrating velocity of the transducer and $n_t$ the normal unit vector. Because the superposition principle applies, $p_{ac}$ in the fluid domain is a linear combination of all the pressure waves emitted by each point source and is thus obtained by integrating Eq. (A9) over $S_t$:

$$p_{ac} = j \frac{\rho c^2 k^2 u_0}{2\pi\omega} \iint_{S_t} \frac{e^{j(\omega t - k \|d\|)}}{\|d\|} dS, \qquad (A10)$$

in which $Q_0$ has been doubled, since each point source radiates waves in the $x > 0$ semi-infinite domain only [37]. Note that the framework of Green's function, used with the image method for plane interface, together with a Neumann boundary condition gives a more rigorous justification for the factor 2 on $Q_0$ (see [64] for a complete derivation). As the transducer vibrates uniformly, the acoustic fields are axisymmetric and thus need to be calculated in a single constant-$\theta$ half-plane only. For simplicity, we shall consider $\theta = 0$, so that the vectorial distance $d$ between a point source at $(r', \theta')$ on $S_t$ and the point at $(x, r)$ where $p_{ac}$ is evaluated (Fig. 10) reduces to:

$$d = x e_x + [r - r' \cos(\theta')] e_r - r' \sin(\theta') e_\theta .$$

An integral expression for $u_{ac}$ is obtained using a similar approach. From Eqs. (A6) and (A9), the acoustic velocity induced by a single point source is

$$u_{ac} = \frac{Q_0}{4\pi d'^2} (1 + jkd') e^{j(\omega t - kd')} e_{d'} .$$

The acoustic velocity radiated by the transducer is then obtained by integrating over $S_t$:

$$u_{ac} = \frac{u_0}{2\pi} \iint_{S_t} \frac{(1 + jk \|d\|)}{\|d\|^2} \frac{d}{\|d\|} e^{j(\omega t - k \|d\|)} dS, \qquad (A11)$$

in which $Q_0$ has again been doubled for the same reason as in Eq. (A10). We would like to emphasize that the Rayleigh integral formalism gives easily integral relations for both $p_{ac}$ and $u_{ac}$, which make the numerical implementation quite efficient. Greens'function formalism on the other hand, cannot give a straightforward derivation of $u_{ac}$ as in Eq.((A11)). From Eq. (A7), it appears that the temporal variations of $p_{ac}$ and $u_{ac}$ do not affect $I$; we can thus drop $e^{j\omega t}$ from Eqs. (A10) and (A11).

### 3. Analytical expressions of the acoustic fields on the transducer axis

Except at particular points in the fluid domain, the integral expressions for $p_{ac}$ and $u_{ac}$ derived in the previous section cannot be computed analytically in the entire $x > 0$ domain, and shall thus be evaluated numerically. It is nevertheless possible to obtain analytical expressions for $p_{ac}$, $u_{ac}$ and $I$ on the acoustic axis; these expressions prove useful for checking the numerical computations of Eqs. (A10) and (A11).

We shall first evaluate $p_{ac}$ at $r = 0$, where the attenuated Rayleigh integral (A10) is:

$$p_{ac}(x, r = 0) = j \frac{\rho c^2 k^2 u_0}{2\pi\omega} e^{j\omega t} \int_{\theta'=0}^{2\pi} \int_{r'=0}^{D_s/2} \frac{e^{-jk\sqrt{x^2+r'^2}}}{\sqrt{x^2+r'^2}} r' dr' d\theta' .$$

Setting $X = \sqrt{x^2 + r'^2}$, the above integral can be computed analytically to give:

$$\begin{aligned} p_{ac}(x, r = 0) &= \frac{\rho c^2 k u_0}{\omega} e^{j\omega t} \left( e^{-jkx} - e^{-jk\sqrt{x^2 + \frac{D_s^2}{4}}} \right) \\ &= \frac{\rho c^2 k u_0}{\omega} e^{j\omega t} e^{-j\frac{k}{2}\left(\sqrt{x^2 + \frac{D_s^2}{4}} + x\right)} \left[ e^{\frac{j\frac{k}{2}\left(\sqrt{x^2 + \frac{D_s^2}{4}} - x\right)} - e^{-j\frac{k}{2}\left(\sqrt{x^2 + \frac{D_s^2}{4}} - x\right)} \right] \\ &= j \frac{2\rho c^2 k u_0}{\omega} e^{j\omega t} e^{-j\frac{k}{2}\left(\sqrt{x^2 + \frac{D_s^2}{4}} + x\right)} \sin\left[\frac{k}{2}\left(\sqrt{x^2 + \frac{D_s^2}{4}} - x\right)\right]. \end{aligned}$$
Then, substituting the first equality of Eq. (A3) for $k$ in the above equation finally gives:

$$p_{ac}(x, r = 0) = \frac{2\rho c^2 (k_a - j\alpha) u_0 e^{-\frac{\alpha}{2} \left(\sqrt{x^2 + \frac{D_s^2}{4}} + x\right)}}{\omega} e^{j \left[\omega t - \frac{k_a}{2} \left(\sqrt{x^2 + \frac{D_s^2}{4}} + x\right)\right]} \\ \times \left(\sinh\left[\frac{\alpha}{2} \left(\sqrt{x^2 + \frac{D_s^2}{4}} - x\right)\right] \cos\left[\frac{k_a}{2} \left(\sqrt{x^2 + \frac{D_s^2}{4}} - x\right)\right]\right. \\ \left. + j \sin\left[\frac{k_a}{2} \left(\sqrt{x^2 + \frac{D_s^2}{4}} - x\right)\right] \cosh\left[\frac{\alpha}{2} \left(\sqrt{x^2 + \frac{D_s^2}{4}} - x\right)\right]\right). \quad (A12)$$

Equation (A12) is an original result accounting for attenuation. To assess its validity, we checked that taking the modulus of Eq. (A12) with $\alpha = 0$ gives the same unattenuated $|p_{ac}|$ as in Ref. [36].

Contrary to $p_{ac}$, Eq. (A11) cannot be calculated analytically even for $r = 0$. However, $u_{ac}(x, r = 0)e_x$ can be computed using Eq. (A6), in which $\nabla p_{ac}$ is calculated from Eq. (A12). The resulting expression of $u_{ac,x}$ can finally be used with Eq. (A7) to get the on-axis acoustic intensity $I_x(x, r = 0)$:

$$I_x(x, r = 0) = \frac{\rho c^2 u_0^2 k_a}{2\omega} e^{-\alpha \left(\sqrt{x^2 + \frac{D_s^2}{4}} + x\right)} \left( \left[1 - \frac{x}{\sqrt{x^2 + \frac{D_s^2}{4}}}\right] \sinh\left[\alpha \left(\sqrt{x^2 + \frac{D_s^2}{4}} - x\right)\right] \right. \\ \left. + \left[\frac{x}{\sqrt{x^2 + \frac{D_s^2}{4}}} + 1\right] \left(\cosh\left[\alpha \left(\sqrt{x^2 + \frac{D_s^2}{4}} - x\right)\right] - \cos\left[k \left(\sqrt{x^2 + \frac{D_s^2}{4}} - x\right)\right]\right) \right. \\ \left. + \frac{\alpha}{k_a} \left(1 - \frac{x}{\sqrt{x^2 + \frac{D_s^2}{4}}}\right) \sin\left[k_a \left(\sqrt{x^2 + \frac{D_s^2}{4}} - x\right)\right] \right). \quad (A13)$$

### 4. Normalization and dimensionless forms of the equations for the acoustic forcing field

We shall finally now provide the dimensionless and normalized version of the equations we just derived.

Setting $\alpha = 0$ and $k_\alpha = \omega/c$ in Eqs. (A12) and (A13) shows that $p_{ac}$ and $I_x$ are bounded respectively by $p_{ac,max} = 2\rho c u_0$ and $I_{max} = 2\rho c u_0^2$ in the unattenuated limit. Since $p_{ac} = \rho c u_{ac}$ for a plane wave, this also implies $u_{ac,max} = 2u_0$. All these quantities can be related to the acoustic power $P_{ac}$ radiated by the transducer through:

$$P_{ac} = \iint_{S_t} \mathbf{I} \cdot \mathbf{e}_x \, dS = \frac{\rho c u_0^2}{8} \pi D_s^2 = \frac{I_{max} \pi D_s^2}{16}, \quad (A14)$$

in which the pressure at the transducer has been taken as $\rho c u_0$. Equations (5) are finally obtained (i) by normalizing Eqs. (A7), (A10) and (A11) by $I_{max}$, $p_{ac,max}$ and $u_{ac,max}$, respectively. (ii) The $\omega\tau \ll 1$ expressions of $k_a$ and $\alpha$ (Eq. (A4)) are then used and (iii) made dimensionless using $D_s$ as length scale. A similar process applied to Eq. (A13) yields the normalized on-axis intensity:

$$\tilde{I}_x(x, r = 0) = \frac{1}{4} e^{-\frac{N}{L} \left(\sqrt{x^2 + \frac{D_s^2}{4}} + x\right)} \left( \left[ \frac{x}{\sqrt{x^2 + \frac{D_s^2}{4}}} + 1 \right] \left( \cosh \left[ \frac{N}{L} \left( \sqrt{x^2 + \frac{1}{4}} - x \right) \right] - \cos \left[ \frac{2.44\pi}{S} \left( \sqrt{x^2 + \frac{1}{4}} - x \right) \right] \right) \right. \\ \left. + \left[ 1 - \frac{x}{\sqrt{x^2 + \frac{D_s^2}{4}}} \right] \sinh \left[ \frac{N}{L} \left( \sqrt{x^2 + \frac{1}{4}} - x \right) \right] + \frac{NS}{2.44\pi L} \left[ 1 - \frac{x}{\sqrt{x^2 + \frac{D_s^2}{4}}} \right] \right. \\ \left. \times \sin \left[ \frac{2.44\pi}{S} \left( \sqrt{x^2 + \frac{1}{4}} - x \right) \right] \right), \quad (A15)$$

which gives the far-field $(1/(4x^2) \ll 1)$ approximation:

$$\tilde{I}_x(x, r = 0) \sim \left(\frac{1.22\pi}{8S}\right)^2 \frac{e^{-\frac{x}{L_\alpha}}}{x^2}. \quad (A16)$$
Equations (A15) and (A16) are respectively plotted as a solid and a dash-dotted lines in Fig. 11 for the parameter values listed in Table I.

Appendix B: Grid sensitivity analysis

<figure><img src="image_13.png" /><figcaption>FIG. 11. Effect of the transducer grid density on the acoustic force calculations. (a) Relative error $e$ between the numerical and analytical values of $\tilde{I}_x(x, r = 0)$ at different points on the axis, as a function of the number $N_s$ of point sources discretizing the transducer. (b) Comparison between the numerical (red, $N_s = 300$) and analytical (solid black curve, given by Eq. (A15)) evaluations of $\tilde{I}_x(x, r = 0)$ over the domain length. The far-field approximation (A16) of $\tilde{I}_x(x, r = 0)$ (dashed black curve) is also displayed, and the vertical purple line locates the Fresnel distance $L_F$. In (b), the numerical computations of $\tilde{I}_x(x, r = 0)$ are made at the collocation points for $N_p = 8$.</figcaption></figure>

We validated the numerical computation of $\tilde{I}$ (and hence, the acoustic force field) and assessed the number $N_s$ of grid points on the transducer to be used for this purpose. The numerical evaluation of $\tilde{I}_x(x, r = 0)$ is compared to its analytical expression in Fig. 11 (a). The error on the numerically-calculated $\tilde{I}_x(x, r = 0)$ decreases rapidly to become lower than 1% for $N_s \ge 250$. Therefore, all intensity fields used in this work were computed with $N_s = 300$ to obtain an efficient balance between accuracy of the forcing field calculations (Fig. 11 (b)) and computational cost.

The sensitivity to the grid discretizing the fluid domain has been assessed by carrying out *p*-refinement. The strong intensity gradients in the near acoustic field are undersampled to keep the numerical calculations tractable (Fig. 11 (b)). Thus, varying $N_p$ does not only act on the quality of the flow structures, it also affects the amount of momentum that is actually injected into the fluid domain.

Table II shows the effect of *p*-refinement on the position $x_p$ and amplitude $u_{x,p}$ of the velocity peak on the jet axis, for the case of the highest forcing magnitude investigated. The errors on both $x_p$ and $u_{x,p}$ fall below 0.1% for $N_p \ge 6$. Despite their non-monotonic evolution with $N_p$, caused by $\tilde{I}$ being undersampled in the acoustic near field, these errors remain extremely small for further-refined grids. The results shown in the remainder of this document were obtained for $N_p = 8$, for which the effects of the acoustic near field undersampling on the flow are proven to be negligible.
<table><caption>TABLE II. Dependence of the steady on-axis velocity peak position $x_p$ and amplitude $u_{x,p}$ with the polynomial degree $N_p$ of the expansion basis. The calculations are made for $Gr_{ac} = 5 \times 10^4$, i.e., for the case of highest forcing magnitude, and the errors are computed relatively to the results obtained for $N_p = 9$. The values of $x_p$ are obtained by finding the point where $\partial u_x / \partial x$ vanishes on the axis, using the Lagrange polynomials of the expansion basis as interpolation functions.</caption><thead><tr><th>N<sub>p</sub></th><th>x<sub>p</sub></th><th>u<sub>x,p</sub></th><th>Relative error on x<sub>p</sub> (%)</th><th>Relative error on u<sub>x,p</sub> (%)</th></tr></thead><tbody><tr><td>4</td><td>23.187</td><td>736.05</td><td>2.2917</td><td>0.4584</td></tr><tr><td>5</td><td>23.702</td><td>732.68</td><td>0.1256</td><td>0.0014</td></tr><tr><td>6</td><td>23.748</td><td>732.45</td><td>0.0691</td><td>0.0334</td></tr><tr><td>7</td><td>23.732</td><td>732.79</td><td>0.0012</td><td>0.0125</td></tr><tr><td>8</td><td>23.729</td><td>732.76</td><td>0.0110</td><td>0.0084</td></tr><tr><td>9</td><td>23.731</td><td>732.70</td><td>—</td><td>—</td></tr></tbody></table>

[1] J. Lighthill, Acoustic streaming, J. Sound Vib. 61, 391 (1978).

[2] C. Eckart, Vortices and streams caused by sound waves, Phys. Rev. 73, 68 (1948).

[3] W. L. Nyborg, Acoustic streaming due to attenuated plane waves, J. Acoust. Soc. Am. 25, 68 (1953).

[4] S. M. Hagsäter, T. Glasdam Jensen, H. Bruus, and J. P. Kutter, Acoustic resonances in microfluidic chips: full-image micro-PIV experiments and numerical simulations, Lab Chip 7 (2007).

[5] T. Frommelt, M. Kostur, M. Wenzel-Schäfer, P. Talkner, P. Hänggi, and A. Wixforth, Microfluidic mixing via acoustically driven chaotic advection, Phys. Rev. Lett. 100, 034502 (2008).

[6] J. Friend and L. Y. Yeo, Microscale acoustofluidics: microfluidics driven via acoustics and ultrasonics, Rev. Mod. Phys. 83, 647 (2011).

[7] P. Muller and H. Bruus, Theoretical study of time-dependent, ultrasound-induced acoustic streaming in microchannels, Phys. Rev. E 92, 063018 (2015).

[8] M. Alghane, Y. Q. Fu, B. X. Chen, Y. Li, M. P. Y. Desmulliez, and A. J. Walton, Frequency effect on streaming phenomenon induced by Rayleigh surface acoustic wave in microdroplets, J. Appl. Phys. 112, 10.1063/1.4758282 (2012).

[9] A. Riaud, M. Baudoin, O. Bou Matar, J.-L. Thomas, and P. Brunet, On the influence of viscosity and caustics on acoustic streaming in sessile droplets: an experimental and a numerical study with a cost-effective method, J. Fluid Mech. 821, 384 (2017).

[10] S. Cleve, M. Guédra, C. Mauger, C. Inserra, and P. Blanc-Benon, Microstreaming induced by acoustically trapped, non-spherically oscillating microbubbles, J. Fluid Mech. 875, 597 (2019).

[11] A. A. Doinikov, G. Regnault, C. Mauger, P. Blanc-Benon, and C. Inserra, Acoustic microstreaming produced by two interacting gas bubbles undergoing axisymmetric shape oscillations, J. Fluid Mech. 931, A19 (2022).

[12] M. Fauconnier, C. Mauger, J.-C. Béra, and C. Inserra, Nonspherical dynamics and microstreaming of a wall-attached microbubble, J. Fluid Mech. 935, A22 (2022).

[13] D. Baresch, J.-L. Thomas, and R. Marchiano, Observation of a single-beam gradient force acoustical trap for elastic particles: Acoustical tweezers, Phys. Rev. Lett. 116, 024301 (2016).

[14] M. Baudoin and J.-L. Thomas, Acoustic tweezers for particle and fluid micromanipulation, Annu. Rev. Fluid Mech. 52, 10 (2019).

[15] J. Li, A. Crivoi, X. Peng, L. Shen, Y. Pu, Z. Fan, and S. A. Cummer, Three dimensional acoustic tweezers with vortex streaming, Commun. Phys. 4 (2021).

[16] O. Pätzold, K. Niemietz, R. Lantzsch, V. Galindo, I. Grants, M. Bellmann, and G. Gerbeth, The use of magnetic fields in vertical bridgman/gradient freeze-type crystal growth, Eur. Phys. J. Spec. Top. 220, 243 (2013).

[17] M. Cablea, K. Zaidat, A. Gagnoud, A. Nouri, G. Chichignoud, and Y. Delannoy, Multi-crystalline silicon solidification under controlled forced convection, J. Cryst. Growth 417, 44 (2015).

[18] M. Chatelain, V. Botton, M. Albaric, D. Pelletier, B. Cariteau, D. Abdo, and M. Borrelli, Mechanical stirring influence on solute segregation during plane front directional solidification, Int. J. Therm. Sci. 126, 252 (2018).

[19] S. Miralles, M. Albaric, V. Brizé, P. Guy, B. Vincent, J.-P. Garandet, J.-C. Willemetz, D. Henry, and V. Botton, In situ ultrasonic interface tracking for photovoltaic silicon directional solidification, J. Cryst. Growth 621, 127362 (2023).

[20] B. Vincent, S. Miralles, D. Henry, V. Botton, and A. Pothérat, Experimental study of a helical acoustic streaming flow, Phys. Rev. Fluids 9, 024101 (2024).

[21] W. Dridi, D. Henry, and H. Ben Hadid, Influence of acoustic streaming on the stability of melt flows in horizontal Bridgman configurations, J. Cryst. Growth 310, 1546 (2008).

[22] S. Absar, P. Pasumarthi, and H. Choi, Numerical and experimental studies about the effect of acoustic streaming on ultrasonic processing of metal matrix nanocomposites (MMNCs), J. Manuf. Process. 28, 515 (2017).

[23] N. El Ghani, S. Miralles, V. Botton, D. Henry, H. Ben Hadid, B. Ter-Ovanessian, and S. Marcelin, Acoustic streaming enhanced mass transfer at a wall, Int. J. Heat Mass Transf. 172 (2021).
[24] M. C. Schenker, M. J. Pourquié, D. G. Eskin, and B. J. Boersma, PIV quantification of the flow induced by an ultrasonic horn and numerical modeling of the flow and related processing times, *Ultrason. Sonochem.* **20**, 502 (2013).

[25] G. Lebon, I. Tzanakis, K. Pericleous, D. Eskin, and P. S. Grant, Ultrasonic liquid metal processing: the essential role of cavitation bubbles in controlling acoustic streaming, *Ultrason. Sonochem.* **55**, 243 (2019).

[26] G. Lebon, G. Salloum-Abou-Jaoude, D. Eskin, I. Tzanakis, K. Pericleous, and P. Jarry, Numerical modelling of acoustic streaming during the ultrasonic melt treatment of direct-chill (DC) casting, *Ultrason. Sonochem.* **54**, 171 (2019).

[27] T. Kamakura, T. Sudo, K. Matsuda, and Y. Kumamoto, Time evolution of acoustic streaming from a planar ultrasound source, *J. Acoust. Soc. Am.* **100**, 132 (1996).

[28] H. Mitome, T. Kozuka, and T. Tuziuti, Effects of nonlinearity in development of acoustic streaming, *Japan. J. Appl. Phys.* **34**, 2584 (1995).

[29] H. Mitome, The mechanism of generation of acoustic streaming, *Electron. Commun. Japan* **81** (1998).

[30] M. B. Dentry, L. Y. Yeo, and J. R. Friend, Frequency effects on the scale and behavior of acoustic streaming, *Phys. Rev. E* **89**, 013203 (2014).

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

[32] J. Benhamou, B. Vincent, S. Miralles, M. Jami, D. Henry, A. Mezrhab, and V. Botton, Application of the lattice Boltzmann method to the study of ultrasound propagation and acoustic streaming in three-dimensional cavities: advantages and limitations, *Th. Comp. Fluid Dyn.* **37**, 725 (2023).

[33] V. Daru, B. Vincent, and M. Baudoin, High-speed and acceleration micrometric jets induced by GHz streaming: A numerical study with direct numerical simulations, *J. Acoust. Soc. Am.* **155**, 2470 (2024).

[34] L. K. Zarembo, Acoustic streaming, in *High-Intensity Ultrasonic Fields*, edited by L. Rozenberg (Springer Science, 1971) Chap. 3, pp. 137 – 199.

[35] O. V. Rudenko and S. I. Soluyan, *Theoretical Foundations of Nonlinear Acoustics* (Springer US, 1977).

[36] L. E. Kinsler, A. R. Frey, A. B. Coppens, and J. V. Sanders, *Fundamentals of Acoustics*, 4th ed. (John Wiley & Sons, Inc., 2000).

[37] D. T. Blackstock, *Fundamentals of Physical Acoustics* (John Wiley & Sons, Inc., 2000).

[38] V. Frenkel, R. Gurka, A. Liberzon, U. Shavit, and E. Kimmel, Preliminary investigations of ultrasound induced acoustic streaming using particle image velocimetry, *Ultrasonics* **39**, 153 (2001).

[39] P. Chassaing, R. A. Antonia, F. Anselmet, L. Joly, and S. Sarkar, *Variable Density Fluid Turbulence*, Vol. 69 (Springer Netherlands, 2002).

[40] F. Coulouvrat, On the equations of nonlinear acoustics, *J. Acoust.* **5**, 321 (1992).

[41] W. Dridi, D. Henry, and H. B. Hadid, Stability of buoyant convection in a layer submitted to acoustic streaming, *Phys. Rev. E* **81**, 056309 (2010).

[42] A. D. Pierce, Wave equation for sound in fluids with unsteady inhomogeneous flow, *J. Acoust. Soc. Am.* **87**, 2292 (1990).

[43] B. Moudjed, V. Botton, D. Henry, S. Millet, J.-P. Garandet, and H. Ben Hadid, Near-field acoustic streaming jet, *Phys. Rev. E* **91**, 033011 (2015).

[44] H. Blackburn, D. Lee, T. Albrecht, and J. Singh, Semtex: a spectral element-Fourier solver for the incompressible Navier-Stokes equations in cylindrical or cartesian coordinates, *Comput. Phys. Commun.* **245**, 106804 (2019).

[45] H. Blackburn and S. Sherwin, Formulation of a Galerkin spectral element-Fourier method for three-dimensional incompressible flows in cylindrical geometries, *J. Comput. Phys.* **197**, 759 (2004).

[46] M. Farrashkhalvat and J. Miles, *Basic Structured Grid Generation with an introduction to unstructured grid generation* (Butterworth Heinemann, 2003).

[47] G. E. Karniadakis, M. Israeli, and S. A. Orszag, High-order splitting methods for the incompressible Navier-Stokes equations, *J. Comput. Phys.* **97**, 414 (1991).

[48] O. Rudenko and A. Sukhorukov, Nonstationary Eckart Streaming and pumping of liquid in an ultrasonic field, *Acoust. Phys.* **44**, 565 (1998).

[49] W. L. Nyborg, Acoustic Streaming, in *Nonlinear Acoustics*, edited by M. F. Hamilton and D. T. Blackstock (Academic Press, 1998).

[50] H. Ben Hadid, W. Dridi, V. Botton, B. Moudjed, and D. Henry, Instabilities in the Rayleigh-Bénard-Eckart problem, *Phys. Rev. E* **86**, 016312 (2012).

[51] A. Bejan, *Convection Heat Transfer* (John Wiley & Sons, 2013).

[52] J. Orosco and J. Friend, Modeling fast acoustic streaming: steady-state and transient flow solutions, *Phys. Rev. E* **106**, 045101 (2022).

[53] H. Schlichting and K. Gersten, *Boundary-Layer Theory*, 8th ed. (Springer Berlin Heidelberg, 2000).

[54] G. Kozhemyakin, Influence of ultrasonic vibrations on the growth of semiconductor single crystals, *Ultrason.* **35**, 599 (1998).

[55] A. A. Busnaina and T. M. Elsawy, Post-CMP cleaning using acoustic streaming, *J. Elec. Mat.* **27**, 1095 (1998).

[56] G. Launay, T. Cambonie, D. Henry, A. Pothérat, and V. Botton, Transition to chaos in an acoustically driven cavity flow, *Phys. Rev. Fluids* **4**, 044401 (2019).

[57] A. Y. Savchenko, N. V. Tabiryan, and B. Y. Zel'dovich, Transfer of momentum and torque from a light beam to a liquid, *Phys. Rev. E* **56**, 4773 (1997).

[58] M. N. Shneider and V. V. Semak, Ponderomotive convection in water induced by a CW laser, *J. Appl. Phys.* **120**, 244902 (2016).

[59] R. Wunenburger, B. Issenmann, E. Brasselet, C. Loussert, V. Hourtane, and J.-P. Delville, Fluid flows driven by light
scattering, J. Fluid Mech. 666, 273 (2011).

[60] G. Bordes, A. Venaille, S. Joubaud, P. Odier, and T. Dauxois, Experimental observation of a strong mean flow induced by internal gravity waves, Phys. Fluids 24, 086602 (2012).

[61] N. Grisouard, M. Leclair, L. Gostiaux, and C. Staquet, Large scale energy transfer from an internal gravity wave reflecting on a simple slope, Procedia IUTAM 8, 119 (2013).

[62] T. Dauxois, S. Joubaud, P. Odier, and A. Venaille, Instabilities of internal gravity wave beams, Annu. Rev. Fluid Mech. 50, 131 (2018).

[63] A. D. Pierce, *Acoustics: An Introduction to Its Physical Principles and Applications* (McGraw-Hill Book Co, 1981).

[64] P. M. Morse and H. Feshbach, *Method of theoretical physics* (McGraw-Hill Book Co, 1953).