# Application of the lattice Boltzmann method to the study of ultrasound propagation and acoustic streaming in three-dimensional cavities: Advantages and limitations

Jaouad Benhamou, Bjarne Vincent, Sophie Miralles, Mohamed Jami, Daniel Henry, Ahmed Mezrhab, Valéry Botton

► To cite this version:

Jaouad Benhamou, Bjarne Vincent, Sophie Miralles, Mohamed Jami, Daniel Henry, et al.. Application of the lattice Boltzmann method to the study of ultrasound propagation and acoustic streaming in three-dimensional cavities: Advantages and limitations. Theoretical and Computational Fluid Dynamics, 2023, 37, pp.725-753. 10.1007/s00162-023-00676-9. hal-03943859v2

HAL Id: hal-03943859

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

Submitted on 16 Nov 2023

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Application of the lattice Boltzmann method to the study of
ultrasound propagation and acoustic streaming in three-
dimensional cavities: Advantages and limitations

Jaouad Benhamou¹*, Bjarne Vincent², Sophie Miralles², Mohammed Jami¹, Daniel Henry²*, Ahmed Mezrhab¹, and Valéry Botton²

¹Laboratoire de Mécanique & Energétique, Département de Physique, Faculté des Sciences, Université Mohammed Premier, 60000 Oujda, Morocco.

²Laboratoire de Mécanique des Fluides et d'Acoustique, CNRS/Université de Lyon, Ecole Centrale de Lyon/Université Lyon 1/INSA Lyon, ECL, 36 Avenue Guy de Collongue, 69134 Ecully Cedex, France.

*Corresponding authors.

E-mail addresses: [benhamou_jaouad1718@ump.ac.ma](mailto:benhamou_jaouad1718@ump.ac.ma) (J. Benhamou), [daniel.henry@ec-lyon.fr](mailto:daniel.henry@ec-lyon.fr) (D. Henry).

**Abstract**

The paper presents a three-dimensional numerical study of the acoustic streaming induced by the dissipation of ultrasounds during their propagation in the air. The waves are generated by a circular acoustic source positioned at the center of the left wall of a parallelepipedic cavity. The simulations are performed with the lattice Boltzmann method associated with the D3Q19 multiple relaxation time model. A validation of this model is first performed by comparing the numerical and analytical acoustic intensities along the central axis of the acoustic source. The main objective of this study is to use two different methods to calculate the acoustic streaming flow. The first method is the direct calculation of the mean velocity fields as the mean values of the instantaneous velocities. The second method is an indirect technique, which first calculates the acoustic streaming force and then injects this force into the numerical code to produce the streaming. A comparison between the results obtained by the two methods was carried out and a good agreement was found between them. These different investigations, rather new in three-dimensional configurations, have allowed us to discuss the advantages and limitations of the lattice Boltzmann approach to simulate real situations of wave propagation and acoustic streaming.

**Keywords:** 3D simulation; acoustic streaming; lattice Boltzmann method; ultrasound.
<table><thead><tr><th colspan="2">Nomenclature</th><th colspan="2">Greek symbols</th></tr></thead><tbody><tr><td>c</td><td>Lattice speed</td><td>α</td><td>Attenuation coefficient (1/m)</td></tr><tr><td>c<sub>i</sub></td><td>Discrete LBM velocities</td><td>γ</td><td>Ratio of specific heats</td></tr><tr><td>c<sub>lbm</sub></td><td>LBM speed of sound, 1/√3</td><td>δt</td><td>LBM time step, δt = 1</td></tr><tr><td>c<sub>ph</sub></td><td>Physical speed of sound</td><td>Δt</td><td>Physical time step (s)</td></tr><tr><td>d</td><td>Diameter</td><td>Δx</td><td>Physical space step (m)</td></tr><tr><td>e</td><td>Kinetic energy</td><td>λ</td><td>Wavelength (m)</td></tr><tr><td>e<sup>2</sup></td><td>Square of the kinetic energy</td><td>ν</td><td>Kinematic viscosity (m<sup>2</sup>/s)</td></tr><tr><td>f<sub>i</sub></td><td>Distribution function</td><td>μ</td><td>Dynamic viscosity (Pa·s)</td></tr><tr><td>F</td><td>Body force</td><td>ρ</td><td>Fluid density (Kg/m<sup>3</sup>)</td></tr><tr><td>g</td><td>Gravitational vector</td><td>ρ<sub>0</sub></td><td>Equilibrium density (Kg/m<sup>3</sup>)</td></tr><tr><td>I</td><td>Intensity</td><td>Ω<sub>i</sub></td><td>Collision operator</td></tr><tr><td>j</td><td>Impulsion</td><td></td><td></td></tr><tr><td>k</td><td>Angular wavenumber</td><td></td><td></td></tr><tr><td>K</td><td>Thermal conductivity</td><td></td><td></td></tr><tr><td>Lx</td><td>Cavity length</td><td></td><td></td></tr><tr><td>Ly</td><td>Cavity width</td><td></td><td></td></tr><tr><td>Lz</td><td>Cavity height</td><td></td><td></td></tr><tr><td>M'</td><td>Dimension of mass</td><td></td><td></td></tr><tr><td>M</td><td>Transformation matrix</td><td></td><td></td></tr><tr><td>M<sup>-1</sup></td><td>Inverse matrix of M</td><td></td><td></td></tr><tr><td>Ma</td><td>Mach number</td><td></td><td></td></tr><tr><td>m<sub>i</sub></td><td>Moment vector</td><td></td><td></td></tr><tr><td>m<sub>i</sub><sup>eq</sup></td><td>Equilibrium moment</td><td></td><td></td></tr><tr><td>p</td><td>Pressure</td><td></td><td></td></tr><tr><td>q</td><td>Energy flux</td><td></td><td></td></tr><tr><td>R<sub>gas</sub></td><td>Specific gas constant</td><td></td><td></td></tr><tr><td>S</td><td>Relaxation time matrix</td><td></td><td></td></tr><tr><td>t</td><td>Time</td><td></td><td></td></tr><tr><td>T</td><td>Period</td><td></td><td></td></tr><tr><td>u</td><td>Velocity along the x-axis</td><td></td><td></td></tr><tr><td>v</td><td>Velocity along the y-axis</td><td></td><td></td></tr><tr><td>w</td><td>Velocity along the z-axis</td><td></td><td></td></tr><tr><td>W<sub>i</sub></td><td>Discretization weights</td><td></td><td></td></tr><tr><td>x<sub>B</sub></td><td>Node boundary</td><td></td><td></td></tr></tbody></table>

# 1. Introduction

Acoustic or sound waves can be defined as pressure fluctuations that can exist in a compressible fluid. Their propagation can give rise to fundamental phenomena such as reflection, interference, and refraction. Generally, they are classified according to their vibration frequencies. This classification concerns audible waves of medium frequency, infrasound of low frequency, and ultrasound of high frequency [1]. The study of their propagation in fluids has been well-developed in the scientific community for many years. These waves have also a major interest in fundamental and applied research. They are used in important applications in various fields, particularly as ultrasound waves. Among the most important applications are their use in medicine for diagnosis and treatments [2–4], the exploitation
of the flows induced by the dissipation of their energy during their propagation in fluids in engineering, for example for the purification of photovoltaic silicon [5] or the separation of species in a mixture [6], their use in daily life for cleaning [7].

In view of these applications, the present work is devoted to the simulation of wave propagation with the objective of computing the flow caused by their propagation in a fluid medium. This phenomenon is known in the literature as acoustic streaming [8]. In general, two types of acoustic streaming flows are considered, Eckart streaming [9] and Rayleigh streaming [10]. In Eckart streaming, the fluid motion is produced by the dissipation of the energy of the traveling waves during their propagation in the core of the viscous fluid. In Rayleigh streaming, the dissipation occurs in the acoustic boundary layers at the solid walls [11]. Acoustic streaming has been studied theoretically, numerically, and experimentally for many years, in connection with many different applications. Different research teams investigated the characteristics of acoustic streaming and the way it is produced [12–14]. Acoustic streaming was also studied as a way to control droplets and bubbles without contact [15,16]. The ability of sound waves to remotely induce fluid flows is also of particular interest in the field of microfluidics, where the strong confinement makes efficient stirring challenging [17,18], and in industrial processes, such as direct-chill casting [19]. Acoustic streaming can also influence heat transfer and convection in heated cavities [20,21]. With the recent discovery that standing sound waves generate much higher velocity flows if the medium under study is inhomogeneous, Qiu et al. [22] focused on an experimental study of the characterization of Rayleigh acoustic streaming in an acoustic resonator with density and compressibility gradients. Finally, acoustic streaming was also recently used to mechanically act on cells in view of intracellular delivery [23,24].

In this paper, the simulation of a three-dimensional acoustic streaming flow is performed with the lattice Boltzmann method (LBM). This method is an alternative technique to more conventional computational fluid dynamics methods for the simulation of fluid flows. Contrary to the traditional approach based on the Navier-Stokes equations, the LBM method consists in discretizing the Boltzmann equation, corresponding to the statistical modeling of the dynamics of the fluid particles [25–27]. Through its mesoscopic nature, it has the high capacity to integrate the laws of microscopic or mesoscopic physics, allowing the reproduction of macroscopic laws at a very reasonable computational cost. As a result, it offers the possibility to simulate complex physical phenomena quite easily and is an object of interest for many researchers in numerical physics. Moreover, its important parallelization capabilities make it attractive for performing fast simulations on parallel computers.

The lattice Boltzmann approach has been well used to study wave propagation and acoustic streaming since the late 1990s. We can mention some of the various works that are more focused on acoustic streaming. Stansell and Created [28] numerically studied the acoustic streaming resulting from the interaction of acoustic waves with no-slip boundaries in a two-dimensional (2D) pipe using the lattice gas automaton fluid modeling method. Haydock and Yeomans [29] performed LBM simulations to determine the Rayleigh acoustic streaming produced by the interaction of an acoustic wave with a solid boundary. They also showed that LBM simulations can be easily used to obtain the streaming flow induced by the attenuation of a traveling wave propagating in the air [30]. Rafat et al. [31] used LBM simulations to study the streaming in a simplified thermoacoustic refrigerator. More recently, Tan et al. [32] coupled the lattice Boltzmann method with the finite difference technique to study the acoustic streaming driven by surface acoustic waves in microchannels and nanochannels. Li et al. [33] used the LBM to study the effect of the acoustic streaming flow on heat transfer in a phase change material. Finally, in previous works, our team also used the lattice Boltzmann technique based on the multiple
relaxation time (MRT) model to study 2D acoustic streaming induced by ultrasound propagation in water [34,35]. Note that all these studies simulate two-dimensional flows.

In this paper, the LBM approach is applied to simulate Eckart acoustic streaming in three dimensions (3D). The main goal of this simulation is to calculate acoustic streaming directly from the calculation of the average values of the instantaneous velocities and then to compare the results with those obtained from a more classical calculation involving the acoustic streaming force. Such studies are rather new in the literature, as we did not find previous works using the LBM to produce acoustic streaming in 3D with the acoustic point source method. Based on these results and those of our previous works [35,36], a discussion on the advantages and limits of the LBM technique to study acoustic streaming is proposed.

The paper is organized into five sections. After the introduction, the second section presents the LBM method. It details the different mathematical formulations of the LBM model used, provides a summary of the wave generation techniques in LBM simulations, explains how a force can be implemented in the LBM code, describes the boundary conditions, and finally gives the relations between the physical units and the LBM units. The third section is the main section devoted to the presentation of the results. After the description of the physical problem and the definition of the acoustic parameters, the ability of our numerical code to accurately simulate wave propagation is demonstrated, thanks to the comparison between the numerically calculated acoustic intensity on the central axis of the source and the corresponding analytical results. The acoustic streaming results, obtained by the two direct and indirect methods previously described, are then presented and discussed. The final section gives the general conclusions of these numerical studies.

## 2. Numerical method

Our simulations of steaming flows are performed with the lattice Boltzmann method. This method can be associated with three different models: SRT, TRT, and MRT models [25,26]. For the study of acoustic waves, the LBM-MRT model presents several advantages over the two other models, first in terms of accuracy [25] and then because of the possible use of small viscosity values [35,36]. Thus, the D3Q19-MRT model will be used in this study to calculate the acoustic density, velocities, pressure, and intensity. The mathematical formulation of the model, the different techniques to generate the waves, the implementation of the acoustic streaming force, and the link between the physical and LBM units are successively discussed in the present section.

### 2.1 D3Q19-MRT lattice Boltzmann model

In this study, the 3D continuous computational domain is discretized according to the D3Q19 model with a 19 points cubic lattice. In this lattice, the particles can propagate towards the centers of the side planes or the corners of the central planes in 18 possible directions (see Fig. 1). During the propagation, collisions can occur. These two phenomena are modeled by the following Boltzmann equation [25,26,36,37]:

$$f_i(x_i + c_i \delta t, t + \delta t) = f_i(x_i, t) + \Omega_i. \quad (1)$$

In the D3Q19 lattice, particles can take the following 19 LBM velocities $c_i = \begin{pmatrix} c_{x,i} \\ c_{y,i} \\ c_{z,i} \end{pmatrix}$ [36,38]:
$$c_i = c \begin{pmatrix} 0 & 1 & -1 & 0 & 0 & 0 & 0 & 1 & -1 & 1 & -1 & 1 & -1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & -1 & 0 & 0 & 1 & 1 & -1 & -1 & 0 & 0 & 0 & 0 & 1 & -1 & 1 & -1 \\ 0 & 0 & 0 & 0 & 0 & 1 & -1 & 0 & 0 & 0 & 0 & 1 & 1 & -1 & -1 & 1 & 1 & -1 & -1 \end{pmatrix} \quad (2)$$

<figure><img src="image_2.png" /><figcaption>Fig. 1. Illustration of the D3Q19-LBM lattice.</figcaption></figure>

LBM models with multiple relaxation times are based on the concept that each physical quantity relaxes to equilibrium in its proper relaxation time. This approach leads to define the collision operator as [25,26]

$$\Omega_i = M^{-1}S[m_i^{eq} - m_i]. \quad (3)$$

The vector $m$ of components $m_i$ is composed of 19 different macroscopic quantities, which are well discussed in references [36,38]. However, typically, researchers are only interested in calculating four quantities corresponding to the density $\rho$ and the momentums $j_x, j_y$, and $j_z$ along the $x, y$, and $z$ axes. $m$ can be represented mathematically as

$$m = (\rho, e, e^2, j_x, q_x, j_y, q_y, j_z, q_z, 3p_{xx}, 3\pi_{xx}, p_{ww}, \pi_{ww}, p_{xy}, p_{yz}, p_{zx}, m_x, m_y, m_z)^T. \quad (4)$$

The matrix $M^{-1}$ can be obtained by calculating the inverse of the transformation matrix $M$, which is defined as [36,38,39]
$$
M = \begin{pmatrix}
1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\
-30 & -11 & -11 & -11 & -11 & -11 & -11 & 8 & 8 & 8 & 8 & 8 & 8 & 8 & 8 & 8 & 8 & 8 & 8 \\
12 & -4 & -4 & -4 & -4 & -4 & -4 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\
0 & 1 & -1 & 0 & 0 & 0 & 0 & 1 & -1 & 1 & -1 & 1 & -1 & 1 & -1 & 0 & 0 & 0 & 0 \\
0 & -4 & 4 & 0 & 0 & 0 & 0 & 1 & -1 & 1 & -1 & 1 & -1 & 1 & -1 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 1 & -1 & 0 & 0 & 1 & 1 & -1 & -1 & 0 & 0 & 0 & 0 & 1 & -1 & 1 & -1 \\
0 & 0 & 0 & -4 & 4 & 0 & 0 & 1 & 1 & -1 & -1 & 0 & 0 & 0 & 0 & 1 & -1 & 1 & -1 \\
0 & 0 & 0 & 0 & 0 & 1 & -1 & 0 & 0 & 0 & 0 & 1 & 1 & -1 & -1 & 1 & 1 & -1 & -1 \\
0 & 0 & 0 & 0 & 0 & -4 & 4 & 0 & 0 & 0 & 0 & 1 & 1 & -1 & -1 & 1 & 1 & -1 & -1 \\
0 & 2 & 2 & -1 & -1 & -1 & -1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & -2 & -2 & -2 & -2 & -2 \\
0 & -4 & -4 & 2 & 2 & 2 & 2 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & -2 & -2 & -2 & -2 & -2 \\
0 & 0 & 0 & 1 & 1 & -1 & -1 & 1 & 1 & 1 & 1 & -1 & -1 & -1 & -1 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & -2 & -2 & 2 & 2 & 1 & 1 & 1 & 1 & -1 & -1 & -1 & -1 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -1 & -1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -1 & -1 & 1 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -1 & -1 & 1 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -1 & 1 & -1 & -1 & 1 & -1 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -1 & 1 & -1 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -1 & 1
\end{pmatrix} \tag{5}
$$

The equilibrium moment vector $m^{eq}$ is related to the equilibrium distribution function $f^{eq}$ by $m^{eq} = M f^{eq}$, with $f^{eq}$ expressed in its discrete form as [25]

$$f_i^{eq} = W_i \rho \left[ 1 + \frac{1}{c_{lbm}^2} \vec{c}_i \cdot \vec{V} + \frac{1}{2c_{lbm}^4} (\vec{c}_i \cdot \vec{V})^2 - \frac{1}{2c_{lbm}^2} |\vec{V}|^2 \right], \quad \text{with } i = 0, \dots, 18 \quad (6)$$

where $\vec{V}$ is the velocity vector ($\vec{V} = (\vec{u}, \vec{v}, \vec{w})$) and $W_i$ represents the discretization weights, which are given as:

$$W_i = \begin{cases} \frac{1}{3} & i = 0 \\ \frac{1}{18} & i = 1, \dots, 6 \\ \frac{1}{36} & i = 7, \dots, 18 \end{cases} \qquad (7)$$

The nineteen elements of $m_i^{eq}$ found after performing the calculations can be given in terms of density and momentums as follows [38,39]:

$$
\begin{align*}
m_0^{eq} &= \rho \\
m_1^{eq} &= -11\rho + 19 \frac{j_x^2 + j_y^2 + j_z^2}{\rho_0} \\
m_2^{eq} &= 3\rho - \frac{11j_x^2 + j_y^2 + j_z^2}{2\rho_0} \\
m_3^{eq} &= j_x \\
m_4^{eq} &= -\frac{2}{3}j_x \\
m_5^{eq} &= j_y \\
m_6^{eq} &= -\frac{2}{3}j_y \\
m_7^{eq} &= j_z \\
m_8^{eq} &= -\frac{2}{3}j_z
\end{align*}
$$
$$ 
\begin{aligned}
m_9^{eq} &= \frac{2j_x^2 - (j_y^2 + j_z^2)}{\rho_0} \\
m_{10}^{eq} &= -\frac{1}{2} \frac{2j_x^2 - (j_y^2 + j_z^2)}{\rho_0}
\end{aligned}
\qquad (8) $$

$$ 
\begin{aligned}
m_{11}^{eq} &= \frac{j_y^2 - j_z^2}{\rho_0} m_{12}^{eq} = -\frac{1}{2} \frac{j_y^2 - j_z^2}{\rho_0} \\
m_{13}^{eq} &= \frac{j_x j_y}{\rho_0} \\
m_{14}^{eq} &= \frac{j_y j_z}{\rho_0} \\
m_{15}^{eq} &= \frac{j_x j_z}{\rho_0} \\
m_{16}^{eq} &= 0 \\
m_{17}^{eq} &= 0 \\
m_{18}^{eq} &= 0
\end{aligned}
$$

The matrix S is diagonal. It defines each relaxation time $s_i$ of the corresponding macroscopic quantities. It can be written as a function of the nineteen $s_i$ as:

$$ S = \mathrm{diag}(s_0, s_1, s_2, s_3, s_4, s_5, s_6, s_7, s_8, s_9, s_{10}, s_{11}, s_{12}, s_{13}, s_{14}, s_{15}, s_{16}, s_{17}, s_{18}) \qquad (9) $$

The choice of $s_i$ values is based in particular on the stability of the D3Q19 model. The values used in this study are chosen according to the literature [36,38,40,41]:

$$ 
\left\{ 
\begin{array}{l}
s_0 = s_3 = s_5 = s_7 = 1.0, \quad s_1 = 1.19, \quad s_2 = s_{10} = s_{12} = 1.4 \\
s_9 = s_{11} = s_{13} = s_{14} = s_{15} = 1/(3\nu + 0.5) \\
s_4 = s_6 = s_8 = 1.2, \quad s_{16} = s_{17} = s_{18} = 1.6
\end{array}
\right. 
\qquad (10) $$

The solution of the lattice Boltzmann equation allows to determine the distribution function $f_i$. The density and velocities can then be calculated as:

$$ \rho = \sum_{i=0}^{18} f_i, \quad \rho \begin{pmatrix} u \\ v \\ w \end{pmatrix} = j = \sum_{i=0}^{18} f_i c_i = \sum_{i=0}^{18} f_i \begin{pmatrix} c_{x,i} \\ c_{y,i} \\ c_{z,i} \end{pmatrix}. \qquad (11) $$

The acoustic intensity is also an interesting quantity in our study. Its instantaneous value in the x-direction can be expressed as [42]

$$ I(t) = u(t) p(t). \qquad (12) $$

Using the formula $p = \rho c_{lbm}^2$, the intensity can then be computed in the LBM units as:

$$ I(t) = c_{lbm}^2 \sum_{i=0}^{18} f_i c_{x,i}. \qquad (13) $$

We will be particularly interested by the average intensity, which will be used to calculate the acoustic streaming force.
## 2.2 Generation of acoustic waves in the LB method

The literature gives two main methods to generate acoustic waves in LBM simulations. The waves can first be produced at the level of the boundary conditions, as proposed by Bouzidi et al. [43]. For a transducer considered to vibrate at a certain angular frequency $\omega$ and with a given velocity amplitude $u_a$, the idea is to write

$$f_i = f_i - 2 W_i c_i c_{bm}^{-2} u_a \sin(\omega t) \quad (14)$$

(where $f_i$ is the distribution function in the opposite direction ($\bar{i} = -i$)) at the boundary points on the transducer surface.

The second possibility is to use the acoustic point source method [36,44], which is mainly based on the approximation of a linear wave (weak oscillations). Pressure, density, or velocity are considered to vibrate around their equilibrium position as a sinusoidal function [36,42]:

$$p = p_0 + p_a \sin(\omega t), \quad (15)$$

$$\rho = \rho_0 + \rho_a \sin(\omega t), \quad (16)$$

$$u = u_0 + u_a \sin(\omega t). \quad (17)$$

where $\rho_0$, $p_0$ and $u_0$ are the density, pressure, and velocity ($u_0 = 0$) at steady state, respectively.

In the case of a single-point source, the application of the condition on the density through Eq. (16) is appropriate to generate the waves, as the waves will travel isotropically in all directions, i.e. with 3D spherical symmetry. In contrast, for a vibrating surface, our previous work [36] has shown that the application of the fluctuation condition on the velocity gives acoustic results closer to the analytical data. In this case, the use of Eq. (17) is appropriate. This condition will then be used in the present study to model the waves generated by a circular acoustic source.

## 2.3 Implementation of a body force in the LB method

In the LB method, different models can be used to implement a body force such as the acoustic streaming force [45]. The most popular models reported in the literature are that of Shan and Chen [46,47] and that suggested by Luo [48]. In the first model, the force is added to the momentum defined with the distribution function, so that we have $j$ expressed as

$$j(x, y, z) = \sum f_i c_i + \frac{\partial t F_{ac}}{2}, \quad (18)$$

the coefficient $1/2$ being chosen so that the Chapman-Enskog expansion used to get the Navier-Stokes equations from the Boltzmann equation will give the Navier-Stokes equations with the right force $F$, in this case [49].

In the second model, the external force is discretized and added to the lattice Boltzmann equation (Eq. (1)). The force is then expressed as [25,35]

$$F_i = W_i c_i F_{ac} c_{bm}^{-2}, \quad (19)$$
and Eq. (1) becomes [37]

$$f_i(x_i + c_i \delta t, t + \delta t) - f_i(x_i, t) = M^{-1}S[m_i^{eq} - m_i] + \delta t W_i c_i F_{ac} c_{bm}^{-2}. \quad (20)$$

The second model proposed by Luo is used to implement the acoustic streaming force into our LBM numerical code.

## 2.4 Boundary conditions

In the case of the lattice Boltzmann method, the boundary conditions are not based on velocity or pressure for example, but rather on the distribution function. In this work, only two types of conditions will be discussed. It concerns the bounce-back (BBC) and non-reflecting (NRBC) boundary conditions.

The bounce-back boundary conditions are mainly used to model stationary solid walls (homogeneous Dirichlet boundary conditions). As the name implies, when a particle reaches the boundary, it bounces in the opposite direction to its arrival. No flow crosses the wall (impermeable wall). Mathematically, these conditions can be implemented as [25–27]

$$f_i(x_B, t) = f_i(x_B, t), \quad (21)$$

where $f_i(x_B, t)$ is the unknown distribution function at the boundary node ($x_B$), which comes from outside into the study domain, and $f_i(x_B, t)$ is the known distribution function in the opposite direction ($\bar{i} = -i$), which comes from the simulation domain towards the outside.

The BBC conditions are the most commonly used for inputs and outputs in LBM simulations. However, these conditions are far from ideal for flows containing acoustic waves, as they will reflect the acoustic waves into the study domain. In our work, we consider situations where an absorber is disposed at the wall opposite to the source in order to avoid the production of standing waves in the cavity. Non-reflecting boundary conditions (NRBC) are then needed and there are two main approaches in LBM simulations [26]. The first approach is based on characteristic boundary conditions (CBC) [50–52]: the fluid equations are decomposed at the boundary nodes in order to cancel the travel velocity of the reflected waves. The second approach is based on absorbing layers (AL) [53,54] which are added beyond the boundary where reflection is to be avoided.

Both approaches have their drawback. The use of CBC conditions requires the solution of some mathematical differential equations at the boundary and the use of AL conditions implies the increase of the grid size, and these changes can become costly in terms of computational time, especially when it concerns 3D simulations. For the purpose of simplicity, other boundary conditions are used in this work to attenuate the waves. These are the modified bounce-back conditions, where the emitted wave arriving at the non-reflecting wall is subtracted from the distribution function at this wall [30]:

$$f_i = f_i - 3 W_i c_i u_a \sin(\omega t - kL_x) e^{-\alpha L_x}, \quad (22)$$

where $L_x$ (the length of the cavity along $x$) is the distance between the source and the absorbing wall. The mathematical expression of the coefficient $\alpha$ is given in the next section. This non-reflecting bounce-back condition was already used in one of our previous papers [36] and shown to be very efficient.

## 2.5 Link between the physical and LBM units
Simulations by LBM lead to results at the mesoscopic scale and expressed in LBM units. It is therefore necessary to clarify the passage from the LBM units to the physical units at the macroscopic scale. The best-known methods in the literature for unit conversion processes between mesoscopic and macroscopic scales are the Buckingham π theorem and the principle of corresponding states [55,56]. The Buckingham π theorem generally consists in mapping dimensionless parameters such as the Reynolds number, Prandtl number, Rayleigh number, and Mach number, to geometric similarities, such as the aspect ratio [56]. Some examples of the application of this theorem in LBM simulations can be found in the references [55,57]. In some cases, the use of the Buckingham π theorem to perform unit conversion becomes difficult, in particular in nonlinear problems studies [55]. Concerning the principle of corresponding states, it consists in exploiting dimensionless reduced numbers to establish proportional relations between physical and LBM units. Finally, more recently, Baakeem et al. [56] propose a dimensional analysis approach to systematically make conversions.

The different conversion techniques mentioned above can be applied to study different physical problems, such as fluid flows and heat and mass transfers. However, for acoustic problems, which involve wave propagation, other conversion techniques have been proposed. The most well-known procedure [36,58–60] is that which allows direct relationships between the LBM and physical acoustic quantities by means of the space step $\Delta x$ (m) and the time step $\Delta t$ (s), which represent the physical distance between the LBM lattice nodes and the physical time needed by a particle to pass from one node to another, respectively.

The physical quantities, most commonly used in acoustics problems, are the speed of sound, the frequency or period of the wave, the wavelength, and the viscosity, responsible for the wave attenuation. For each quantity, it is possible to find a mathematical expression involving $\Delta x$ and $\Delta t$, which can connect their values in the two-unit spaces. For example, the physical speed of sound and viscosity can be related to their LBM counterparts by:

$$c_{ph} = \frac{\Delta x}{\Delta t} c_{lbm}, \qquad (23)$$

$$v_{ph} = \frac{\Delta x^2}{\Delta t} v_{lbm}. \qquad (24)$$

The $ph$ and $lbm$ indices are used to refer to the quantities calculated in the physical scale and the LBM units, respectively.

The quantities $c_{ph}$, $c_{lbm}$, $v_{ph}$ expressed in Eqs. (23) and (24) are known. For air and the D3Q19 model, their values are given in Tables 1 and 2. Concerning the viscosity $v_{lbm}$, the choice of its value generally depends on the problem studied and the LBM model used. It is therefore adjustable. For the SRT model, viscosities higher than 0.1 must be used to stabilize the calculations [25]. In contrast, for the MRT model, it was shown in [36] that small values of $v_{lbm}$ can be used without affecting the stability of the calculations. Once the value of $v_{lbm}$ is chosen, Eqs. (23) and (24) can be used to obtain $\Delta x$ and $\Delta t$:

$$\Delta x = \frac{v_{ph}}{v_{lbm}} \frac{c_{lbm}}{c_{ph}}, \qquad (25)$$

$$\Delta t = \frac{v_{ph}}{v_{lbm}} \left( \frac{c_{lbm}}{c_{ph}} \right)^2 . \qquad (26)$$
Using the values of $\Delta x$ and $\Delta t$, the other quantities such as the physical period ($T_{ph}$) and wavelength ($\lambda_{ph}$) can be easily calculated from those used in the LBM simulations:

$$T_{ph} = \Delta t T_{lbm} \quad \text{and} \quad \lambda_{ph} = \Delta x \lambda_{lbm}. \qquad (27)$$

In the case of the simulation of acoustic streaming flow, additional acoustic quantities are used, such as the acoustic intensity and the acoustic force. The use of equation (12) and the formula $p = \rho c^2$ first enables to express the acoustic intensity as $I_{ac} = \rho u c^2$. In physical and LBM units, this intensity is given by $I_{ac,ph} = \rho_{ph} u_{ph} c_{ph}^2$ and $I_{ac,lbm} = \rho_{lbm} u_{lbm} c_{lbm}^2$. Dividing the first equation by the second and using $u_{ph}/u_{lbm} = c_{ph}/c_{lbm} = \Delta x/\Delta t$ allows us to express $I_{ac,ph}$ as a function of $I_{ac,lbm}$:

$$I_{ac,ph} = \frac{\rho_{ph}}{\rho_{lbm}} \left(\frac{\Delta x}{\Delta t}\right)^3 I_{ac,lbm}. \qquad (28)$$

The acoustic force can then be defined from the intensity $I_{ac}$ as [5,61]:

$$F_{ac} = 2\alpha \frac{l_{ac}}{c}. \qquad (29)$$

Using (28) and $\alpha_{ph}/\alpha_{lbm} = 1/\Delta x$, we can write:

$$F_{ac,ph} = \frac{\rho_{ph}}{\rho_{lbm}} \frac{\Delta x}{\Delta t^2} F_{ac,lbm}. \qquad (30)$$

The conversion described above makes it possible to relate directly the physical quantities ($\phi_{ph}$) to the LBM quantities ($\phi_{lbm}$) in the form $\phi_{ph} = a \phi_{lbm}$, where the parameter $a$ is the conversion factor which generally depends on $\Delta x$ and $\Delta t$. Table 3 lists the various physical quantities used in acoustics in relation to those obtained by LBM numerical simulations.

Details on the determination of the parameters used in this work will be given in section 4.

<table><caption>**Table 1.** LBM parameters used in the simulations of acoustic waves propagation.</caption><thead><tr><th>Δx<sub>lbm</sub></th><th>Δt<sub>lbm</sub></th><th>T<sub>lbm</sub></th><th>λ<sub>lbm</sub></th><th>c<sub>lbm</sub></th><th>ν<sub>lbm</sub></th><th>ρ<sub>lbm</sub></th><th>u<sub>a</sub></th></tr></thead><tbody><tr><td>1</td><td>1</td><td>34.64</td><td>20</td><td>1/√3</td><td>0.015</td><td>1</td><td>0.01 – 0.1</td></tr></tbody></table>

<table><caption>**Table 2.** Physical characteristics associated with the LBM simulations of acoustic waves propagation in air and properties of air considered (at 20°C).</caption><thead><tr><th colspan="4">Physical characteristics</th><th colspan="3">Properties of air</th></tr><tr><th>Δx (μm)</th><th>Δt(ns)</th><th>f<sub>ph</sub> (MHz)</th><th>λ<sub>ph</sub> (μm)</th><th>c<sub>ph</sub> (m/s)</th><th>ν<sub>ph</sub> (m<sup>2</sup>/s)</th><th>ρ<sub>ph</sub> (kg/m<sup>3</sup>)</th></tr></thead><tbody><tr><td><strong>1.698</strong></td><td>2.883</td><td>10.011</td><td>33.961</td><td>340</td><td>1.5 10<sup>-5</sup></td><td>1.204</td></tr></tbody></table>
<table><caption>**Table 3.** Factors *a* used to convert acoustic quantities from the LBM scale to the physical scale ($\phi_{ph} = a \phi_{lbm}$).</caption><thead><tr><th>Physical quantity</th><th>SI unit</th><th>Dimensions</th><th>Conversion factor from LBM to physical scale</th></tr></thead><tbody><tr><td>Distance</td><td>m</td><td>L</td><td>Δx</td></tr><tr><td>Time</td><td>s</td><td>t</td><td>Δt</td></tr><tr><td>Velocity</td><td>m/s</td><td>Lt<sup>-1</sup></td><td>Δx/Δt</td></tr><tr><td>Attenuation coefficient</td><td>1/m</td><td>L<sup>-1</sup></td><td>1/Δx</td></tr><tr><td>Mass</td><td>kg</td><td>M'</td><td>ρ<sub>ph</sub>/ρ<sub>lbm</sub> Δx<sup>3</sup></td></tr><tr><td>Pressure</td><td>Pa</td><td>M'L<sup>-1</sup> t<sup>-2</sup></td><td>ρ<sub>ph</sub>/ρ<sub>lbm</sub> (Δx/Δt)<sup>2</sup></td></tr><tr><td>Intensity</td><td>W/m<sup>2</sup></td><td>M't<sup>-3</sup></td><td>ρ<sub>ph</sub>/ρ<sub>lbm</sub> (Δx/Δt)<sup>3</sup></td></tr><tr><td>Force per unit volume</td><td>N/m<sup>3</sup></td><td>M'L<sup>-2</sup> t<sup>-2</sup></td><td>ρ<sub>ph</sub>/ρ<sub>lbm</sub> Δx/Δt<sup>2</sup></td></tr><tr><td>Power</td><td>W</td><td>M'L<sup>2</sup> t<sup>-3</sup></td><td>ρ<sub>ph</sub>/ρ<sub>lbm</sub> Δx<sup>5</sup>/Δt<sup>3</sup></td></tr></tbody></table>

### 3. Results and Discussion

The acoustic results obtained by the LBM simulations are presented in this section. We first describe the physical problem, then study the propagation of the waves, and finally calculate the acoustic streaming. We also give the advantages and limitations of the LBM approach for the study of acoustic waves and streaming.

#### 3.1 Description of the physical problem and definition of the acoustic parameters

We consider a parallelepiped cavity filled with air (see Fig.2). A circular acoustic source is positioned at the center of the left vertical wall ($x = 0$) in order to produce acoustic waves. The spatial dimensions of the cavity are $L_x$, $L_y$, and $L_z$ corresponding to length, width, and height, respectively ($L_y = L_z$). The size of the acoustic source is characterized by its diameter which is defined by $d = L_z/3$. The aspect ratio of the cavity is $L_x/L_z = 4/3$. The right vertical wall opposite the source is generally assumed to be absorbing, whereas the other walls are considered reflective, with a usual bounce-back boundary condition.

To ensure that the waves are well represented in the LBM grid, it is necessary to choose LBM wavelengths that are defined over a fairly large number of grid nodes. For example, to deal with acoustic streaming in 2D, Haydock and Yeomans [30] chose a wavelength defined by 50 nodes ($\lambda_{lbm} = 50$). They also used a large size mesh to discretize their 2D domain (2000 × 400 intervals). This choice is possible in 2D studies, but rather smaller size meshes are used in 3D studies, due to the limited capabilities of the computers. In our previous work on 3D wave simulation [36], we chose $\lambda_{lbm} = 23.12$ and a mesh size of 320 × 240 × 240 intervals, and the computations were performed on short times (small number of iterations). In this work, however, as we will see later, it is necessary to use a
very large number of iterations to directly simulate acoustic streaming. Thus, we decided to keep meshes
of similar size and to use a wavelength $\lambda_{lbm} = 20$. This choice leads to a LBM period $T_{lbm} =$
$\lambda_{lbm}/c_{lbm} = 34.64$.

To find the physical characteristics associated with our simulations in air, we used the direct conversion
described in section 2. The properties of air, sound speed $c_{ph}$, kinematic viscosity $\nu_{ph}$, density $\rho_{ph}$, are
given in Table 2. The physical values of $\Delta x$ and $\Delta t$, which can be obtained from Eqs. (23) and (24), then
only depend on the LBM viscosity. If we note that the ratio $c_{lbm}/c_{ph}$ is small and the kinematic
viscosity $\nu_{ph}$ is a small quantity (Tables 1 and 2), the only way to get a not too small $\Delta x$ from Eq. (23)
is to take the smallest possible value of $\nu_{lbm}$. In this study, we chose to fix $\nu_{lbm}$ at 0.015, a value
allowed by our MRT model, but not too small to avoid instability problems. This value is used to define
the relaxation times $s_9, s_{11}, s_{13}, s_{14}$ and $s_{15}$, as expressed in Eq. (10). This choice will be discussed and
justified in more detail below. Using $\nu_{lbm} = 0.015$, we get $\Delta x = 1.698 \ \mu m$ and $\Delta t = 2.883 \ ns$, and
Eq. (27) then gives a physical frequency of 10.011 MHz and a wavelength $\lambda_{ph} = 33.961 \ \mu m$.

The choice of the acoustic wave amplitude is constrained by two limitations of the LB method.
The first limitation is associated with the compressibility of the fluid, quantified by the Mach
number $Ma = u_a/c_{lbm}$, which must remain weak. It is a consequence of the simplified
isoentropic equation of state $p = \rho c^2$ used by the model. The LBM approach indeed is more
stable for simulations in incompressible or weakly compressible mediums [25,26]. Therefore,
for stability reasons, it is generally recommended to use values of $Ma < 0.3$ [62,63], which
gives values of LBM velocity amplitude $u_a$ lower than 0.173. The second limitation is the linear
wave assumption of the acoustic model described in Eqs. (15-17). Indeed, this model is only
valid for cases of weak oscillations, i.e. oscillations with an amplitude that is small compared
to the equilibrium value of the forced quantity. For example, in the case of density used to
produce the waves (Eq. (16)), Viggen [58] recommends a value of $\rho_a$ less than twice $\rho_0$, i.e.
$\rho_a < 2$ in LBM units. In our study where the waves are generated by the fluctuation of the
velocity around its equilibrium value (Eq. (17)), using $u_a = \rho_a c_{lbm}$ [36], we can express the
constraint in terms of $u_a$: $u_a < 2 c_{lbm} \approx 1.15$. The velocity amplitudes chosen in this study,
$u_a = 0.01$ and 0.1, fulfill well both these constraints. The corresponding physical values for
the proposed cases in air can be determined from the reference velocity $V_{ref} = \Delta x/\Delta t \approx$
589 m/s. We obtain 5.89 and 58.9 m/s for $u_a = 0.01$ and 0.1, respectively.

It should be noted that all the quantities that are presented in the figures in the following are
expressed in LBM units. What can be then obtained in terms of physical quantities will depend
on the fluid considered and its properties. In practice, these properties will allow to prescribe
the values of $\Delta x$ (25) and $\Delta t$ (26) (as proposed for air in Table 2) and then to obtain the
conversion factors (described in subsection 2.5) allowing to switch to the physical units for the
different quantities. The indications on the physical quantities will be given in the text, generally
in the case of air, sometimes for water, as in sections 3.3 and 3.4.3.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Geometry of the physical problem studied.</figcaption></figure>

## 3.2 Waves generated by a circular acoustic source in a cavity

In the cavity filled with air presented in Fig. 2, we first study the propagation of the waves generated by a circular acoustic source at the left vertical wall. We consider two different meshes with 240 × 180 × 180 and 560 × 420 × 420 nodes. If the physical grid space is fixed at the value of Δx presented in Table 2, the two meshes will correspond to cavities and sound sources of different sizes. This will give two test cases with the following characteristics in air:

**Case 1:** $L_{x,ph} = 407.541 \mu m = 12 \lambda_{ph}$, $L_{y,ph} = L_{z,ph} = 305.656 \mu m = 9 \lambda_{ph}$, $d_{ph} = \frac{L_{z,ph}}{3} = 101.885 \mu m$.

**Case 2:** $L_{x,ph} = 950.92 \mu m = 28 \lambda_{ph}$, $L_{y,ph} = L_{z,ph} = 713.19 \mu m = 21 \lambda_{ph}$, $d_{ph} = \frac{L_{z,ph}}{3} = 237.73 \mu m$.

In the two test cases, the acoustic source vibrates with a period $T_{lbm} = 34.64$, which gives a frequency of 10.011 MHz in air. The amplitude $u_a$ of the emitted acoustic velocity wave is fixed at 0.01 and the wall opposite the source is considered absorbing.

Fig. 3 shows the instantaneous acoustic velocity field obtained in case 1 after 1400 iterations, which corresponds to a physical time $t_{ph} = 1400 \Delta t = 4.037 \mu s$ for the propagation in air. Fig. 3(A) gives a 3D view, whereas Fig. 3(B) gives a vertical cross-section at $y = L_y/2$. Each point source of the discretized source (with 61 points along the diameter) emits spherical waves. The resulting waves in the cavity take the shape of discs, rather flat, with curvatures on the sides, in the near field, and slightly curved in the far field. The 3D view in Fig. 3(A) shows that the velocity field is not perfectly axisymmetric around the central x-axis (the symmetry axis of the source). This can be explained by the discretization of the source which is not axisymmetric and by the square section of the cavity. As about 420 iterations are necessary for the wave to reach the end wall, opposite to the source, the wave can be considered as already well-established in the cavity. As expected, we observe 12 oscillations along the horizontal direction. These waves are quite intense close to the source, but their intensity strongly
decreases when they move towards the end wall. The absorbing condition applied at this wall is efficient, as no reflection of the wave is observed: the wave remains a traveling wave, without being affected by interferences. The reflection at the side walls also appears very weak, justifying that absorbing conditions were not needed at these walls.

The instantaneous acoustic velocity field obtained in case 2 after 2000 iterations ($t_{ph} = 5.767\mu s$ in air) is presented in Fig. 4. The cavity is larger in this case and the circular source has 141 points along the diameter. The waves travel towards the end wall, but, due to the longer propagation distance along $x$ ($L_x = 28\lambda$), they are strongly attenuated in the right end zone. Moreover, the change of the source diameter changes the characteristics of the beam, in particular the transverse position of the maximum and minimum acoustic velocities.

Fig. 5 shows the longitudinal velocity profiles plotted along the axis of the transducer (longitudinal centerline of the cavity) for both cases 1 and 2. We see that the waves fluctuate differently in the two cases. This is first due to the change in source diameter, which also implies a change in the Fresnel length (generally given by $L_F = d^2/(4\lambda)$), the limit between the acoustic near field and far field: $L_{F,lbm} = 45$ in case 1 and $L_{F,lbm} = 245$ in case 2. For case 1, Fig. 5(A) shows that the amplitude of the wave starts with a value of 0.005, reaches a maximum value of 0.0062 at the third fluctuation, at about the Fresnel length, and then decreases exponentially with distance. The wave presents 12 oscillations before reaching the wall opposite the source where it is absorbed. In case 2, the wave amplitude evolution along the 28 successive oscillations is more complex. After the initial amplitude of 0.0041, the wave keeps a rather strong amplitude during 4 oscillations, then strongly decreases, before a last increase and the final exponential decrease in the far field zone. Note that the last amplitude maximum occurs clearly before the Fresnel length in this case. From these results, it can be noted that the wave attenuation is very large in both test cases. This is explained by the high value of the vibration frequency because the attenuation coefficient is proportional to the square of this frequency.

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Longitudinal acoustic velocity field $u$ produced after 1400 iterations (4.037$\mu$s in the air) by a circular acoustic source vibrating with a period $T_{lbm} = 34.64$ (10.011 MHz frequency in the air) with a LBM amplitude $u_a = 0.01$ in the microcavity defined in case 1: (A) 3D representation: the 3D yellow and green surfaces correspond to the smallest absolute isovalues in the scale, i.e. 0.0005 and</figcaption></figure>
-0.0005, respectively. An isovalue is also shown in the plane at x = 0; (B) 2D view in the vertical section at y = L_y/2.

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Longitudinal acoustic velocity field *u* produced after 2000 iterations (5.767 µs in the air) by a circular acoustic source vibrating with a period *T*<sub>lbm</sub> = 34.64 (10.011 MHz frequency in the air) with an amplitude *u*<sub>a</sub> = 0.01 in the microcavity defined in case 2: (A) 3D representation: the 3D yellow and green surfaces correspond to the smallest absolute isovvalues in the scale, i.e. 0.0005 and -0.0005, respectively. An isovalue is also shown in the plane at *x* = 0; (B) 2D view in the vertical section at *y* = *L*<sub>y</sub>/2.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Profiles along the central *x*-axis for the longitudinal acoustic velocity *u* produced by a circular source vibrating with a period *T*<sub>lbm</sub> = 34.64 (10.011 MHz frequency in the air) with an amplitude *u*<sub>a</sub> = 0.01 in two microcavities of different sizes: (A) case 1, as in Fig. 3; (B) case 2, as in Fig. 4.</figcaption></figure>

To validate our numerical calculations, we propose to compute the acoustic intensity and compare it
with analytical results. As the direction of wave propagation is along the *x*-axis, only the acoustic
intensity along this direction is calculated. At an instant $t$, it is defined as the product of the velocity $u(t)$ and the acoustic pressure $p(t)$ (Eq. (12)) and calculated with the LBM approach according to Eq. (13). The average acoustic intensity $I_{ac} = \langle I(t) \rangle$ is then obtained by averaging $I(t)$ over a large number of iterations, 400 in case 1 (iteration 1000 to 1400) and 600 in case 2 (iteration 1400 to 2000). For reasons of symmetry and better visualization, this intensity is plotted in Fig. 6 with a 3D view in only half of the cavity. In the two cases, the intensity structure is complex, with minima and maxima alternating in both longitudinal and transverse directions ahead of the acoustic source. The complexity is increased in case 2, when the diameter of the source is larger.

<figure><img src="image_7.png" /><figcaption>Fig. 6. 3D half view of the acoustic intensity field $I_{ac}$ produced by a circular source vibrating with a period $T_{lbm} = 34.64$ (10.011 MHz frequency in the air) with an amplitude $u_a = 0.01$. Intensity is averaged over 400 iterations in case 1 (A) and 600 iterations in case 2 (B). The 3D blue surface corresponds to the smallest isovalue in the scale, i.e. $2 \times 10^{-7}$ in (A) and $2 \times 10^{-8}$ in (B). In both figures, the isovalues given in the scale are also shown in the central plane at $y = L_y/2$.</figcaption></figure>

To validate the acoustic intensity presented in Fig. 6, we will compare it to the analytical intensity defined on the central $x$-axis of a circular acoustic source. We first recall that, without attenuation, the acoustic pressure $p(x)$ on this central axis [36,42] can be expressed as:

$$p(x) = p_{max} \left| \sin \left( \frac{1}{2} kx \left( \sqrt{1 + \left(\frac{d}{2x}\right)^2} - 1 \right) \right) \right|. \quad (31)$$

In (31), $p_{max}$ is the maximum pressure without attenuation, $x$ refers to the distance along the centerline of the transducer, and $k$ is the wave vector. $p_{max}$ can be expressed in LBM units as a function of the
equilibrium density $\rho_0$, sound speed $c_{lbm}$, and velocity amplitude $u_a$. The general expression, given for a baffled acoustic source, is $p_{max} = 2\rho_0 c_{lbm} u_a$. Factor 2 takes into account the fact that, with a baffled source, the radiation is restricted to a half-space, ahead of the source, and is then doubled [64]. In a previous work [36], we have shown that in our case where the radiation comes from point sources on the boundary, the amplitude of the waves is not doubled and, therefore, factor 2 must be removed. $p_{max}$ is then given by:

$$p_{max} = \rho_0 c_{lbm} u_a. \qquad (32)$$

To find an expression of the pressure in the case of attenuation, it is better to adapt the mathematical expression for the complex and unattenuated acoustic pressure radiated by a baffled piston on its central x-axis given in [42] to account for acoustic losses. We obtain

$$\bar{p} = -\frac{j\pi p_{max}}{\lambda (\alpha + jk)} e^{j(\omega t)} \left[ e^{-(\alpha + jk)\sqrt{x^2 + \frac{d^2}{4}}} - e^{-(\alpha + jk)x} \right]. \qquad (33)$$

Concerning the acoustic intensity, it is obtained in the general case as

$$I_{ac} = \langle p(t)u(t) \rangle = \mathrm{Re}\{\bar{p} \times \bar{u}^*\}/2, \qquad (34)$$

or in the often-used plane wave approximation [42,61] as

$$I_{ac} = \frac{\langle p(t)^2 \rangle}{\rho_0 c_{lbm}} = \frac{\mathrm{Re}\{\bar{p} \times \bar{p}^*\}}{2\rho_0 c_{lbm}}, \qquad (35)$$

where superscript * refers to the complex conjugate operator and $Re$ denotes the real part. To get the acoustic velocity involved in (34), we use the linear Euler equation [42]:

$$\rho_0 \frac{\partial u(t)}{\partial t} = - \frac{\partial p(t)}{\partial x}. \qquad (36)$$

Using the complex notation for signals varying in time as $e^{j(\omega t)}$, the velocity $\bar{u}$ can then be written as

$$\bar{u} = \frac{j}{\rho_0 \omega} \frac{\partial \bar{p}}{\partial x}, \qquad (37)$$

and the intensity given by (34) can then be expressed as a function of the complex pressure as

$$I_{ac} = \frac{1}{2\rho_0 \omega} \mathrm{Re} \left\{ \bar{p} \times \left( j \frac{\partial \bar{p}}{\partial x} \right)^* \right\}. \qquad (38)$$

The expression (33) of $\bar{p}$ can be introduced in either Eq. (38) or Eq. (35) to get the acoustic intensity along the x-axis. We obtain the general expression
$$ I_{ac} = I_m \left[ \frac{\alpha}{k} \left( \frac{x}{\sqrt{x^2 + \frac{d^2}{4}}} - 1 \right) e^{-\alpha \left( x + \sqrt{x^2 + \frac{d^2}{4}} \right)} \sin \left( k \left( \sqrt{x^2 + \frac{d^2}{4}} - x \right) \right) \right. \\ \left. - \left( \frac{x}{\sqrt{x^2 + \frac{d^2}{4}}} + 1 \right) e^{-\alpha \left( x + \sqrt{x^2 + \frac{d^2}{4}} \right)} \cos \left( k \left( \sqrt{x^2 + \frac{d^2}{4}} - x \right) \right) \right. \\ \left. + e^{-2\alpha x} + \frac{x}{\sqrt{x^2 + \frac{d^2}{4}}} e^{-2\alpha \sqrt{x^2 + \frac{d^2}{4}}} \right] \quad (39) $$

and a simpler expression in the plane wave approximation

$$ I_{ac} = I_m \left( e^{-2\alpha x} + e^{-2\alpha\sqrt{x^2+\frac{d^2}{4}}} - 2e^{-\alpha\left(x+\sqrt{x^2+\frac{d^2}{4}}\right)} \cos\left(k\left(\sqrt{x^2+\frac{d^2}{4}}-x\right)\right) \right) \quad (40) $$

with

$$ I_m = \frac{p_{max}^2}{8\rho_0 c_{lbm}} \quad (41) $$

Note that these expressions use the hypothesis that the attenuation length, related to $\alpha$, is large compared with the wavelength $\lambda$, more precisely that $\alpha\lambda/\pi \ll 1$.

To calculate $I_{ac}$ given by Eqs. (39) and (40), an estimation of the acoustic attenuation coefficient $\alpha$ is needed. Typically, the attenuation of an acoustic wave is due to two physical effects: geometric spreading and dissipation [60,65]. The first effect results from the fact that, as the sound moves away from the source, due to the spreading of the acoustic beam by diffraction, the area covered by the acoustic energy increases, and thus the intensity of the sound decreases. The second effect, related to dissipation, can be explained by three mechanisms: the viscous effects, the thermal effects, and the effects associated with internal molecular processes (rotational and vibrational relaxation effects of the molecules constituting the fluid) [8,42,65]. In the case of waves in gases such as air, the different dissipation effects are taken into account through the following attenuation coefficient [66,67]:

$$ \alpha = \frac{\omega^2}{2\rho_0 c_{lbm}^3} \left( \frac{4}{3}\mu + \frac{(\gamma - 1)^2 K}{\gamma R_{gas}} + \mu_B \right) \quad (42) $$

where $\gamma$ is the ratio of specific heats, $K$ is the thermal conductivity, $R_{gas}$ is the specific gas constant, and $\mu_B$ is the bulk dynamic viscosity.

The estimation of the two last terms that appear in Eq. (42) remains difficult. The bulk viscosity is sometimes obtained from experimental measurements, but they are few. It is often expressed as a function of the viscosity of the fluid used. For the case of water, different studies take this bulk viscosity as $\mu_B = 3\mu$ [5,42]. Concerning thermal effects, they are typically neglected when simulating wave propagation in water [5], as they are considered small compared to the viscous effects. In contrast, in
the case of liquid metals, these effects have an important contribution to the attenuation of waves and cannot be neglected.

In the present study, where the LBM method is used to simulate ultrasound propagation in air, attenuation is involved through the different relaxation times. A discussion on these relaxation times and their effects on wave propagation was proposed in our previous work [36]. Here, in our calculations, the value $v_{lbm} = 0.015$ (related to $v_{ph}$ of air by Eqs. (23) and (24)) is introduced in the relaxation times $s_9, s_{11}, s_{13}, s_{14}$ and $s_{15}$, as mentioned in the previous section. The thermal effects are not taken into account in our LBM method and then cannot be considered in our calculations. Finally, the bulk kinematic viscosity ($\nu_B = \mu_B/\rho$) is involved in the relaxation time $s_1$ defined as $s_1 = 2/(9\nu_{B,lbm} + 1)$. To stabilize the MRT-D3Q19 model, it is recommended to set $s_1$ to 1.19 [38]. Using this value leads to a bulk kinematic viscosity $\nu_{B,lbm} = 7.563 \times 10^{-2}$.

Using Eq. (42) in LBM units with $v_{lbm} = 0.015$ and $\nu_{B,lbm} = 0.07563$ and without thermal effects, we obtain $\alpha_{lbm} = 8.173 \times 10^{-3}$. It is this value of the attenuation coefficient which is used in the comparisons of our numerical results with the analytical expressions (40) and (39) (for such value, $\alpha\lambda/\pi \approx 0.05$). These comparisons are shown in Fig. 7 with the plots of the longitudinal profiles of the acoustic intensity along the central x-axis, which is the axis of the source. In case 1 (Fig. 7(A)), the acoustic intensity obtained numerically starts with a peak value of $1.25 \times 10^{-5}$, drops to about 0 and then increases towards the last peak with a maximum value of about $1.26 \times 10^{-5}$. Beyond this peak, the intensity rapidly decreases with distance in the far-field zone. In case 2 (Fig. 7(B)), the situation is quite different as three peaks are now found in the near field zone before the last peak and the final decrease in the far field zone. The intensity starts with a peak value of about $1.1 \times 10^{-5}$ and has then three peaks with decreasing amplitudes, the last peak being rather weak. The comparison with the analytical data shows that, in the two cases, the different peaks found numerically are well predicted by the analytical expressions, with practically the same position. The amplitude of the peaks, however, is differently predicted by the two analytical expressions, expression (40) involving the plane wave approximation giving larger amplitudes. The interesting result is that it is the more accurate analytical expression (39) which compares rather well with the numerical intensity profiles, at the peaks and more globally all along the x-axis. The comparison is particularly striking in case 2 with the four intensity peaks. These comparisons show that our LBM code is well adapted to study the propagation of ultrasound acoustic waves and that the plane wave approximation is not appropriate for accurately describing the waves in the near field zone.

Some comments can be done on the intensity profiles shown in Fig. 7. As we will discuss the number of peaks and their positions (and not their amplitudes), we will refer to the plane wave approximation. We have first quoted the different number of peaks between case 1 and case 2. In fact, according to the simple expression without attenuation (31), the number of peaks in the near field zone is given by $N = \text{Ent}(\frac{d}{2\lambda} - \frac{1}{2})$, with Ent the entire part, and there is also the last peak which marks the transition to the far field zone. In LBM units, the wavelength $\lambda$ is 20 and the source diameter $d$ is 60 in case 1 and 140 in case 2. We then obtain $N = 1$ in case 1 and $N = 3$ in case 2, which corresponds well to what is shown in Fig. 7, the first peak being exactly at $x = 0$ in the two cases. Concerning the last peak, its position without attenuation is usually given by the Fresnel length $L_F = d^2/(4\lambda)$, or $L_F = \frac{d^2}{4\lambda} - \frac{\lambda}{4}$ if the wavelength is not assumed small compared to the source diameter. These two expressions give $L_F = 45$ and 40 in case 1 and $L_F = 245$ and 240 in case 2, values which do not fit well with the results shown in Fig. 7, particularly in case 2. In fact, we have to take into account the acoustic attenuation which will
decrease the intensity of the peaks, but also displace them towards smaller $x$. This effect will particularly affect the last peak. Its position in case of attenuation cannot be expressed analytically, but can be calculated, from the intensity expression (40) for example, by a Newton method. In this case, we get $L_F = 34.65$ and $170.85$ in cases 1 and 2, respectively. These values now fit rather well with the position of the last peaks in Fig. 7.

It is interesting to calculate the physical value of the attenuation coefficient to see if it is appropriate for ultrasound waves in the air. This value $\alpha_{ph}$ can be obtained as $\alpha_{ph} = \alpha_{lbm}/\Delta x$. With a space step $\Delta x = 1.698 \ \mu m$ obtained in our case with $v_{ph}$ of air and $v_{lbm} = 0.015$ (Tables 1 and 2), we get $\alpha_{ph} \approx 4800 \ m^{-1}$. This value of the acoustic attenuation coefficient seems quite large, but it is of the same order of magnitude as the theoretical values that can be extrapolated from the results of Lin et al. [67] for $f_{ph} = 10$ MHz. Finally, despite the neglect of the thermal effects in our LBM approach, the choice of a bulk viscosity defined from $\eta_1 = 1.19$ allows to get a reasonable estimate of the acoustic attenuation in the air in our microcavities.

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Profiles of acoustic intensity $I_{ac}$ along the central $x$-axis of the cavity for a circular source vibrating with a period $T_{lbm} = 34.64$ (10.011 MHz frequency in the air) with an amplitude $u_a = 0.01$: case 1 (A) and case 2 (B). The numerical profiles are compared with the analytical expressions (40) involving plane wave approximation and (39) derived with more general hypotheses. The usual position of the last peak at the Fresnel length $L_F = d^2/(4\lambda)$ (orange vertical line) is compared to its position taking into account attenuation and obtained from (40) (black vertical line).</figcaption></figure>

### 3.3 Discussion on the choice of the LBM viscosity

As shown before, the MRT model allows flexibility in the choice of the LBM viscosity. In the present calculations, we have chosen a value $v_{lbm} = 0.015$, but it could be interesting to discuss the possible choice of values greater or smaller.

For a larger value of the LBM viscosity, the space step $\Delta x$ and time step $\Delta t$ will decrease, in the same way, according to Eqs. (25) and (26). For example, the choice of $v_{lbm} = 0.033$, the value used by
Haydock and Yeomans [30], gives $\Delta x = 0.772 \ \mu m$ and $\Delta t = 1.31 \ ns$. For a given LBM period, it would also decrease the physical period and increase the physical frequency, giving a frequency $f_{ph} \approx 22$ MHz. Using $v_{lbm} = 0.1$ will still increase these effects leading, for example, to $f_{ph} \approx 66.7$ MHz and $\Delta x = 0.254 \ \mu m$. Thus, the use of large values of LBM viscosity leads to simulations figuring very high frequency waves in very small geometries. This is particularly interesting because the consideration of high-frequency ultrasound (frequencies ranging from 50 MHz to a few GHz) in microgeometries is a new topic in the literature, connected with important applications of the acoustic streaming flows generated by such ultrasounds. Among these applications, we can find intracellular delivery [68], micromixing [69], and microscale particle [70] and nanoscale particle [71] manipulation.

Decreasing the LBM viscosity to some extent is also a possibility for the MRT model. For example, using $v_{lbm} = 0.001$ gives a relatively larger space step ($\Delta x = 25.47 \ \mu m$) and a fairly low frequency ($f_{ph} = 0.66$ MHz) compared to the previous cases. Thus, decreasing the viscosity would allow to deal with waves at a smaller frequency in geometries of increasing size. However, the treatment of configurations with a decimetric size as in the case of reference [61] remains difficult. For example, for a cavity with a characteristic length of 25 cm and with the space step just given, about $10^6$ nodes would be required in each space direction, which is beyond the present computing possibilities.

Similar difficulties were found by Salomons et al. [60] for the simulation of sound propagation in air. In their simulation with 2000 nodes in the propagation direction, if, as they propose, a lattice spacing $\Delta x = 0.1 \ m$ is chosen, the time step is $\Delta t = 1.7 \ 10^{-4} \ s$, and, for a period $T_{lbm} = 40$, the wave frequency is then $f_{ph} = 147.2 \ Hz$. Such simulation is performed with $v_{lbm} = 0.06$ (limitations due to the SRT model used) and corresponds then to $v_{ph} = 3.5 \ m^2/s$, whereas the real value of the viscosity for air is about $1.5 \ 10^{-5} \ m^2/s$. In this case, the interpretation of the LBM simulation as a real size experiment leads to a non-realistic value of the physical viscosity and so to a bad prediction of the wave attenuation and the streaming induced by this attenuation. This confirms the difficulties for the LBM approach to accurately simulate the propagation of low and medium-frequency waves in realistic fluids.

The present calculations with $v_{lbm} = 0.015$ could also be interpreted as an experiment in water. The choice of a bulk kinematic viscosity $\nu_{B,lbm} = 7.563 \ 10^{-2}$ is a little too high, as previous studies rather take this bulk viscosity as $\mu_B = 3\mu$ [5,42] and thermal effects are negligible [5], so that wave attenuation will be a little increased. Nevertheless, using the values $c_{ph} = 1480 \ m/s$, $\nu_{ph} = 10^{-6} \ m^2/s$ for water, we get $\Delta x \approx 26 \ nm$, $\Delta t \approx 10^{-11} \ s$, and, for a wavelength $\lambda_{lbm} = 20$, the wave frequency is then $f_{ph} = 2.845 \ GHz$. It will then correspond to the propagation of a gigahertz ultrasound wave in a micrometric cavity, a configuration typically in the microgeometries issue just mentioned.

## 3.4 Acoustic streaming produced by a circular acoustic source in a cavity

After validation of the wave calculation, the LBM method is now used to calculate streaming. As they could be very costly in terms of computing time, the streaming calculations are performed in case 1, using the smallest mesh with 240 × 180 × 180 nodes. The characteristics are then those given previously in case 1, with in particular the acoustic period $T_{lbm} = 34.64$. Two different methods will be used. The direct method is based on the long-time evolution of the waves until a stabilized oscillatory state is reached, and the streaming velocities (denoted with capital letters, e.g. $U$ in the $x$ direction) are then obtained by averaging the instantaneous acoustic velocities over a large number of acoustic periods. The indirect method is based on the calculation of the acoustic streaming force from the acoustic intensity $I_{ac}$, which, according to (13) and (29), can be written as:
$$F_{ac} = 2\alpha c_{lbm} \left\langle \sum_{i=0}^{18} f_i c_{x,i} \right\rangle, \qquad (43)$$

the average being taken over 400 iterations, as indicated for $I_{ac}$ in case 1 in the previous subsection. The force is then introduced in the LBM code to produce the streaming.

Note that the main objective here is the direct calculation of the streaming, in a way that mimics the real development of the streaming flow from the absorbed acoustic waves. The indirect calculation is rather used for comparison. The force introduced in this indirect calculation is here obtained from the wave propagation. In many studies, the body force is obtained as the Rayleigh integral summing the emitted contributions of the discretized acoustic source, without solving acoustic propagation [72].

Two cases will be considered. The first case corresponds to a rather small wave amplitude $u_a = 0.01$ and the second case to a larger amplitude $u_a = 0.1$. As the time step used is the same in both direct and indirect methods, the CPU times are rather similar for both methods. They depend on the number of iterations performed and correspond to almost five days of computations for 100 000 iterations in the cases we have performed.

### 3.4.1 Vibration amplitude $u_a = 0.01$

#### a) Direct method

In order to estimate the time evolution of the waves necessary to get a stabilized oscillatory state for $u_a = 0.01$, we plot the time variation of the instantaneous acoustic velocity $u$ in two different points along the central x-axis in Fig. 8. For the point at mid-length ($x = L_x/2$, black curve), we observe that the oscillations, first centered on zero, progressively evolve around a non-zero value, which eventually stabilizes, giving the streaming velocity at this point. From this plot, we see that 80 000 to 100 000 iterations are needed to have a stabilized oscillatory state. This time is very long compared to the oscillation period which corresponds to only about 35 iterations (see the zoomed view). For the second point at $x = 9L_x/10$ (blue curve), the long-time variation is weaker, indicating a smaller streaming velocity at this point. From these observations, we decided to calculate the streaming velocities by averaging the instantaneous velocities over 40 000 iterations, from iteration 100 000 to 140 000.
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Time variation of the acoustic velocity *u* at different positions along the central *x*-axis for a circular acoustic source vibrating with an amplitude *u*<sub>a</sub> = 0.01: *x* = L<sub>x</sub>/2 (black) and *x* = 9L<sub>x</sub>/10 (blue).</figcaption></figure>

The streaming velocity field thus obtained is presented by different 2D views in Figs. 9 and 10. The longitudinal component *U* of the streaming is plotted in the longitudinal plane at *y* = *L*<sub>y</sub>/2 and the transverse plane at *x* = *L*<sub>x</sub>/2 in Fig. 9, whereas velocity vector plots and streamlines in the plane at *y* = *L*<sub>y</sub>/2 are shown in Fig. 10. We see that, for *u*<sub>a</sub> = 0.01, the strongest velocities are located in a kind of jet close to the source. They determine a main toroidal roll centered in the first half of the cavity. Streamlines can penetrate farther inside the cavity, but they are associated with smaller velocities and will correspond to slow trajectories. Note that, due to the square transverse section of the cavity, the return flow occurs preferentially in the four corners, where confinement is less than in the center of the edges, but the return flow velocities remain weak compared to the jet velocities.
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** 2D views of the streaming velocity *U* generated by a circular acoustic source vibrating with an amplitude *u*<sub>a</sub> = 0.01: (A) longitudinal section at *y* = *L*<sub>y</sub>/2 and (B) transverse section at *x* = *L*<sub>x</sub>/2.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Velocity vectors (A) and streamlines (B) in the central plane at *y* = *L*<sub>y</sub>/2 for the streaming velocity generated by a circular acoustic source vibrating with an amplitude *u*<sub>a</sub> = 0.01. The streamlines are used to depict the structure of the flow, but do not give information on its intensity.</figcaption></figure>

**b) Indirect method for validation**

The acoustic force $F_{ac}$, which is injected in the LBM code in this case, is obtained by averaging the
wave propagation from iteration 1000 to 1400. The 3D acoustic force field is presented in Fig. 11. The
force is maximum along the central x-axis, close to the source, and the intensity of the force strongly
decreases for larger x and outside the acoustic beam. The time variation of the streaming velocity $U$ in
two different points along the central x-axis is shown in Fig. 12. For the two points, the time to reach a
steady state is similar to the time necessary to get a stabilized oscillatory state in the direct method (Fig.
8). The streaming velocity field created by this force field inside the cavity is almost identical to the
streaming presented in Figs. 9 and 10, which was obtained as a result of the direct calculation. We then
only show some comparisons between the streaming fields obtained by the two methods through
characteristic profiles in Fig. 13. For the longitudinal profile of $U$ along the central x-axis as well as its
transverse profile along the central y-axis, the comparison is very good, with differences not exceeding
a few percent. Fig. 13(A) shows that the maximum velocity along the central x-axis is reached quite quickly, at only a quarter of the cavity. Fig. 13(B) points out the weakness of the return flow compared to the central flow.

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** 3D half view of the acoustic force field $F_{ac}$ produced by a circular source vibrating with an amplitude $u_a = 0.01$. The force is averaged over 400 iterations. The 3D blue surface corresponds to the smallest isovalue in the scale, i.e. $1 \times 10^{-8}$. The isovalues given in the scale are also shown in the central plane at $y = L_y/2$.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Time variation of the streaming velocity *U* obtained by the LBM code with an imposed acoustic streaming force calculated for a source velocity amplitude *u*<sub>a</sub> = 0.01. The force is depicted in Fig. 11. The streaming velocity is given at two positions along the central *x*-axis: *x* = L<sub>x</sub>/2 (black) and *x* = 9L<sub>x</sub>/10 (blue).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Longitudinal profile along the central *x*-axis (A) and transverse profile along the central *y*-axis (B) for the streaming velocity *U* generated by a circular acoustic source vibrating with an amplitude *u*<sub>a</sub> = 0.01. Comparison between the results obtained by an indirect method involving the calculation of the acoustic force and by a direct method based on the long-time evolution of the acoustic waves.</figcaption></figure>
### 3.4.2 Vibration amplitude $u_a = 0.1$

#### a) Direct method

For this larger velocity vibration amplitude $u_a = 0.1$, the time variation of the instantaneous acoustic velocity $u$ in two different points along the central $x$-axis is plotted in Fig. 14. We see that for both points (at $x = L_x/2$ and $x = 9L_x/10$), the time evolution of the waves necessary to get a stabilized oscillatory state is strongly reduced compared with the case at $u_a = 0.01$. Only 20 000 iterations are now needed, and this time is even less for regions with higher velocities as at $x = L_x/2$. In both cases, the oscillations first decrease in amplitude while almost centered on zero, and are then very suddenly shifted towards the streaming velocity around which the oscillation amplitude eventually settles. The final oscillation amplitude appears smaller than the initial amplitude. Finally, to get the streaming velocities, we chose to average the instantaneous velocities over 40 000 iterations, from iteration 60000 to 100000.

The streaming velocity field thus obtained is presented by different 2D views in the longitudinal plane at $y = L_y/2$ and the transverse plane at $x = L_x/2$ in Figs. 15 and 16. For $u_a = 0.1$, the flow in the cavity appears quite different from that obtained previously for a smaller oscillation amplitude $u_a = 0.01$. The longitudinal component $U$ of the streaming has now strong values up to the end wall, indicating the presence of a jet all along the $x$-axis (Fig. 15). The jet spreads on this end wall, creating transverse velocities, and then returns towards the left wall with decreasing longitudinal velocities (Fig. 16). The streaming flow in the cavity is a large toroidal roll with a center of rotation close to the end wall. Here also, due to the square transverse section of the cavity, the return flow occurs preferentially in the four corners.
<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Time variation of the acoustic velocity *u* at different positions along the central *x*-axis for a circular acoustic source vibrating with an amplitude *u*<sub>a</sub> = 0.1: *x* = L<sub>x</sub>/2 (black) and *x* = 9L<sub>x</sub>/10 (blue).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 15.** 2D views of the streaming velocity *U* generated by a circular acoustic source vibrating with an amplitude *u*<sub>a</sub> = 0.1: (A) longitudinal section at *y* = L<sub>y</sub>/2 and (B) transverse section at *x* = L<sub>x</sub>/2.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Fig. 16.** Velocity vectors (A) and streamlines (B) in the central plane at $y = L_y/2$ for the streaming velocity generated by a circular acoustic source vibrating with an amplitude $u_a = 0.1$. The streamlines are used to depict the structure of the flow, but do not give information on its intensity.</figcaption></figure>

b) Indirect method for validation

The acoustic force $F_{ac}$, which is injected in the LBM code in this case, is the same as that calculated for $u_a = 0.01$ and presented in Fig. 11. Only its intensity has to be multiplied by 100 (see Fig. 17). The time variation of the streaming velocity $U$ in two different points along the central x-axis is shown in Fig. 18. Here also, for the two points, the time to reach a steady state is similar to the time necessary to get a stabilized oscillatory state in the direct method (Fig. 14). The streaming velocity field created by the force field inside the cavity is almost identical to the streaming presented in Figs. 15 and 16 and obtained as a result of the direct calculation. The comparisons between the streaming fields obtained by the two methods are shown in Fig. 19. For $u_a = 0.1$, the comparisons are still very good, with differences not exceeding a few percent, except in the near-source region. In this region, the indirect method gives positive longitudinal velocities, increasing from 0 at the source wall, which is coherent with the positive applied body force and the bounce-back boundary condition at the wall. In contrast, the direct method, in which the velocities are obtained as the mean acoustic velocities, gives negative longitudinal velocities in this near-source region. The reason for this behavior is not clear to us. Note that it occurs in the complex near-field zone of the acoustic field, and that, farther inside the cavity, the two calculations give very similar results. A similar (but weaker) effect was in fact already observed in Fig. 13 for $u_a = 0.01$.
<figure><img src="image_18.png" /><figcaption>**Fig. 17.** 3D half view of the acoustic force field $F_{ac}$ produced by a circular source vibrating with an amplitude $u_a = 0.1$. The force is averaged over 400 iterations. The 3D blue surface corresponds to the smallest isovalue in the scale, i.e. $3 \times 10^{-7}$. The isovalues given in the scale are also shown in the central plane at $y = L_y/2$.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Fig. 18.** Time variation of the streaming velocity *U* obtained by the LBM code with an imposed acoustic streaming force calculated for a source velocity amplitude *u*<sub>a</sub> = 0.1. The force is depicted in Fig. 17. The streaming velocity is given at two positions along the central *x*-axis: *x* = L<sub>x</sub>/2 (black) and *x* = 9L<sub>x</sub>/10 (blue).</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. 19.** Longitudinal profile along the central *x*-axis (A) and transverse profile along the central *y*-axis (B) for the streaming velocity *U* generated by a circular acoustic source vibrating with an amplitude *u*<sub>a</sub> = 0.1. Comparison between the results obtained by an indirect method involving the calculation of the acoustic force and by a direct method based on the long-time evolution of the acoustic waves.</figcaption></figure>
### 3.4.3 Discussion

The simulations presented in the previous subsections gave the streaming flow structure obtained for two values of the velocity oscillation amplitude $u_a$. They also indicate the intensity of the streaming flow, which, according to Figs. 9 and 15, has maximum values around $u_{lbm,max} = 1.5 \ 10^{-3}$ and $4 \ 10^{-2}$ for $u_a = 0.01$ and 0.1, respectively. Using the reference velocity given by $c_{ph}/c_{lbm}$, indications on the physical streaming velocities thus created can be obtained.

If we consider the experiment in air, $c_{ph}/c_{lbm} \approx 589 \ m/s$, and the maximum physical velocities will be $u_{ph,max} = 0.88$ and $3.84 \ m/s$ for $u_a = 0.01$ and 0.1, respectively. In the case of an experiment in water, $c_{ph}/c_{lbm} \approx 2563 \ m/s$, and the maximum physical velocities will then be $u_{ph,max} = 23.56$ and $102.5 \ m/s$ for $u_a = 0.01$ and 0.1, respectively. We then see that the experiments in water with high wave frequencies and in micrometric geometries will generate stronger streaming flows than the experiments in the air with smaller wave frequencies and in, still small, but larger geometries.

## 4. Conclusion

This paper was devoted to the simulation of ultrasound waves propagation and the study of the acoustic streaming flow produced by the dissipation of the acoustic energy of these waves during their propagation in a parallelepiped cavity filled with air. The waves are generated by a circular acoustic source, which vibrates at about 10 MHz. The numerical approach used is the lattice Boltzmann method (LBM) based on the D3Q19 multiple relaxation time model. The ability of this model to simulate wave propagation was shown by comparing the numerical and analytical acoustic intensities along the central axis of the acoustic source. A very good comparison is obtained when the analytical model does not use the plane wave approximation and involves the right viscosities responsible for the attenuation. In particular, the simulations were thus shown to be able to reproduce the different peaks of the acoustic intensity profiles, with the right positions and amplitudes.

The LBM technique has then been used to study the acoustic streaming induced by the waves in the 3D cavity. The acoustic streaming was simulated using two different methods. The first method is based on the direct computation of the mean velocity fields obtained as the mean values of the instantaneous acoustic velocities. The second method is an indirect technique, which has been used in previous simulations with conventional CFD methods, as in [72]. In this method, the acoustic force is first calculated (analytically or, as here, by a wave propagation simulation) and is then injected into the numerical code to produce the acoustic streaming flow. Both methods need a very large number of iterations due to the very small-time step, necessary for wave propagation simulation, compared with the relatively long time needed for acoustic streaming establishment. A comparison between the results obtained by the two methods was carried out and a good agreement was found between them. The acoustic streaming flows thus obtained correspond to a toroidal roll, with the fluid moving in the wave propagation direction around the transducer axis and returning along the side walls, preferentially in the corners. Depending on the velocity vibration amplitude $u_a$, the toroidal roll is rather developed close to the transducer (small $u_a$) or in the whole cavity, with the main central jet spreading on the end wall (larger $u_a$).

These different investigations have allowed us to discuss the advantages and limitations of the LBM approach to simulate real situations of wave propagation and acoustic streaming. The LBM simulations are performed in a non-dimensional approach and they have then to be interpreted in real dimensions to see the concrete situations they represent. If our interest in the wave propagation simulations is the good
estimate of the streaming flow intensities, it is important to consider the real values of the fluid
viscosities that are involved in the streaming generation. The fact that these values are small for usual
fluids such as air or water and that $v_{lm}$ cannot be chosen very small leads to small values of the space
and time steps $\Delta x$ and $\Delta t$ (Eqs. (25-26)). As a consequence, the LBM is mainly able to simulate the
propagation of ultrasound and the streaming generation at very high frequencies (MHz-GHz) and in
microgeometries. The simulations in geometries of relatively large sizes (cm to dm) would require the
discretization of the domain into millions of points, and this would exceed the capacity of the best
computers, particularly in 3D situations. Similarly, the consideration of low and medium-frequency
waves is difficult as this would require many time steps in a period and many space steps in a
wavelength. The progress of LBM to simulate wave propagation and acoustic streaming seems then to
be strongly connected with the improvement of computer capacities.

Finally, it should be noted that the limitations due to the use of relatively high LBM viscosity values are associated with standard LBM methods such as the MRT and SRT models mentioned in the article. However, there are other sophisticated LBM models in which very low viscosity values can be achieved with reasonable stability [73]. This could be interesting for the simulation of acoustic streaming by the lattice Boltzmann method in larger-size cavities and could be the subject of further works.

## Acknowledgments

This collaborative work was supported by the PHC Maghreb Partnership Program no. 36951NG and the Phase 2 PHC Toubkal Program. A grant to J.B. and the support from the P2CHPD of the University Claude Bernard Lyon 1 and the PMCS2I of Ecole Centrale de Lyon for the numerical calculations are also gratefully acknowledged.

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.

## Declarations

### Ethical Approval

This article does not contain any experimental studies on humans and/or animals.

### Competing interests

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Authors' contributions

**Jaouad Benhamou:** Conceptualization, Investigation, Methodology, Formal analysis, Software, Validation, Writing – original draft. **Bjarne Vincent:** Conceptualization, Investigation, Formal analysis, Validation, Writing – original draft. **Sophie Miralles:** Visualization, Investigation, Formal analysis, Writing – review & editing. **Mohammed Jami:** Supervision, Visualization, Methodology, Investigation, Formal analysis, Writing – review & editing. **Daniel Henry:** Conceptualization, Supervision, Visualization, Methodology, Investigation, Formal analysis, Writing – review & editing. **Ahmed Mezrhab:** Supervision, Investigation, Visualization, Methodology, Formal analysis, Writing –
review & editing. **Valéry Botton:** Supervision, Investigation, Visualization, Formal analysis, Writing –
review & editing.

Funding

No special funding was received for this study.

Availability of data and materials

No data or specific materials were used for the research described in the article.

References

[1] D. Ensminger, L.J. Bond, Ultrasonics: Fundamentals, technologies, and applications, Third Edition, CRC Press, Boca Raton, 2011. https://doi.org/10.1201/b11173.

[2] A.K.W. Wood, C.M. Sehgal, A review of low-intensity ultrasound for cancer therapy, Ultrasound Med. Biol. 41 (2015) 905–928. https://doi.org/10.1016/j.ultrasmedbio.2014.11.019.

[3] A. Carovac, F. Smajlovic, D. Junuzovic, Application of ultrasound in medicine, Acta Inform. Medica. 19 (2011) 168. https://doi.org/10.5455/aim.2011.19.168-171.

[4] A.P. Sarvazyan, M.W. Urban, J.F. Greenleaf, Acoustic Waves in medical imaging and diagnostics, Ultrasound Med. Biol. 39 (2013) 1133–1146. https://doi.org/10.1016/j.ultrasmedbio.2013.02.006.

[5] B. Moudjed, Caractérisation expérimentale et théorique des écoulements entraînés par ultrasons. Perspectives d'utilisation dans les procédés de solidification du Silicium Photovoltaïque, Ph.D. thesis, INSA Lyon, France, 2013.

[6] M.C. Charrier-Mojtabi, X. Jacob, T. Dochy, A. Mojtabi, Species separation of a binary mixture under acoustic streaming, Eur. Phys. J. E. 42 (2019) 1–8. https://doi.org/10.1140/epje/i2019-11824-9.

[7] F.J. Fuchs, Ultrasonic cleaning and washing of surfaces, in: Power Ultrason. Appl. High-Intensity Ultrasound, Woodhead Publishing, 2015: pp. 577–609. https://doi.org/10.1016/B978-1-78242-028-6.00019-3.

[8] J. Lighthill, Acoustic streaming, J. Sound Vib. 61 (1978) 391–418. https://doi.org/10.1016/0022-460X(78)90388-7.

[9] C. Eckart, Vortices and streams caused by sound, Phys. Rev. 73 (1948) 68–76. https://doi.org/10.1103/PhysRev.73.68.

[10] L. Rayleigh, On the circulation of air observed in Kundt's tubes, and on some allied acoustical problems, Philos. Trans. R. Soc. London. 175 (1948) 1–21. www.jstor.org/stable/109434.

[11] N. Riley, Acoustic streaming, Theor. Comput. Fluid Dyn. 10 (1998) 349–356. https://doi.org/10.1007/s001620050068.

[12] Y. Fang, T. Yamamoto, S. Komarov, Cavitation and acoustic streaming generated by different sonotrode tips, Ultrason. Sonochem. 48 (2018) 79–87. https://doi.org/10.1016/j.ultsonch.2018.05.011.
[13] J. Lei, G. Zheng, Z. Yao, Z. Huang, Outer acoustic streaming flow driven by asymmetric acoustic resonances, Micromachines. 13 (2021) 65. https://doi.org/10.3390/MI13010065.

[14] X. Zhang, Y. Cui, F. Qu, H. Li, Acoustic streaming in water induced by an asymmetric dielectric-barrier-discharge plasma actuator at the initiation stage, Phys. Fluids. 34 (2022) 017113. https://doi.org/10.1063/5.0076231.

[15] P. Zhang, C. Chen, X. Su, J. Mai, Y. Gu, Z. Tian, H. Zhu, Z. Zhong, H. Fu, S. Yang, K. Chakrabarty, T.J. Huang, Acoustic streaming vortices enable contactless, digital control of droplets, Sci. Adv. 6 (2020) eaba0606. https://doi.org/10.1126/sciadv.aba0606.

[16] D.-S. Guo, X.-B. Li, H.-N. Zhang, F.-C. Li, P.-J. Ming, Numerical Study on the vapor bubble removal by acoustic streaming in downward-facing pool boiling, Int. J. Heat Mass Transf. 196 (2022) 123264. https://doi.org/10.2139/ssrn.4137986.

[17] J. Wu, Acoustic streaming and its applications, Fluids. 3 (2018) 108. https://doi.org/10.3390/fluids3040108.

[18] Pp.K. Das, V.R. Bhethanabotla, Extra stress-mediated acoustic streaming in a surface acoustic wave driven microchannel filled with second-order fluids, Phys. Rev. Fluids. 7 (2022) 074404. https://doi.org/10.1103/PhysRevFluids.7.074404.

[19] G.B. Lebon, G. Salloum-Abou-Jaoude, D. Eskin, I. Tzanakis, K. Pericleous, P. Jarry, Numerical modelling of acoustic streaming during the ultrasonic melt treatment of direct-chill (DC) casting, Ultrason. Sonochem. 54 (2019) 171–182. https://doi.org/10.1016/j.ultsonch.2019.02.002.

[20] G. Michel, C. Gissinger, Cooling by baroclinic acoustic streaming, Phys. Rev. Appl. 16 (2021) L051003. https://doi.org/10.1103/physrevapplied.16.1051003.

[21] F. Nocera, D. Eskin, B. Zej Baran, K. Machaj, Z. Malecha, K. Tomczuk, Numerical study of baroclinic acoustic streaming phenomenon for various flow parameters, Energies. 15 (2022) 854. https://doi.org/10.3390/en15030854.

[22] W. Qiu, J.T. Karlsen, H. Bruus, P. Augustsson, Experimental characterization of acoustic streaming in gradients of density and compressibility, Phys. Rev. Appl. 11 (2019) 024018. https://doi.org/10.1103/PhysRevApplied.11.024018.

[23] U. Farooq, X. Liu, W. Zhou, M. Hassan, L. Niu, L. Meng, Cell lysis induced by nanowire collision based on acoustic streaming using surface acoustic waves, Sensors Actuators, B Chem. 345 (2021) 130335. https://doi.org/10.1016/j.snb.2021.130335.

[24] S. He, W. Pang, X. Wu, Y. Yang, W. Li, H. Qi, K. Yang, X. Duan, Y. Wang, Bidirectional regulation of cell mechanical motion via a gold nanorods-acoustic streaming system, ACS Nano. 16 (2022) 8427–8439. https://doi.org/10.1021/acsnano.2c02980.

[25] A.A. Mohamad, Lattice Boltzmann method: Fundamentals and engineering applications with computer codes, Second Edition, Springer London, 2019. https://doi.org/10.1007/978-1-4471-7423-3.

[26] T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, S. Goncalo, E.M. Viggen, The lattice Boltzmann method: Principles and practice, Springer Cham, 2017. https://doi.org/10.1007/978-3-319-44649-3.

[27] Z. Guo, C. Shu, Lattice Boltzmann Method and its Applications in Engineering, World Scientific Publishing, 2013. https://doi.org/10.1142/8806.
[28] P. Stansell, C.A. Greated, Lattice gas automaton simulation of acoustic streaming in a two-dimensional pipe, Phys. Fluids. 9 (1997) 3288-3299. https://doi.org/10.1063/1.869443.

[29] D. Haydock, J.M. Yeomans, Lattice Boltzmann simulations of acoustic streaming, J. Phys. A. Math. Gen. 34 (2001) 5201. https://doi.org/10.1088/0305-4470/34/25/304.

[30] D. Haydock, J.M. Yeomans, Lattice Boltzmann simulations of attenuation-driven acoustic streaming, J. Phys. A. Math. Gen. 36 (2003) 5683-5694. https://doi.org/10.1088/0305-4470/36/20/322.

[31] Y. Rafat, K. Habibi, L. Mongeau, Direct numerical simulations of acoustic streaming in standing wave tubes using the lattice Boltzmann method, in: Proc. Meet. Acoust., Acoustical Society of America, 2013: p. 045006. https://doi.org/10.1121/1.4800937.

[32] M.K. Tan, L.Y. Yeo, Hybrid finite-difference/lattice Boltzmann simulations of microchannel and nanochannel acoustic streaming driven by surface acoustic waves, Phys. Rev. Fluids. 3 (2018) 044202. https://doi.org/10.1103/PhysRevFluids.3.044202.

[33] X. Li, X. Li, W. Cui, T. Ma, L. Lu, Q. Wang, Effect of acoustic streaming on heat transfer of porous composite phase change material by using lattice Boltzmann simulation, J. Heat Transfer. 143 (2021) 092701. https://doi.org/10.1115/1.4051506.

[34] J. Benhamou, M. Jami, A. Mezrhab, Direct modeling method of acoustic streaming using lattice Boltzmann approach, in: 4th Int. Conf. Networking, Inf. Syst. Secur., ACM, New York, NY, USA, 2021: p. 6. https://doi.org/10.1145/3454127.3458771.

[35] J. Benhamou, M. Jami, A. Mezrhab, D. Henry, V. Botton, Numerical simulation study of acoustic waves propagation and streaming using MRT-lattice Boltzmann method, Int. J. Comput. Methods Eng. Sci. Mech. (2022) 1-14. https://doi.org/10.1080/15502287.2022.2050844.

[36] J. Benhamou, S. Channouf, M. Jami, A. Mezrhab, D. Henry, V. Botton, Three-dimensional lattice Boltzmann model for acoustic waves emitted by a source, Int. J. Comut. Fluid Dyn. 35 (2021) 850-871. https://doi.org/10.1080/10618562.2021.2019226.

[37] J. Benhamou, E.B. Lahmer, M. Jami, Three-dimensional simulation of conjugate heat transfer using the hybrid lattice Boltzmann-finite difference method, Int. Commun. Heat Mass Transf. 139 (2022) 106486. https://doi.org/10.1016/j.icheatmasstransfer.2022.106486.

[38] D. D'Humières, I. Ginzburg, M. Krafczyk, P. Lallemand, L.S. Luo, Multiple-relaxation-time lattice Boltzmann models in three dimensions, Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 360 (2002) 437-451. https://doi.org/10.1098/rsta.2001.0955.

[39] J. Benhamou, M. Jami, Three-dimensional numerical study of heat transfer enhancement by sound waves using mesoscopic and macroscopic approaches, Heat Transf. 51 (2022) 3892-391. https://doi.org/10.1002/htj.22482.

[40] Z. Li, M. Yang, Y. Zhang, Lattice Boltzmann method simulation of 3-D natural convection with double MRT model, Int. J. Heat Mass Transf. 94 (2016) 222-238. https://doi.org/10.1016/j.ijheatmasstransfer.2015.11.042.

[41] K.N. Premnath, J. Abraham, Three-dimensional multi-relaxation time (MRT) lattice-Boltzmann models for multiphase flow, J. Comput. Phys. 224 (2007) 539-559. https://doi.org/10.1016/j.jcp.2006.10.023.

[42] L.E. Kinsler, A.R. Frey, A.B. Coppens, J.V. Sanders, Fundamentals of Acoustics, 4th Edition, John Wiley & Sons, New York, 2000.
[43] M. Bouzidi, M. Firdaouss, P. Lallemand, Momentum transfer of a Boltzmann-lattice fluid with boundaries, Phys. Fluids. 13 (2001) 3452-3459. https://doi.org/10.1063/1.1399290.

[44] J. Benhamou, M. Jami, A. Mezrhab, V. Botton, D. Henry, Numerical study of natural convection and acoustic waves using the lattice Boltzmann method, Heat Transf. 49 (2020) 3779-3796. https://doi.org/10.1002/htj.21800.

[45] G. V. Krivovichev, Stability analysis of body force action models used in the single-relaxation-time single-phase lattice Boltzmann method, Appl. Math. Comput. 348 (2019) 25-41. https://doi.org/10.1016/j.amc.2018.11.056.

[46] X. Shan, H. Chen, Lattice Boltzmann model for simulating flows with multiple phases and components, Phys. Rev. E. 47 (1993) 1815. https://doi.org/10.1103/PhysRevE.47.1815.

[47] X. Shan, H. Chen, Simulation of nonideal gases and liquid-gas phase transitions by the lattice Boltzmann equation, Phys. Rev. E. 49 (1994) 2941. https://doi.org/10.1103/PhysRevE.49.2941.

[48] L.S. Luo, Unified theory of lattice Boltzmann models for nonideal gases, Phys. Rev. Lett. 81 (1998) 1618. https://doi.org/10.1103/PhysRevLett.81.1618.

[49] Z. Guo, C. Zheng, B. Shi, Discrete lattice effects on the forcing term in the lattice Boltzmann method, Phys. Rev. E. 65 (2002) 046308. https://doi.org/10.1103/PhysRevE.65.046308.

[50] N. Jung, H.W. Seo, C.S. Yoo, Two-dimensional characteristic boundary conditions for open boundaries in the lattice Boltzmann methods, J. Comput. Phys. 302 (2015) 191-199. https://doi.org/10.1016/j.jcp.2015.08.044.

[51] D. Heubes, A. Bartel, M. Ehrhardt, Characteristic boundary conditions in the lattice Boltzmann method for fluid and gas dynamics, J. Comput. Appl. Math. 262 (2014) 51-61. https://doi.org/10.1016/j.cam.2013.09.019.

[52] G. Wissocq, N. Gourdain, O. Malaspinas, A. Eyssartier, Regularized characteristic boundary conditions for the Lattice-Boltzmann methods at high Reynolds number flows, J. Comput. Phys. 331 (2017) 1-18. https://doi.org/10.1016/j.jcp.2016.11.037.

[53] A. Najafi-Yazdi, L. Mongeau, An absorbing boundary condition for the lattice Boltzmann method based on the perfectly matched layer, Comput. Fluids. 68 (2012) 203-218. https://doi.org/10.1016/j.compfluid.2012.07.017.

[54] H. Xu, P. Sagaut, Analysis of the absorbing layers for the weakly-compressible lattice Boltzmann methods, J. Comput. Phys. 245 (2013) 14-42. https://doi.org/10.1016/j.jcp.2013.02.051.

[55] J. Huang, C. Bao, Z. Jiang, X. Zhang, A general approach of unit conversion system in lattice Boltzmann method and applications for convective heat transfer in tube banks, Int. J. Heat Mass Transf. 135 (2019) 873-884. https://doi.org/10.1016/j.ijheatmasstransfer.2019.02.014.

[56] S.S. Baakeem, S.A. Bawazeer, A.A. Mohamad, A novel approach of unit conversion in the lattice Boltzmann method, Appl. Sci. 11 (2021) 6386. https://doi.org/10.3390/app11146386.

[57] X.L. H. Huang, M. Sukop, Multiphase lattice Boltzmann methods: Theory and application, John Wiley & Sons, 2015. https://doi.org/10.1002/9781118971451.

[58] E.M. Viggen, The lattice Boltzmann method with applications in acoustics, Master's thesis, Norges teknisk-naturvitenskapelige universitet, Fakultet for naturvitenskap og teknologi, Institutt for fysikk, 2009.
[59] E.M. Viggen, The lattice Boltzmann method: Fundamentals and acoustics, Ph.D. thesis, Norwegian University of Science and Technology, Trondheim, Norway, 2014.

[60] E.M. Salomons, W.J.A. Lohman, H. Zhou, Simulation of sound waves using the lattice Boltzmann method for fluid flow: Benchmark cases for outdoor sound propagation, PLoS One. 11 (2016) e0147206. https://doi.org/10.1371/journal.pone.0147206.

[61] 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. https://doi.org/10.1063/1.4895518.

[62] J. Wang, D. Wang, P. Lallemand, L.S. Luo, Lattice Boltzmann simulations of thermal convective flows in two dimensions, Comput. Math. with Appl. 65 (2013) 262-286. https://doi.org/10.1016/j.camwa.2012.07.001.

[63] A.A. Mohamad, A. Kuzmin, A critical evaluation of force term in lattice Boltzmann method, natural convection problem, Int. J. Heat Mass Transf. 53 (2010) 990-996. https://doi.org/10.1016/j.ijheatmasstransfer.2009.11.014.

[64] D.T. Blackstock, Fundamentals of physical acoustics, John Wiley & Sons, New York, 2000.

[65] E. M. Salomons, Computational atmospheric acoustics, Springer Dordrecht, 2001. https://doi.org/10.1007/978-94-010-0660-6.

[66] X.D. Li, Z.M. Hu, Z.L. Jiang, Continuum perspective of bulk viscosity in compressible fluids, J. Fluid Mech. 812 (2017) 966-990. https://doi.org/10.1017/jfm.2016.834.

[67] J. Lin, C. Scalo, L. Hesselink, Bulk viscosity model for near-equilibrium acoustic wave attenuation, ArXiv. (2017) 1707.05876. http://arxiv.org/abs/1707.05876.

[68] S. Pan, T. Jeon, D.C. Luther, X. Duan, V.M. Rotello, Cytosolic delivery of functional proteins in vitro through tunable gigahertz acoustics, ACS Appl. Mater. Interfaces. 12 (2020) 15823-15829. https://doi.org/10.1021/acsami.9b211131.

[69] W. Cui, W. Pang, Y. Yang, T. Li, X. Duan, Theoretical and experimental characterizations of gigahertz acoustic streaming in microscale fluids, Nanotechnol. Precis. Eng. 2 (2019) 15-22. https://doi.org/10.1016/j.npe.2019.03.004.

[70] X. Guo, J. Gu, S. Guo, Z. Xu, C. Yang, S. Liu, L. Cheng, K. Huang, 3D object detection and tracking based on streaming data, in: Proc. - IEEE Int. Conf. Robot. Autom., 2020: pp. 8376-8382. https://doi.org/10.1109/ICRA40945.2020.9197183.

[71] X. Guo, M. Sun, Y. Yang, H. Xu, J. Liu, S. He, Y. Wang, L. Xu, W. Pang, X. Duan, Controllable cell deformation using acoustic streaming for membrane permeability modulation, Adv. Sci. 8 (2021) 2002489. https://doi.org/10.1002/advs.202002489.

[72] B. Moudjed, V. Botton, D. Henry, S. Millet, J.P. Garandet, H. Ben Hadid, Near-field acoustic streaming jet, Phys. Rev. E Phys. 91 (2015) 033011. https://doi.org/10.1103/PhysRevE.91.033011.

[73] F. Renard, G. Wissocq, J.F. Boussuge, P. Sagaut, A linear stability analysis of compressible hybrid lattice Boltzmann methods, J. Comput. Phys. 446 (2021) 110649. https://doi.org/10.1016/j.jcp.2021.110649.