Jaouad Benhamou
email: benhamou_jaouad1718@ump.ac.ma

Bjarne Vincent

Sophie Miralles

Mohamed Jami

Daniel Henry
email: daniel.henry@ec-lyon.fr

Ahmed Mezrhab

Valéry Botton

Mohammed Jami

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

Keywords: 3D simulation, acoustic streaming, lattice Boltzmann method, ultrasound

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.

Nomenclature

Greek symbols 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 [START_REF] Moudjed | Caractérisation expérimentale et théorique des écoulements entraînés par ultrasons[END_REF] or the separation of species in a mixture [START_REF] Charrier-Mojtabi | Species separation of a binary mixture under acoustic streaming[END_REF], their use in daily life for cleaning [START_REF] Fuchs | Ultrasonic cleaning and washing of surfaces[END_REF].

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 [START_REF] Lighthill | Acoustic streaming[END_REF]. In general, two types of acoustic streaming flows are considered, Eckart streaming [START_REF] Eckart | Vortices and streams caused by sound[END_REF] and Rayleigh streaming [START_REF] Rayleigh | On the circulation of air observed in Kundt's tubes, and on some allied acoustical problems[END_REF]. 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 [START_REF] Riley | Acoustic streaming[END_REF]. 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 [START_REF] Fang | Cavitation and acoustic streaming generated by different sonotrode tips[END_REF][START_REF] Lei | Outer acoustic streaming flow driven by asymmetric acoustic resonances[END_REF][START_REF] Zhang | Acoustic streaming in water induced by an asymmetric dielectricbarrier-discharge plasma actuator at the initiation stage[END_REF]. Acoustic streaming was also studied as a way to control droplets and bubbles without contact [START_REF] Zhang | Acoustic streaming vortices enable contactless, digital control of droplets[END_REF][START_REF] Guo | Numerical Study on the vapor bubble removal by acoustic streaming in downward-facing pool boiling[END_REF]. 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 [START_REF] Wu | Acoustic streaming and its applications[END_REF][START_REF] Pp | Extra stress-mediated acoustic streaming in a surface acoustic wave driven microchannel filled with second-order fluids[END_REF], and in industrial processes, such as direct-chill casting [START_REF] Lebon | Numerical modelling of acoustic streaming during the ultrasonic melt treatment of direct-chill (DC) casting[END_REF]. Acoustic streaming can also influence heat transfer and convection in heated cavities [START_REF] Michel | Cooling by baroclinic acoustic streaming[END_REF][START_REF] Nocera | Numerical study of baroclinic acoustic streaming phenomenon for various flow parameters[END_REF]. With the recent discovery that standing sound waves generate much higher velocity flows if the medium under study is inhomogeneous, Qiu et al. [START_REF] Qiu | Experimental characterization of acoustic streaming in gradients of density and compressibility[END_REF] 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 [START_REF] Farooq | Cell lysis induced by nanowire collision based on acoustic streaming using surface acoustic waves[END_REF][START_REF] He | Bidirectional regulation of cell mechanical motion via a gold nanorods-acoustic streaming system[END_REF].

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 [START_REF] Mohamad | Lattice Boltzmann method: Fundamentals and engineering applications with computer codes[END_REF][START_REF] Krüger | The lattice Boltzmann method: Principles and practice[END_REF][START_REF] Guo | Lattice Boltzmann Method and its Applications in Engineering[END_REF]. 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 Greated [START_REF] Stansell | Lattice gas automaton simulation of acoustic streaming in a twodimensional pipe[END_REF] 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 [START_REF] Haydock | Lattice Boltzmann simulations of acoustic streaming[END_REF] 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 [START_REF] Haydock | Lattice Boltzmann simulations of attenuation-driven acoustic streaming[END_REF]. Rafat et al. [START_REF] Rafat | Direct numerical simulations of acoustic streaming in standing wave tubes using the lattice Boltzmann method[END_REF] used LBM simulations to study the streaming in a simplified thermoacoustic refrigerator. More recently, Tan et al. [START_REF] Tan | Hybrid finite-difference/lattice Boltzmann simulations of microchannel and nanochannel acoustic streaming driven by surface acoustic waves[END_REF] 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. [START_REF] Li | Effect of acoustic streaming on heat transfer of porous composite phase change material by using lattice Boltzmann simulation[END_REF] 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 [START_REF] Benhamou | Direct modeling method of acoustic streaming using lattice Boltzmann approach[END_REF][START_REF] Benhamou | Numerical simulation study of acoustic waves propagation and streaming using MRT-lattice Boltzmann method[END_REF]. 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 [START_REF] Benhamou | Numerical simulation study of acoustic waves propagation and streaming using MRT-lattice Boltzmann method[END_REF][START_REF] Benhamou | Three-dimensional lattice Boltzmann model for acoustic waves emitted by a source[END_REF], 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.

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 [START_REF] Mohamad | Lattice Boltzmann method: Fundamentals and engineering applications with computer codes[END_REF][START_REF] Krüger | The lattice Boltzmann method: Principles and practice[END_REF]. For the study of acoustic waves, the LBM-MRT model presents several advantages over the two other models, first in terms of accuracy [START_REF] Mohamad | Lattice Boltzmann method: Fundamentals and engineering applications with computer codes[END_REF] and then because of the possible use of small viscosity values [START_REF] Benhamou | Numerical simulation study of acoustic waves propagation and streaming using MRT-lattice Boltzmann method[END_REF][START_REF] Benhamou | Three-dimensional lattice Boltzmann model for acoustic waves emitted by a source[END_REF]. 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.

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 [START_REF] Mohamad | Lattice Boltzmann method: Fundamentals and engineering applications with computer codes[END_REF][START_REF] Krüger | The lattice Boltzmann method: Principles and practice[END_REF][START_REF] Benhamou | Three-dimensional lattice Boltzmann model for acoustic waves emitted by a source[END_REF][START_REF] Benhamou | Three-dimensional simulation of conjugate heat transfer using the hybrid lattice Boltzmann-finite difference method[END_REF]:

𝑓 𝑖 (𝑥 𝑖 + 𝑐 𝑖 𝑡, 𝑡 + 𝑡) = 𝑓 𝑖 (𝑥 𝑖 , 𝑡) +  𝑖 . (1)
In the D3Q19 lattice, particles can take the following 19 LBM velocities 𝑐 𝑖 = ( 𝑐 𝑥,𝑖 𝑐 𝑦,𝑖 𝑐 𝑧,𝑖 ) [START_REF] Benhamou | Three-dimensional lattice Boltzmann model for acoustic waves emitted by a source[END_REF][START_REF] Humières | Multiple-relaxation-time lattice Boltzmann models in three dimensions[END_REF]:

𝑐 𝑖 = 𝑐 ( 0 0 0 1 0 0 -1 0 0 0 1 0 0 -1 0 0 0 1 0 0 -1 1 1 0 -1 1 0 1 -1 0 -1 -1 0 1 0 1 -1 0 1 1 0 -1 -1 0 -1 0 1 1 0 -1 1 0 1 -1 0 -1 -1 ) (2)
Fig. 1. Illustration of the D3Q19-LBM lattice.

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 [START_REF] Mohamad | Lattice Boltzmann method: Fundamentals and engineering applications with computer codes[END_REF][START_REF] Krüger | The lattice Boltzmann method: Principles and practice[END_REF]]

 𝑖 = 𝑀 -1 𝑆[𝑚 𝑖 𝑒𝑞 -𝑚 𝑖 ]. (3)
The vector 𝑚 of components 𝑚 𝑖 is composed of 19 different macroscopic quantities, which are well discussed in references [START_REF] Benhamou | Three-dimensional lattice Boltzmann model for acoustic waves emitted by a source[END_REF][START_REF] Humières | Multiple-relaxation-time lattice Boltzmann models in three dimensions[END_REF]. However, typically, researchers are only interested in calculating four quantities corresponding to the density 𝜌 and the momentums 𝑗 𝑥 , 𝑗 𝑦 , and 𝑗 𝑧 along the 𝑥, 𝑦, and 𝑧 axes.

𝑚 can be represented mathematically as

𝑚 = (𝜌, 𝑒, 𝑒
The matrix 𝑀 -1 can be obtained by calculating the inverse of the transformation matrix 𝑀, which is defined as [START_REF] Benhamou | Three-dimensional lattice Boltzmann model for acoustic waves emitted by a source[END_REF][START_REF] Humières | Multiple-relaxation-time lattice Boltzmann models in three dimensions[END_REF][START_REF] Benhamou | Three-dimensional numerical study of heat transfer enhancement by sound waves using mesoscopic and macroscopic approaches[END_REF]]

𝑀 = ( 1 -30 1 -11 -4 1 -4 0 0 0 0 2 -4 0 0 0 0 0 0 0 0 1 -11 -4 -1 4 0 0 0 0 2 -4 0 0 0 0 0 0 0 0 1 -11 -4 0 0 1 -4 0 0 -1 2 1 -2 0 0 0 0 0 0 1 -11 -4 0 0 -1 4 0 0 -1 2 1 -2 0 0 0 0 0 0 1 -11 -4 0 0 0 0 1 -4 -1 2 -1 2 0 0 0 0 0 0 1 -11 -4 0 0 0 0 -1 4 -1 2 -1 2 0 0 0 0 0 0 1 8 1 1 1 1 1 0 0 1 1 1 1 1 0 0 1 -1 0 1 8 1 -1 -1 1 1 0 0 1 1 1 1 -1 0 0 -1 -1 0 1 8 1 1 1 -1 -1 0 0 1 1 1 1 -1 0 0 1 1 0 1 8 1 -1 -1 -1 -1 0 0 1 1 1 1 1 0 0 -1 1 0 1 8 1 1 1 0 0 1 1 1 1 -1 -1 0 0 1 -1 0 1 1 8 1 -1 -1 0 0 1 1 1 1 -1 -1 0 0 -1 1 0 1 1 8 1 1 1 0 0 -1 -1 1 1 -1 -1 0 0 -1 -1 0 -1 1 8 1 -1 -1 0 0 -1 -1 1 1 -1 -1 0 0 1 1 0 -1 1 8 1 0 0 1 1 1 1 -2 -2 0 0 0 1 0 0 1 -1 1 8 1 0 0 -1 -1 1 1 -2 -2 0 0 0 -1 0 0 -1 -1 1 8 1 0 0 1 1 -1 -1 -2 -2 0 0 0 -1 0 0 1 1 1 8 1 0 0 -1 -1 -1 -1 -2 -2 0 0 0 1 0 0 -1 1 ) (5)
The equilibrium moment vector 𝑚 𝑒𝑞 is related to the equilibrium distribution function 𝑓 𝑒𝑞 by 𝑚 𝑒𝑞 = 𝑀 𝑓 𝑒𝑞 , with 𝑓 𝑒𝑞 expressed in its discrete form as [START_REF] Mohamad | Lattice Boltzmann method: Fundamentals and engineering applications with computer codes[END_REF]

𝑓 𝑖 𝑒𝑞 = 𝑊 𝑖  [1 + 1 𝑐 𝑙𝑏𝑚 2 𝑐 ⃗ 𝑖 . 𝑉 ⃗⃗ + 1 2𝑐 𝑙𝑏𝑚 4 (𝑐 ⃗ 𝑖 . 𝑉 ⃗⃗ ) 2 - 1 2𝑐 𝑙𝑏𝑚 2 |𝑉 ⃗⃗ | 2 ], with 𝑖 = 0, … , 18 (6)
where 𝑉 ⃗⃗ is the velocity vector (𝑉 ⃗⃗ = (𝑢 ⃗⃗, 𝑣 ⃗, 𝑤 ⃗⃗⃗)) and 𝑊 𝑖 represents the discretization weights, which are given as:

𝑊 𝑖 = { 1 3 𝑖 = 0 1 18 𝑖 = 1, … , 6 1 36 𝑖 = 7, … , 18 (7)
The nineteen elements of 𝑚 𝑖 𝑒𝑞 found after performing the calculations can be given in terms of density and momentums as follows [START_REF] Humières | Multiple-relaxation-time lattice Boltzmann models in three dimensions[END_REF][START_REF] Benhamou | Three-dimensional numerical study of heat transfer enhancement by sound waves using mesoscopic and macroscopic approaches[END_REF]:

𝑚 0 𝑒𝑞 =  𝑚 1 𝑒𝑞 = -11 + 19 𝑗 𝑥 2 + 𝑗 𝑦 2 + 𝑗 𝑧 2  0 𝑚 2 𝑒𝑞 = 3 - 11 2 𝑗 𝑥 2 + 𝑗 𝑦 2 + 𝑗 𝑧 2  0 𝑚 3 𝑒𝑞 = 𝑗 𝑥 𝑚
The choice of 𝑠 𝑖 values is based in particular on the stability of the D3Q19 model. The values used in this study are chosen according to the literature [START_REF] Benhamou | Three-dimensional lattice Boltzmann model for acoustic waves emitted by a source[END_REF][START_REF] Humières | Multiple-relaxation-time lattice Boltzmann models in three dimensions[END_REF][START_REF] Li | Lattice Boltzmann method simulation of 3-D natural convection with double MRT model[END_REF][START_REF] Premnath | Three-dimensional multi-relaxation time (MRT) lattice-Boltzmann models for multiphase flow[END_REF]: ).

{
The acoustic intensity is also an interesting quantity in our study. Its instantaneous value in the 𝑥 -direction can be expressed as [START_REF] Kinsler | Fundamentals of Acoustics[END_REF] 𝐼(𝑡) = 𝑢(𝑡) 𝑝(𝑡). (

) 12
Using the formula 𝑝 =  𝑐 𝑙𝑏𝑚 2 , the intensity can then be computed in the LBM units as:

𝐼(𝑡) = 𝑐 𝑙𝑏𝑚 2 ∑ 𝑓 𝑖 18 𝑖=0 𝑐 𝑥,𝑖 . (13)
We will be particularly interested by the average intensity, which will be used to calculate the acoustic streaming force.

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. [START_REF] Bouzidi | Momentum transfer of a Boltzmann-lattice fluid with boundaries[END_REF]. For a transducer considered to vibrate at a certain angular frequency  and with a given velocity amplitude 𝑢 𝑎 , the idea is to write

𝑓 𝑖̅ = 𝑓 𝑖 -2 𝑊 𝑖 𝑐 𝑖 𝑐 𝑙𝑏𝑚 -2 𝑢 𝑎 𝑠𝑖𝑛(𝑡) (14)
(where 𝑓 𝑖̅ is the distribution function in the opposite direction (𝑖̅ = -𝑖)) at the boundary points on the transducer surface.

The second possibility is to use the acoustic point source method [START_REF] Benhamou | Three-dimensional lattice Boltzmann model for acoustic waves emitted by a source[END_REF][START_REF] Benhamou | Numerical study of natural convection and acoustic waves using the lattice Boltzmann method[END_REF], 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 [START_REF] Benhamou | Three-dimensional lattice Boltzmann model for acoustic waves emitted by a source[END_REF][START_REF] Kinsler | Fundamentals of Acoustics[END_REF]:

𝑝 = 𝑝 0 + 𝑝 𝑎 𝑠𝑖𝑛(𝑡), (15)
 =  0 +  𝑎 𝑠𝑖𝑛(𝑡), (16)
𝑢 = 𝑢 0 + 𝑢 𝑎 𝑠𝑖𝑛(𝑡). (17)
where  0 , 𝑝 0 and 𝑢 0 are the density, pressure, and velocity (𝑢 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 [START_REF] Benhamou | Three-dimensional lattice Boltzmann model for acoustic waves emitted by a source[END_REF] 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.

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 [START_REF] Krivovichev | Stability analysis of body force action models used in the single-relaxationtime single-phase lattice Boltzmann method[END_REF]. The most popular models reported in the literature are that of Shan and Chen [START_REF] Shan | Lattice Boltzmann model for simulating flows with multiple phases and components[END_REF][START_REF] Shan | Simulation of nonideal gases and liquid-gas phase transitions by the lattice Boltzmann equation[END_REF] and that suggested by Luo [START_REF] Luo | Unified theory of lattice Boltzmann models for nonideal gases[END_REF]. In the first model, the force is added to the momentum defined with the distribution function, so that we have 𝑗 expressed as

𝑗(𝑥, 𝑦, 𝑧) = ∑ 𝑓 𝑖 𝑐 𝑖 + 𝑡 𝐹 𝑎𝑐 2 , (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 𝐹, in this case [START_REF] Guo | Discrete lattice effects on the forcing term in the lattice Boltzmann method[END_REF].

In the second model, the external force is discretized and added to the lattice Boltzmann equation (Eq.

(1)). The force is then expressed as [START_REF] Mohamad | Lattice Boltzmann method: Fundamentals and engineering applications with computer codes[END_REF][START_REF] Benhamou | Numerical simulation study of acoustic waves propagation and streaming using MRT-lattice Boltzmann method[END_REF]]

𝐹 𝑖 = 𝑊 𝑖 𝑐 𝑖 𝐹 𝑎𝑐 𝑐 𝑙𝑏𝑚 -2 , ( 19
)
and Eq. ( 1) becomes [START_REF] Benhamou | Three-dimensional simulation of conjugate heat transfer using the hybrid lattice Boltzmann-finite difference method[END_REF] 𝑓 𝑖 (𝑥 𝑖 + 𝑐 𝑖 𝑡, 𝑡 + 𝑡) -𝑓 𝑖 (𝑥 𝑖 , 𝑡)

= 𝑀 -1 𝑆[𝑚 𝑖 𝑒𝑞 -𝑚 𝑖 ] + 𝑡 𝑊 𝑖 𝑐 𝑖 𝐹 𝑎𝑐 𝑐 𝑙𝑏𝑚 -2 . ( 20
)
The second model proposed by Luo is used to implement the acoustic streaming force into our LBM numerical code.

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 [START_REF] Mohamad | Lattice Boltzmann method: Fundamentals and engineering applications with computer codes[END_REF][START_REF] Krüger | The lattice Boltzmann method: Principles and practice[END_REF][START_REF] Guo | Lattice Boltzmann Method and its Applications in Engineering[END_REF] 𝑓 𝑖 (𝑥 𝐵 , 𝑡) = 𝑓 𝑖̅ (𝑥 𝐵 , 𝑡),

where 𝑓 𝑖 (𝑥 𝐵 , 𝑡) is the unknown distribution function at the boundary node (𝑥 𝐵 ), which comes from outside into the study domain, and 𝑓 𝑖̅ (𝑥 𝐵 , 𝑡) is the known distribution function in the opposite direction (𝑖̅ = -𝑖), 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. Nonreflecting boundary conditions (NRBC) are then needed and there are two main approaches in LBM simulations [START_REF] Krüger | The lattice Boltzmann method: Principles and practice[END_REF]. The first approach is based on characteristic boundary conditions (CBC) [START_REF] Jung | Two-dimensional characteristic boundary conditions for open boundaries in the lattice Boltzmann methods[END_REF][START_REF] Heubes | Characteristic boundary conditions in the lattice Boltzmann method for fluid and gas dynamics[END_REF][START_REF] Wissocq | Regularized characteristic boundary conditions for the Lattice-Boltzmann methods at high Reynolds number flows[END_REF]: 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) [START_REF] Najafi-Yazdi | An absorbing boundary condition for the lattice Boltzmann method based on the perfectly matched layer[END_REF][START_REF] Xu | Analysis of the absorbing layers for the weakly-compressible lattice Boltzmann methods[END_REF] 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 [START_REF] Haydock | Lattice Boltzmann simulations of attenuation-driven acoustic streaming[END_REF]:

𝑓 𝑖̅ = 𝑓 𝑖 -3 𝑊 𝑖 𝑐 𝑖 𝑢 𝑎 𝑠𝑖𝑛(𝑡 -𝑘𝐿 𝑥 ) 𝑒 - 𝐿 𝑥 , (22)
where 𝐿 𝑥 (the length of the cavity along 𝑥) is the distance between the source and the absorbing wall.

The mathematical expression of the coefficient  is given in the next section. This non-reflecting bounce-back condition was already used in one of our previous papers [START_REF] Benhamou | Three-dimensional lattice Boltzmann model for acoustic waves emitted by a source[END_REF] and shown to be very efficient.

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 [START_REF] Huang | A general approach of unit conversion system in lattice Boltzmann method and applications for convective heat transfer in tube banks[END_REF][START_REF] Baakeem | A novel approach of unit conversion in the lattice Boltzmann method[END_REF]. 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 [START_REF] Baakeem | A novel approach of unit conversion in the lattice Boltzmann method[END_REF]. Some examples of the application of this theorem in LBM simulations can be found in the references [START_REF] Huang | A general approach of unit conversion system in lattice Boltzmann method and applications for convective heat transfer in tube banks[END_REF][START_REF] Huang | Multiphase lattice Boltzmann methods: Theory and application[END_REF]. In some cases, the use of the Buckingham  theorem to perform unit conversion becomes difficult, in particular in nonlinear problems studies [START_REF] Huang | A general approach of unit conversion system in lattice Boltzmann method and applications for convective heat transfer in tube banks[END_REF]. 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. [START_REF] Baakeem | A novel approach of unit conversion in the lattice Boltzmann method[END_REF] 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 [START_REF] Benhamou | Three-dimensional lattice Boltzmann model for acoustic waves emitted by a source[END_REF][START_REF] Viggen | The lattice Boltzmann method with applications in acoustics[END_REF][START_REF] Viggen | The lattice Boltzmann method: Fundamentals and acoustics[END_REF][START_REF] Salomons | Simulation of sound waves using the lattice Boltzmann method for fluid flow: Benchmark cases for outdoor sound propagation[END_REF] is that which allows direct relationships between the LBM and physical acoustic quantities by means of the space step ∆𝑥 (𝑚) and the time step ∆𝑡 (𝑠), 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 Δ𝑥 and Δ𝑡, 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:

𝑐 𝑝ℎ = ∆𝑥 ∆𝑡 𝑐 𝑙𝑏𝑚 , (23)
 𝑝ℎ = ∆𝑥 2 ∆𝑡  𝑙𝑏𝑚 . (24)
The 𝑝ℎ and 𝑙𝑏𝑚 indices are used to refer to the quantities calculated in the physical scale and the LBM units, respectively.

The quantities 𝑐 𝑝ℎ , 𝑐 𝑙𝑏𝑚,  𝑝ℎ expressed in Eqs. ( 23) and ( 24) are known. For air and the D3Q19 model, their values are given in Tables 1 and2. Concerning the viscosity  𝑙𝑏𝑚 , 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 [START_REF] Mohamad | Lattice Boltzmann method: Fundamentals and engineering applications with computer codes[END_REF]. In contrast, for the MRT model, it was shown in [START_REF] Benhamou | Three-dimensional lattice Boltzmann model for acoustic waves emitted by a source[END_REF] that small values of  𝑙𝑏𝑚 can be used without affecting the stability of the calculations. Once the value of  𝑙𝑏𝑚 is chosen, Eqs. ( 23) and ( 24) can be used to obtain ∆𝑥 and ∆𝑡:

∆𝑥 =  𝑝ℎ  𝑙𝑏𝑚 𝑐 𝑙𝑏𝑚 𝑐 𝑝ℎ , ( 25
) ∆𝑡 =  𝑝ℎ  𝑙𝑏𝑚 ( 𝑐 𝑙𝑏𝑚 𝑐 𝑝ℎ ) 2 . ( 26
)
Using the values of ∆𝑥 and ∆𝑡, the other quantities such as the physical period (𝑇 𝑝ℎ ) and wavelength (λ 𝑝ℎ ) can be easily calculated from those used in the LBM simulations:

𝑇 𝑝ℎ = ∆𝑡 𝑇 𝑙𝑏𝑚 and λ 𝑝ℎ = ∆𝑥 λ 𝑙𝑏𝑚 . ( 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 𝑝 = 𝑐 2 first enables to express the acoustic intensity as 𝐼 𝑎𝑐 = 𝑢𝑐 2 . In physical and LBM units, this intensity is given by 𝐼 𝑎𝑐,𝑝ℎ =  𝑝ℎ 𝑢 𝑝ℎ 𝑐 𝑝ℎ 2 and 𝐼 𝑎𝑐,𝑙𝑏𝑚 =  𝑙𝑏𝑚 𝑢 𝑙𝑏𝑚 𝑐 𝑙𝑏𝑚 2 . Dividing the first equation by the second and using 𝑢 𝑝ℎ /𝑢 𝑙𝑏𝑚 = 𝑐 𝑝ℎ /𝑐 𝑙𝑏𝑚 = ∆𝑥/∆𝑡 allows us to express 𝐼 𝑎𝑐,𝑝ℎ as a function of 𝐼 𝑎𝑐,𝑙𝑏𝑚 :

𝐼 𝑎𝑐,𝑝ℎ =  𝑝ℎ  𝑙𝑏𝑚 ( ∆𝑥 ∆𝑡 ) 3 𝐼 𝑎𝑐,𝑙𝑏𝑚 . ( 28
)
The acoustic force can then be defined from the intensity 𝐼 𝑎𝑐 as [START_REF] Moudjed | Caractérisation expérimentale et théorique des écoulements entraînés par ultrasons[END_REF][START_REF] Moudjed | Scaling and dimensional analysis of acoustic streaming jets[END_REF]:

𝐹 𝑎𝑐 = 2𝛼 𝐼 𝑎𝑐 𝑐 . ( 29
)
Using ( 28) and  𝑝ℎ / 𝑙𝑏𝑚 = 1/∆𝑥, we can write:

𝐹 𝑎𝑐,𝑝ℎ =  𝑝ℎ  𝑙𝑏𝑚 ∆𝑥 ∆𝑡 2 𝐹 𝑎𝑐,𝑙𝑏𝑚. (30)
The conversion described above makes it possible to relate directly the physical quantities ( 𝑝ℎ ) to the LBM quantities ( 𝑙𝑏𝑚 ) in the form  𝑝ℎ = 𝑎  𝑙𝑏𝑚 , where the parameter 𝑎 is the conversion factor which generally depends on ∆𝑥 and ∆𝑡. 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.

1/𝑚 𝐿 -1 1 ∆𝑥 Mass 𝑘𝑔 𝑀 ′  𝑝ℎ  𝑙𝑏𝑚 𝑥 3 Pressure 𝑃𝑎 𝑀 ′ 𝐿 -1 𝑡 -2  𝑝ℎ  𝑙𝑏𝑚 ( ∆𝑥 ∆𝑡 ) 2 Intensity 𝑊/𝑚 2 𝑀 ′ 𝑡 -3  𝑝ℎ  𝑙𝑏𝑚 ( ∆𝑥 ∆𝑡 ) 3
Force per unit volume

𝑁/𝑚 3 𝑀 ′ 𝐿 -2 𝑡 -2  𝑝ℎ  𝑙𝑏𝑚 𝑥 𝑡 2 Power 𝑊 𝑀 ′ 𝐿 2 𝑡 -3  𝑝ℎ  𝑙𝑏𝑚 𝑥 5 𝑡 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.

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 (𝑥 = 0) in order to produce acoustic waves. The spatial dimensions of the cavity are L 𝑥 , L 𝑦 , and L 𝑧 corresponding to length, width, and height, respectively (L 𝑦 = L 𝑧 ). The size of the acoustic source is characterized by its diameter which is defined by 𝑑 = L 𝑧 /3. The aspect ratio of the cavity is L 𝑥 /L 𝑧 = 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 [START_REF] Haydock | Lattice Boltzmann simulations of attenuation-driven acoustic streaming[END_REF] chose a wavelength defined by 50 nodes ( 𝑙𝑏𝑚 = 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 [START_REF] Benhamou | Three-dimensional lattice Boltzmann model for acoustic waves emitted by a source[END_REF], we chose  𝑙𝑏𝑚  = 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  𝑙𝑏𝑚  = 20 . This choice leads to a LBM period 𝑇 𝑙𝑏𝑚 =  𝑙𝑏𝑚 /𝑐 𝑙𝑏𝑚 = 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 𝑐 𝑝ℎ , kinematic viscosity 𝜈 𝑝ℎ , density 𝜌 𝑝ℎ , are given in Table 2. The physical values of ∆𝑥 and ∆𝑡, which can be obtained from Eqs. ( 23) and ( 24), then only depend on the LBM viscosity. If we note that the ratio 𝑐 𝑙𝑏𝑚 /𝑐 𝑝ℎ is small and the kinematic viscosity 𝜈 𝑝ℎ is a small quantity (Tables 1 and2), the only way to get a not too small ∆𝑥 from Eq. ( 23) is to take the smallest possible value of  𝑙𝑏𝑚 . In this study, we chose to fix  𝑙𝑏𝑚 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  𝑙𝑏𝑚 = 0.015, we get ∆𝑥 = 1.698 𝑚 and ∆𝑡 = 2.883 𝑛𝑠, and Eq. ( 27) then gives a physical frequency of 10.011 MHz and a wavelength 𝜆 𝑝ℎ = 33.961 𝑚.

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 𝑀𝑎 = 𝑢 𝑎 /𝑐 𝑙𝑏𝑚 , which must remain weak. It is a consequence of the simplified isoentropic equation of state 𝑝 =  𝑐 2 used by the model. The LBM approach indeed is more stable for simulations in incompressible or weakly compressible mediums [START_REF] Mohamad | Lattice Boltzmann method: Fundamentals and engineering applications with computer codes[END_REF][START_REF] Krüger | The lattice Boltzmann method: Principles and practice[END_REF]. Therefore, for stability reasons, it is generally recommended to use values of 𝑀𝑎 < 0.3 [START_REF] Wang | Lattice Boltzmann simulations of thermal convective flows in two dimensions[END_REF][START_REF] Mohamad | A critical evaluation of force term in lattice Boltzmann method, natural convection problem[END_REF], which gives values of LBM velocity amplitude 𝑢 𝑎 lower than 0.173. The second limitation is the linear wave assumption of the acoustic model described in Eqs. [START_REF] Zhang | Acoustic streaming vortices enable contactless, digital control of droplets[END_REF][START_REF] Guo | Numerical Study on the vapor bubble removal by acoustic streaming in downward-facing pool boiling[END_REF][START_REF] Wu | Acoustic streaming and its applications[END_REF]. 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 [START_REF] Viggen | The lattice Boltzmann method with applications in acoustics[END_REF] recommends a value of  𝑎 less than twice  0 , i.e.  𝑎 < 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 𝑢 𝑎 =  𝑎 𝑐 𝑙𝑏𝑚 [START_REF] Benhamou | Three-dimensional lattice Boltzmann model for acoustic waves emitted by a source[END_REF], we can express the constraint in terms of 𝑢 𝑎 : 𝑢 𝑎 < 2 𝑐 𝑙𝑏𝑚 ≈ 1.15. The velocity amplitudes chosen in this study, 𝑢 𝑎 = 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 𝑉 𝑟𝑒𝑓 = ∆𝑥/∆𝑡 ≈ 589 𝑚/𝑠. We obtain 5.89 and 58.9 𝑚/𝑠 for 𝑢 𝑎 = 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 Δ𝑥 (25) and Δ𝑡 (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.

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 𝑥 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: In the two test cases, the acoustic source vibrates with a period 𝑇 𝑙𝑏𝑚 = 34.64, which gives a frequency of 10.011 MHz in air. The amplitude 𝑢 𝑎 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 𝑡 𝑝ℎ = 1400 𝑡 = 4.037𝑠 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 𝑥-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 (𝑡 𝑝ℎ = 5.767𝑠 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 𝑥 (L x = 28 𝜆), 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 𝐿 𝐹 = 𝑑 2 /(4 𝜆)), the limit between the acoustic near field and far field: 𝐿 𝐹,𝑙𝑏𝑚 = 45 in case 1 and 𝐿 𝐹,𝑙𝑏𝑚 = 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. 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 𝑥-axis, only the acoustic intensity along this direction is calculated. At an instant 𝑡, it is defined as the product of the velocity 𝑢(𝑡) and the acoustic pressure 𝑝(𝑡) (Eq. ( 12)) and calculated with the LBM approach according to Eq. ( 13). The average acoustic intensity 𝐼 𝑎𝑐 = 〈𝐼(𝑡)〉 is then obtained by averaging 𝐼(𝑡) 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. To validate the acoustic intensity presented in Fig. 6, we will compare it to the analytical intensity defined on the central 𝑥-axis of a circular acoustic source. We first recall that, without attenuation, the acoustic pressure 𝑝(𝑥) on this central axis [START_REF] Benhamou | Three-dimensional lattice Boltzmann model for acoustic waves emitted by a source[END_REF][START_REF] Kinsler | Fundamentals of Acoustics[END_REF] can be expressed as:

𝑝(𝑥) = 𝑝 𝑚𝑎𝑥 |𝑠𝑖𝑛 ( 1 2 𝑘𝑥 ( √ 1 + ( 𝑑 2𝑥 ) 2 -1) ) |. (31)
In [START_REF] Rafat | Direct numerical simulations of acoustic streaming in standing wave tubes using the lattice Boltzmann method[END_REF], 𝑝 𝑚𝑎𝑥 is the maximum pressure without attenuation, 𝑥 refers to the distance along the centerline of the transducer, and 𝑘 is the wave vector. 𝑝 𝑚𝑎𝑥 can be expressed in LBM units as a function of the equilibrium density  0 , sound speed 𝑐 𝑙𝑏𝑚 , and velocity amplitude 𝑢 𝑎 . The general expression, given for a baffled acoustic source, is 𝑝 𝑚𝑎𝑥 = 2 0 𝑐 𝑙𝑏𝑚 𝑢 𝑎 . 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 [START_REF] Blackstock | Fundamentals of physical acoustics[END_REF]. In a previous work [START_REF] Benhamou | Three-dimensional lattice Boltzmann model for acoustic waves emitted by a source[END_REF], 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. 𝑝 𝑚𝑎𝑥 is then given by:

𝑝 𝑚𝑎𝑥 =  0 𝑐 𝑙𝑏𝑚 𝑢 𝑎 . ( 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 𝑥-axis given in [START_REF] Kinsler | Fundamentals of Acoustics[END_REF] to account for acoustic losses. We obtain

𝑝̅ = - 𝑗 𝑝 𝑚𝑎𝑥  ( + 𝑗𝑘) 𝑒 𝑗(𝑡) [𝑒 -( +𝑗𝑘) √ 𝑥 2 + 𝑑 2 4 -𝑒 -( +𝑗𝑘)𝑥 ]. ( 33
)
Concerning the acoustic intensity, it is obtained in the general case as

𝐼 𝑎𝑐 = 〈𝑝(𝑡)𝑢(𝑡)〉 = 𝑅𝑒{ 𝑝̅ × 𝑢 ̅ * }/2, ( 34
)
or in the often-used plane wave approximation [START_REF] Kinsler | Fundamentals of Acoustics[END_REF][START_REF] Moudjed | Scaling and dimensional analysis of acoustic streaming jets[END_REF] as

𝐼 𝑎𝑐 = 〈 𝑝(𝑡) 2 〉  0 𝑐 𝑙𝑏𝑚 = 𝑅𝑒{ 𝑝̅ × 𝑝̅ * } 2 0 𝑐 𝑙𝑏𝑚 , (35)
where superscript * refers to the complex conjugate operator and 𝑅𝑒 denotes the real part. To get the acoustic velocity involved in [START_REF] Benhamou | Direct modeling method of acoustic streaming using lattice Boltzmann approach[END_REF], we use the linear Euler equation [START_REF] Kinsler | Fundamentals of Acoustics[END_REF]:

 0 𝜕𝑢(𝑡) 𝜕𝑡 = - 𝜕𝑝(𝑡) 𝜕𝑥 . ( 36
)
Using the complex notation for signals varying in time as 𝑒 𝑗(𝑡) , the velocity 𝑢 ̅ can then be written as

𝑢 ̅ = 𝑗  0  𝜕𝑝̅ 𝜕𝑥 , (37)
and the intensity given by (34) can then be expressed as a function of the complex pressure as

𝐼 𝑎𝑐 = 1 2 0  𝑅𝑒 { 𝑝̅ × (𝑗 𝜕𝑝̅ 𝜕𝑥 ) * }. ( 38
)
The expression [START_REF] Li | Effect of acoustic streaming on heat transfer of porous composite phase change material by using lattice Boltzmann simulation[END_REF] of 𝑝̅ can be introduced in either Eq. [START_REF] Humières | Multiple-relaxation-time lattice Boltzmann models in three dimensions[END_REF] or Eq. ( 35) to get the acoustic intensity along the 𝑥-axis. We obtain the general expression

𝐼 𝑎𝑐 =
and a simpler expression in the plane wave approximation

𝐼 𝑎𝑐 = 𝐼 𝑚 (𝑒 -2𝛼𝑥 + 𝑒 -2 √ 𝑥 2 + 𝑑 2 4 -2 𝑒 -(𝑥+ √ 𝑥 2 + 𝑑 2 4 ) cos (𝑘 ( √ 𝑥 2 + 𝑑 2 4 -𝑥))) (40)
with

𝐼 𝑚 = 𝑝 𝑚𝑎𝑥 2 8 0 𝑐 𝑙𝑏𝑚 . ( 41
)
Note that these expressions use the hypothesis that the attenuation length, related to 𝛼, is large compared with the wavelength 𝜆, more precisely that 𝛼𝜆/𝜋 ≪ 1.

To calculate 𝐼 𝑎𝑐 given by Eqs. ( 39) and ( 40), an estimation of the acoustic attenuation coefficient 𝛼 is needed. Typically, the attenuation of an acoustic wave is due to two physical effects: geometric spreading and dissipation [START_REF] Salomons | Simulation of sound waves using the lattice Boltzmann method for fluid flow: Benchmark cases for outdoor sound propagation[END_REF][START_REF] Salomons | Computational atmospheric acoustics[END_REF]. 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) [START_REF] Lighthill | Acoustic streaming[END_REF][START_REF] Kinsler | Fundamentals of Acoustics[END_REF][START_REF] Salomons | Computational atmospheric acoustics[END_REF]. In the case of waves in gases such as air, the different dissipation effects are taken into account through the following attenuation coefficient [START_REF] Li | Continuum perspective of bulk viscosity in compressible fluids[END_REF][START_REF] Lin | Bulk viscosity model for near-equilibrium acoustic wave attenuation[END_REF]:

𝛼 =  2 2 0 𝑐 𝑙𝑏𝑚 3 ( 4 3
 + ( -1) 2   𝑅 𝑔𝑎𝑠 +  𝐵 ) ( 42
)
where 𝛾 is the ratio of specific heats, 𝐾 is the thermal conductivity, 𝑅 𝑔𝑎𝑠 is the specific gas constant, and  𝐵 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  𝐵 = 3 [START_REF] Moudjed | Caractérisation expérimentale et théorique des écoulements entraînés par ultrasons[END_REF][START_REF] Kinsler | Fundamentals of Acoustics[END_REF]. Concerning thermal effects, they are typically neglected when simulating wave propagation in water [START_REF] Moudjed | Caractérisation expérimentale et théorique des écoulements entraînés par ultrasons[END_REF], 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 [START_REF] Benhamou | Three-dimensional lattice Boltzmann model for acoustic waves emitted by a source[END_REF]. Here, in our calculations, the value  𝑙𝑏𝑚 = 0.015 (related to  𝑝ℎ 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 ( 𝐵 =  𝐵 /) is involved in the relaxation time 𝑠 1 defined as 𝑠 1 = 2/(9 𝜈 𝐵,𝑙𝑏𝑚 + 1). To stabilize the MRT-D3Q19 model, it is recommended to set 𝑠 1 to 1.19 [START_REF] Humières | Multiple-relaxation-time lattice Boltzmann models in three dimensions[END_REF]. Using this value leads to a bulk kinematic viscosity 𝜈 𝐵,𝑙𝑏𝑚 = 7.563 10 -2 .

Using Eq. ( 42) in LBM units with  𝑙𝑏𝑚 = 0.015 and  𝐵,𝑙𝑏𝑚 = 0.07563 and without thermal effects, we obtain 𝛼 𝑙𝑏𝑚 = 8.173 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, 𝛼𝜆/𝜋 ≈ 0.05). These comparisons are shown in Fig. 7 with the plots of the longitudinal profiles of the acoustic intensity along the central 𝑥-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 10 -5 , drops to about 0 and then increases towards the last peak with a maximum value of about 1. 26 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 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 those 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 [START_REF] Li | Lattice Boltzmann method simulation of 3-D natural convection with double MRT model[END_REF] involving the plane wave approximation giving larger amplitudes. The interesting result is that it is the more accurate analytical expression [START_REF] Benhamou | Three-dimensional numerical study of heat transfer enhancement by sound waves using mesoscopic and macroscopic approaches[END_REF] which compares rather well with the numerical intensity profiles, at the peaks and more globally all along the 𝑥-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 [START_REF] Rafat | Direct numerical simulations of acoustic streaming in standing wave tubes using the lattice Boltzmann method[END_REF], the number of peaks in the near field zone is given by 𝑁

= 𝐸𝑛𝑡( 𝑑 2 𝜆 - 1
2 ), with 𝐸𝑛𝑡 the entire part, and there is also the last peak which marks the transition to the far field zone. In LBM units, the wavelength 𝜆 is 20 and the source diameter 𝑑 is 60 in case 1 and 140 in case 2. We then obtain 𝑁 = 1 in case 1 and 𝑁 = 3 in case 2, which corresponds well to what is shown in Fig. 7, the first peak being exactly at 𝑥 = 0 in the two cases. Concerning the last peak, its position without attenuation is usually given by the Fresnel length 𝐿 𝐹 = 𝑑 2 /(4 𝜆) , or 𝐿 𝐹 = 𝑑 2 4 𝜆 -𝜆 4 if the wavelength is not assumed small compared to the source diameter. These two expressions give 𝐿 𝐹 = 45 and 40 in case 1 and 𝐿 𝐹 = 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 𝑥. 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 [START_REF] Li | Lattice Boltzmann method simulation of 3-D natural convection with double MRT model[END_REF] for example, by a Newton method. In this case, we get 𝐿 𝐹 = 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 𝛼 𝑝ℎ can be obtained as 𝛼 𝑝ℎ = 𝛼 𝑙𝑏𝑚 /∆𝑥. With a space step ∆𝑥 = 1.698 𝑚 obtained in our case with 𝜈 𝑝ℎ of air and 𝜈 𝑙𝑏𝑚 = 0.015 (Tables 1 and2), we get 𝛼 𝑝ℎ ≈ 4800 𝑚 -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. [START_REF] Lin | Bulk viscosity model for near-equilibrium acoustic wave attenuation[END_REF] for 𝑓 𝑝ℎ = 10 MHz. Finally, despite the neglect of the thermal effects in our LBM approach, the choice of a bulk viscosity defined from 𝑠 1 = 1.19 allows to get a reasonable estimate of the acoustic attenuation in the air in our microcavities.

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 𝜈 𝑙𝑏𝑚 = 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 ∆𝑥 and time step ∆𝑡 will decrease, in the same way, according to Eqs. ( 25) and [START_REF] Krüger | The lattice Boltzmann method: Principles and practice[END_REF]. For example, the choice of  𝑙𝑏𝑚 = 0.033, the value used by Haydock and Yeomans [START_REF] Haydock | Lattice Boltzmann simulations of attenuation-driven acoustic streaming[END_REF], gives ∆𝑥 = 0.772 𝑚 and ∆𝑡 = 1.31 𝑛𝑠. For a given LBM period, it would also decrease the physical period and increase the physical frequency, giving a frequency 𝑓 𝑝ℎ ≈ 22 MHz. Using  𝑙𝑏𝑚 = 0.1 will still increase these effects leading, for example, to 𝑓 𝑝ℎ ≈ 66.7 MHz and ∆𝑥 = 0.254 𝑚. 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 [START_REF] Pan | Cytosolic delivery of functional proteins in vitro through tunable gigahertz acoustics[END_REF], micromixing [START_REF] Cui | Theoretical and experimental characterizations of gigahertz acoustic streaming in microscale fluids[END_REF], and microscale particle [START_REF] Guo | 3D object detection and tracking based on streaming data[END_REF] and nanoscale particle [START_REF] Guo | Controllable cell deformation using acoustic streaming for membrane permeability modulation[END_REF] manipulation.

Decreasing the LBM viscosity to some extent is also a possibility for the MRT model. For example, using  𝑙𝑏𝑚 = 0.001 gives a relatively larger space step (∆𝑥 = 25.47 𝑚) and a fairly low frequency (𝑓 𝑝ℎ = 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 [START_REF] Moudjed | Scaling and dimensional analysis of acoustic streaming jets[END_REF] remains difficult. For example, for a cavity with a characteristic length of 25 𝑐𝑚 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. [START_REF] Salomons | Simulation of sound waves using the lattice Boltzmann method for fluid flow: Benchmark cases for outdoor sound propagation[END_REF] 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 ∆𝑥 = 0.1 𝑚 is chosen, the time step is ∆𝑡 = 1.7 10 -4 𝑠 , and, for a period 𝑇 𝑙𝑏𝑚 = 40 , the wave frequency is then 𝑓 𝑝ℎ = 147.2 Hz. Such simulation is performed with  𝑙𝑏𝑚 = 0.06 (limitations due to the SRT model used) and corresponds then to  𝑝ℎ = 3.5 𝑚 2 /𝑠, whereas the real value of the viscosity for air is about 1.5 10 -5 𝑚 2 /𝑠. 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 𝜈 𝑙𝑏𝑚 = 0.015 could also be interpreted as an experiment in water. The choice of a bulk kinematic viscosity 𝜈 𝐵,𝑙𝑏𝑚 = 7.563 10 -2 is a little too high, as previous studies rather take this bulk viscosity as  𝐵 = 3 [START_REF] Moudjed | Caractérisation expérimentale et théorique des écoulements entraînés par ultrasons[END_REF][START_REF] Kinsler | Fundamentals of Acoustics[END_REF] and thermal effects are negligible [START_REF] Moudjed | Caractérisation expérimentale et théorique des écoulements entraînés par ultrasons[END_REF], so that wave attenuation will be a little increased. Nevertheless, using the values 𝑐 𝑝ℎ = 1480 𝑚/𝑠, 𝜈 𝑝ℎ = 10 -6 𝑚 2 /𝑠 for water, we get ∆𝑥 ≈ 26 𝑛𝑚, ∆𝑡 ≈ 10 -11 𝑠, and, for a wavelength 𝜆 𝑙𝑏𝑚 = 20, the wave frequency is then 𝑓 𝑝ℎ = 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.

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 𝑇 𝑙𝑏𝑚 = 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. 𝑈 in the 𝑥 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 𝐼 𝑎𝑐 , which, according to [START_REF] Lei | Outer acoustic streaming flow driven by asymmetric acoustic resonances[END_REF] and [START_REF] Haydock | Lattice Boltzmann simulations of acoustic streaming[END_REF], can be written as:

𝐹 𝑎𝑐 = 2𝛼 𝑐 𝑙𝑏𝑚 〈∑ 𝑓 𝑖 18 𝑖=0
𝑐 𝑥,𝑖 〉, [START_REF] Bouzidi | Momentum transfer of a Boltzmann-lattice fluid with boundaries[END_REF] the average being taken over 400 iterations, as indicated for 𝐼 𝑎𝑐 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 [START_REF] Moudjed | Near-field acoustic streaming jet[END_REF].

Two cases will be considered. The first case corresponds to a rather small wave amplitude 𝑢 𝑎 = 0.01 and the second case to a larger amplitude 𝑢 𝑎 = 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.

Vibration amplitude 𝒖 𝒂 = 𝟎. 𝟎𝟏 a) Direct method

In order to estimate the time evolution of the waves necessary to get a stabilized oscillatory state for 𝑢 𝑎 = 0.01, we plot the time variation of the instantaneous acoustic velocity 𝑢 in two different points along the central 𝑥-axis in Fig. 8. For the point at mid-length (𝑥 = L 𝑥 /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 𝑥 = 9L 𝑥 /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. 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.

b) Indirect method for validation

The acoustic force 𝐹 𝑎𝑐 , 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 𝑥-axis, close to the source, and the intensity of the force strongly decreases for larger 𝑥 and outside the acoustic beam. The time variation of the streaming velocity 𝑈 in two different points along the central 𝑥-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 𝑈 along the central 𝑥-axis as well as its transverse profile along the central 𝑦-axis, the comparison is very good, with differences not exceeding a few percent. Fig. 13(A) shows that the maximum velocity along the central 𝑥-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.

Vibration amplitude 𝒖 𝒂 = 𝟎. 𝟏 a) Direct method

For this larger velocity vibration amplitude 𝑢 𝑎 = 0.1, the time variation of the instantaneous acoustic velocity 𝑢 in two different points along the central 𝑥-axis is plotted in Fig. 14. We see that for both points (at 𝑥 = L 𝑥 /2 and 𝑥 = 9L 𝑥 /10), the time evolution of the waves necessary to get a stabilized oscillatory state is strongly reduced compared with the case at 𝑢 𝑎 = 0.01. Only 20 000 iterations are now needed, and this time is even less for regions with higher velocities as at 𝑥 = L 𝑥 /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 𝑦 = L 𝑦 /2 and the transverse plane at 𝑥 = L 𝑥 /2 in Figs. 15 and16. For 𝑢 𝑎 = 0.1, the flow in the cavity appears quite different from that obtained previously for a smaller oscillation amplitude 𝑢 𝑎 = 0.01 . The longitudinal component 𝑈 of the streaming has now strong values up to the end wall, indicating the presence of a jet all along the 𝑥-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.

b) Indirect method for validation

The acoustic force 𝐹 𝑎𝑐 , which is injected in the LBM code in this case, is the same as that calculated for 𝑢 𝑎 = 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 𝑈 in two different points along the central 𝑥-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 𝑢 𝑎 = 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 𝑢 𝑎 = 0.01.

Discussion

The simulations presented in the previous subsections gave the streaming flow structure obtained for two values of the velocity oscillation amplitude 𝑢 𝑎 . They also indicate the intensity of the streaming flow, which, according to Figs. 9 and 15, has maximum values around 𝑢 𝑙𝑏𝑚,𝑚𝑎𝑥 = 1.5 10 -3 and 4 10 -2 for 𝑢 𝑎 = 0.01 and 0.1, respectively. Using the reference velocity given by 𝑐 𝑝ℎ /𝑐 𝑙𝑏𝑚 , indications on the physical streaming velocities thus created can be obtained.

If we consider the experiment in air, 𝑐 𝑝ℎ /𝑐 𝑙𝑏𝑚 ≈ 589 𝑚/𝑠, and the maximum physical velocities will be 𝑢 𝑝ℎ,𝑚𝑎𝑥 = 0.88 and 3.84 𝑚/𝑠 for 𝑢 𝑎 = 0.01 and 0.1, respectively. In the case of an experiment in water, 𝑐 𝑝ℎ /𝑐 𝑙𝑏𝑚 ≈ 2563 𝑚/𝑠, and the maximum physical velocities will then be 𝑢 𝑝ℎ,𝑚𝑎𝑥 = 23.56 and 102.5 𝑚/𝑠 for 𝑢 𝑎 = 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.

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 [START_REF] Moudjed | Near-field acoustic streaming jet[END_REF]. 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 𝑢 𝑎 , the toroidal roll is rather developed close to the transducer (small 𝑢 𝑎 ) or in the whole cavity, with the main central jet spreading on the end wall (larger 𝑢 𝑎 ). 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  𝑙𝑏𝑚 cannot be chosen very small leads to small values of the space and time steps ∆𝑥 and ∆𝑡 (Eqs. [START_REF] Mohamad | Lattice Boltzmann method: Fundamentals and engineering applications with computer codes[END_REF][START_REF] Krüger | The lattice Boltzmann method: Principles and practice[END_REF]). 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 [START_REF] Renard | A linear stability analysis of compressible hybrid lattice Boltzmann methods[END_REF]. 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.

review & editing. Valéry Botton: Supervision, Investigation, Visualization, Formal analysis, Writingreview & editing.

𝑠 0 = 6 ( 10 )

 0610 𝑠 3 = 𝑠 5 = 𝑠 7 = 1.0, 𝑠 1 = 1.19, 𝑠 2 = 𝑠 10 = 𝑠 12 = 1.4 s 9 = s 11 = s 13 = s 14 = s 15 = 1/(3ν + 0.5) 𝑠 4 = 𝑠 6 = 𝑠 8 = 1.2, 𝑠 16 = 𝑠 17 = 𝑠 18 = 1.The solution of the lattice Boltzmann equation allows to determine the distribution function 𝑓 𝑖 . The density and velocities can then be calculated as:

Fig. 2 .

 2 Fig. 2. Geometry of the physical problem studied.

Case 1 :LCase 2 :L

 12 L 𝑥,𝑝ℎ = 407.541 𝑚 = 12 𝜆 𝑝ℎ , L 𝑦,𝑝ℎ = L 𝑧,𝑝ℎ = 305.656 𝑚 = 9 𝜆 𝑝ℎ , 𝑑 𝑝ℎ = L 𝑥,𝑝ℎ = 950.92 𝑚 = 28 𝜆 𝑝ℎ , L 𝑦,𝑝ℎ = L 𝑧,𝑝ℎ = 713.19 𝑚 = 21 𝜆 𝑝ℎ , 𝑑 𝑝ℎ =

Fig. 3 .Fig. 4 .

 34 Fig. 3. Longitudinal acoustic velocity field 𝑢 produced after 1400 iterations (4.037𝑠 in the air) by a circular acoustic source vibrating with a period 𝑇 𝑙𝑏𝑚 = 34.64 (10.011 MHz frequency in the air) with a LBM amplitude 𝑢 𝑎 = 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

Fig. 5 .

 5 Fig. 5. Profiles along the central 𝑥-axis for the longitudinal acoustic velocity 𝑢 produced by a circular source vibrating with a period 𝑇 𝑙𝑏𝑚 = 34.64 (10.011 MHz frequency in the air) with an amplitude 𝑢 𝑎 = 0.01 in two microcavities of different sizes: (A) case 1, as in Fig. 3; (B) case 2, as in Fig. 4.

Fig. 6 .

 6 Fig. 6. 3D half view of the acoustic intensity field 𝐼 𝑎𝑐 produced by a circular source vibrating with a period 𝑇 𝑙𝑏𝑚 = 34.64 (10.011 MHz frequency in the air) with an amplitude 𝑢 𝑎 = 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 × 10 -7 in (A) and 2 × 10 -8 in (B). In both figures, the isovalues given in the scale are also shown in the central plane at 𝑦 = L 𝑦 /2.

Fig. 7 .

 7 Fig. 7. Profiles of acoustic intensity 𝐼 𝑎𝑐 along the central 𝑥-axis of the cavity for a circular source vibrating with a period 𝑇 𝑙𝑏𝑚 = 34.64 (10.011 MHz frequency in the air) with an amplitude 𝑢 𝑎 = 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 𝐿 𝐹 = 𝑑 2 /(4 𝜆) (orange vertical line) is compared to its position taking into account attenuation and obtained from (40) (black vertical line).

Fig. 8 .

 8 Fig. 8. Time variation of the acoustic velocity 𝑢 at different positions along the central 𝑥-axis for a circular acoustic source vibrating with an amplitude 𝑢 𝑎 = 0.01: 𝑥 = L 𝑥 /2 (black) and 𝑥 = 9L 𝑥 /10 (blue).

Fig. 9 .

 9 Fig. 9. 2D views of the streaming velocity 𝑈 generated by a circular acoustic source vibrating with an amplitude 𝑢 𝑎 = 0.01: (A) longitudinal section at 𝑦 = L 𝑦 /2 and (B) transverse section at 𝑥 = L 𝑥 /2.

Fig. 10 .

 10 Fig. 10. Velocity vectors (A) and streamlines (B) in the central plane at 𝑦 = L 𝑦 /2 for the streaming velocity generated by a circular acoustic source vibrating with an amplitude 𝑢 𝑎 = 0.01. The streamlines are used to depict the structure of the flow, but do not give information on its intensity.

Fig. 11 .

 11 Fig. 11. 3D half view of the acoustic force field 𝐹 𝑎𝑐 produced by a circular source vibrating with an amplitude 𝑢 𝑎 = 0.01. The force is averaged over 400 iterations. The 3D blue surface corresponds to the smallest isovalue in the scale, i.e. 1 × 10 -8 . The isovalues given in the scale are also shown in the central plane at 𝑦 = L 𝑦 /2.

Fig. 12 .

 12 Fig. 12. Time variation of the streaming velocity 𝑈 obtained by the LBM code with an imposed acoustic streaming force calculated for a source velocity amplitude 𝑢 𝑎 = 0.01. The force is depicted in Fig. 11. The streaming velocity is given at two positions along the central 𝑥-axis: 𝑥 = L 𝑥 /2 (black) and 𝑥 = 9L 𝑥 /10 (blue).

Fig. 13 .

 13 Fig. 13. Longitudinal profile along the central 𝑥-axis (A) and transverse profile along the central 𝑦axis (B) for the streaming velocity 𝑈 generated by a circular acoustic source vibrating with an amplitude 𝑢 𝑎 = 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.

Fig. 14 .

 14 Fig. 14. Time variation of the acoustic velocity 𝑢 at different positions along the central 𝑥-axis for a circular acoustic source vibrating with an amplitude 𝑢 𝑎 = 0.1: 𝑥 = L 𝑥 /2 (black) and 𝑥 = 9L 𝑥 /10 (blue).

Fig. 15 .

 15 Fig. 15. 2D views of the streaming velocity 𝑈 generated by a circular acoustic source vibrating with an amplitude 𝑢 𝑎 = 0.1: (A) longitudinal section at 𝑦 = L 𝑦 /2 and (B) transverse section at 𝑥 = L 𝑥 /2.

Fig. 16 .

 16 Fig. 16. Velocity vectors (A) and streamlines (B) in the central plane at 𝑦 = L 𝑦 /2 for the streaming velocity generated by a circular acoustic source vibrating with an amplitude 𝑢 𝑎 = 0.1. The streamlines are used to depict the structure of the flow, but do not give information on its intensity.

Fig. 17 .

 17 Fig. 17. 3D half view of the acoustic force field 𝐹 𝑎𝑐 produced by a circular source vibrating with an amplitude 𝑢 𝑎 = 0.1. The force is averaged over 400 iterations. The 3D blue surface corresponds to the smallest isovalue in the scale, i.e. 3 × 10 -7 . The isovalues given in the scale are also shown in the central plane at 𝑦 = L 𝑦 /2.

Fig. 18 .

 18 Fig. 18. Time variation of the streaming velocity 𝑈 obtained by the LBM code with an imposed acoustic streaming force calculated for a source velocity amplitude 𝑢 𝑎 = 0.1. The force is depicted in Fig. 17. The streaming velocity is given at two positions along the central 𝑥-axis: 𝑥 = L 𝑥 /2 (black) and 𝑥 = 9L 𝑥 /10 (blue).

  Fig. 18. Time variation of the streaming velocity 𝑈 obtained by the LBM code with an imposed acoustic streaming force calculated for a source velocity amplitude 𝑢 𝑎 = 0.1. The force is depicted in Fig. 17. The streaming velocity is given at two positions along the central 𝑥-axis: 𝑥 = L 𝑥 /2 (black) and 𝑥 = 9L 𝑥 /10 (blue).

Fig. 19 .

 19 Fig. 19. Longitudinal profile along the central 𝑥-axis (A) and transverse profile along the central 𝑦axis (B) for the streaming velocity 𝑈 generated by a circular acoustic source vibrating with an amplitude 𝑢 𝑎 = 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.

  2 , 𝑗 𝑥 , 𝑞 𝑥 , 𝑗 𝑦 , 𝑞 𝑦 , 𝑗 𝑧 , 𝑞 𝑧 , 3𝑝 𝑥𝑥 , 3𝜋 𝑥𝑥 , 𝑝 𝑤𝑤 , 𝜋 𝑤𝑤 , 𝑝 𝑥𝑦 , 𝑝 𝑦𝑧 , 𝑝 𝑧𝑥 , 𝑚 𝑥 , 𝑚 𝑦 , 𝑚 𝑧 )

𝑇

.

  The matrix 𝑆 is diagonal. It defines each relaxation time 𝑠 𝑖 of the corresponding macroscopic quantities. It can be written as a function of the nineteen 𝑠 𝑖 as:𝑆 = diag(𝑠 0 , 𝑠 1 , 𝑠2 , 𝑠 3 , 𝑠 4 , 𝑠 5 , 𝑠 6 , 𝑠 7 , 𝑠 8 , 𝑠 9 , 𝑠 10 , 𝑠 11 , 𝑠 12 , 𝑠 13 , 𝑠 14 , 𝑠 15 , 𝑠 16 , 𝑠 17 , 𝑠 18 )

	𝑚 9 𝑒𝑞 =	2𝑗 𝑥 2 -(𝑗 𝑦 2 + 𝑗 𝑧 2 )  0
	𝑚 10 𝑒𝑞 = -	1 2	2𝑗 𝑥 2 -(𝑗 𝑦 2 + 𝑗 𝑧 2 )  0	(8)
	𝑚 11 𝑒𝑞 =	𝑗 𝑦 2 -𝑗 𝑧 2  0	𝑚 12 𝑒𝑞 = -	1 2	𝑗 𝑦 2 -𝑗 𝑧 2  0
		𝑚 13 𝑒𝑞 =	𝑗 𝑥 𝑗 𝑦  0
		𝑚 14 𝑒𝑞 =	𝑗 𝑦 𝑗 𝑧  0
		𝑚 15 𝑒𝑞 = 𝑒𝑞 = 0 𝑗 𝑥 𝑗 𝑧  0 𝑚 16 𝑒𝑞 = 0 𝑚 17 𝑚 18 𝑒𝑞 = 0
			4 𝑒𝑞 = -𝑚 5 𝑒𝑞 = 𝑗 𝑦 2 3	𝑗 𝑥
		𝑚 6 𝑒𝑞 = -𝑚 7 𝑒𝑞 = 𝑗 𝑧 2 3	𝑗 𝑦
		𝑚 8 𝑒𝑞 = -	2 3	𝑗 𝑧

Table 1 .

 1 LBM parameters used in the simulations of acoustic waves propagation.

	𝒙 𝒍𝒃𝒎 𝒕 𝒍𝒃𝒎	𝑻 𝒍𝒃𝒎	𝝀 𝒍𝒃𝒎	𝒄 𝒍𝒃𝒎	 𝒍𝒃𝒎	 𝒍𝒃𝒎	𝒖 𝒂
	1	1	34.64	20	1/√3	0.015	1	0.01 -0.1

Table 2 .

 2 Physical characteristics associated with the LBM simulations of acoustic waves propagation in air and properties of air considered (at 20 0 C).

		Physical characteristics			Properties of air
	𝟏. 𝟔𝟗𝟖	2.883	10.011	33.961	340	1.5 10 -5	1.204

𝒙 (𝒎) 𝑡(𝑛𝑠) 𝑓 𝑝ℎ (𝑀𝐻𝑧) 𝜆 𝑝ℎ (𝑚) 𝑐 𝑝ℎ (m/s)  𝑝ℎ (𝑚 2 /s)  𝑝ℎ (kg/𝑚 3 )

Table 3 .

 3 Factors 𝑎 used to convert acoustic quantities from the LBM scale to the physical scale ( 𝑝ℎ = 𝑎  𝑙𝑏𝑚 ).

	Physical	SI unit	Dimensions	Conversion factor from
	quantity			LBM to physical scale
	Distance	𝑚	𝐿	∆𝑥
	Time	𝑠	𝑡	∆𝑡
	Velocity	𝑚/𝑠	𝐿𝑡 -1	∆𝑥
				∆𝑡
	Attenuation
	coefficient

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.

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.

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