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

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.


F
Body force (N/m 3 ) g Gravitational vector (m/s  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][3][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][13][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][26][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 1990 s.We can mention some of the various works that are more focused on acoustic streaming.Stansell and Greated [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.

Numerical method
Our simulations of streaming 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.

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] In the D3Q19 lattice, particles can take the following 19 LBM velocities 36,38]: 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] ( 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 ρ and the momentums j x , j y , and j z along the x, y, and z axes.m can be represented mathematically as m = ρ, e, e 2 , j x , q x , j y , q y , j z , q z , 3 p x x , 3π x x , p ww , π ww , p xy , p yz , p zx , m x ,m y , m z T . (4) The matrix M −1 can be obtained by calculating the inverse of the transformation matrix M, which is defined as [36,38,39] 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] where V is the velocity vector ( V = ( u, v, w)) and W i represents the discretization weights, which are given as: The nineteen elements of m eq i found after performing the calculations can be given in terms of density and momentums as follows [38,39]: 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]: 19, s 2 = s 10 = s 12 = 1.4 s 9 = s 11 =s 13 =s 14 =s 15 = 1/(3ν + 0.5) s 4 = s 6 = s 8 = 1.2, s 16 = s 17 = s 18 = 1.6 (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: The acoustic intensity is also an interesting quantity in our study.Its instantaneous value in the x−direction can be expressed as [42] Using the formula p = ρc 2 lbm , the intensity can then be computed in the LBM units as: 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. [43].For a transducer considered to vibrate at a certain angular frequency ω and with a given velocity amplitude u a , the idea is to write (where f ī is the distribution function in the opposite direction ( ī = −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]: where ρ 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.

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 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] (19) and Eq. ( 1) becomes [37] 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 [25][26][27] 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 ī (x B , t) is the known distribution function in the opposite direction ( ī = −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][51][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]: 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 α 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.

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][59][60] is that which allows direct relationships between the LBM and physical acoustic quantities by means of the space step x (m) and the time step 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 x and 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: 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, ν 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 ν 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 ν lbm can be used without affecting the stability of the calculations.Once the value of ν lbm is chosen, Eqs. ( 23) and ( 24) can be used to obtain x and t: Using the values of x and t, the other quantities such as the physical period (T ph ) and wavelength (λ ph ) can be easily calculated from those used in the LBM simulations: 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 = ρc 2 first enable to express the acoustic intensity as I ac = ρuc 2 .In physical and LBM units, this intensity is given by I ac, ph = ρ ph u ph c 2 ph and I ac,lbm = ρ lbm u lbm c 2 lbm .Dividing the first equation by the second and using u ph /u lbm = c ph /c lbm = x/ t allows us to express I ac, ph as a function of I ac,lbm : The acoustic force can then be defined from the intensity I ac as [5,61]: Using (28) and α ph /α lbm = 1/ x we can write: The conversion described above makes it possible to relate directly the physical quantities (φ ph ) to the LBM quantities (φ lbm ) in the form φ ph = a φ lbm , where the parameter a is the conversion factor which generally depends on x and 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 Sect. 4.

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 (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 (λ 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 λ 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 λ lbm = 20.This choice leads to a LBM period T lbm = λ lbm /c lbm = 34.64.
To find the physical characteristics associated with our simulations in air, we used the direct conversion described in Sect. 2. The properties of air, sound speed c ph , kinematic viscosity ν ph , density ρ ph , are given in Table 2.The physical values of x and 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 ν ph is a small quantity (Tables 1 and 2), the only way to get a not too small x from Eq. ( 23) is to take the smallest possible value of ν lbm .In this study, we chose to fix ν 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 ν lbm = 0.015, we get x = 1.698 μm and t = 2.883 ns, and Eq. ( 27) then gives a physical frequency of 10.011 MHz and a wavelength λ ph = 33.961μ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 = ρ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)(16)(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 ρ a less than twice ρ 0 , i.e. ρ 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 = ρ a c lbm [36], we can express the constraint in terms of u a : u a < 2 c lbm ≈ 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 re f = x/ t ≈ 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 x (25) and t (26) (as proposed for air in Table 2) and then to obtain the conversion factors (described in Sect.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 Sects.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 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:  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.
Figure 3 shows the instantaneous acoustic velocity field obtained in case 1 after 1400 iterations, which corresponds to a physical time t ph = 1400 t = 4.037μs for the propagation in air. Figure 3A gives a 3D view, whereas Fig. 3B 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. 3A 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μ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 λ), 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.
Figure 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 λ)), 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. 5A 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 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  3; B case 2, as in Fig. 4 p(t) (Eq.( 12)) and calculated with the LBM approach according to Eq. ( 13).The average acoustic intensity I ac = I (t) 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.
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: 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 ρ 0 , sound speed c lbm , and velocity amplitude u a .The general expression, given for a baffled acoustic source, is p max = 2ρ 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: 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 p= − jπ p max λ (α+ jk) e j(ωt) e −(α+ jk) Concerning the acoustic intensity, it is obtained in the general case as or in the often-used plane wave approximation [42,61] as 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]: Using the complex notation for signals varying in time as e j(ωt) , the velocity ū can then be written as and the intensity given by (34) can then be expressed as a function of the complex pressure as The expression (33) of p can be introduced in either Eqs.(38) or (35) to get the acoustic intensity along the x-axis.We obtain the general expression and a simpler expression in the plane wave approximation with 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 I ac 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 [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]: where γ is the ratio of specific heats, K is the thermal conductivity, R gas is the specific gas constant, and μ 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 μ B = 3μ [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 ν lbm = 0.015 (related to ν 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 (ν B = μ B /ρ) is involved in the relaxation time s 1 defined as s 1 = 2/(9 ν 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 ν B,lbm = 7.563 10 −2 .
Using Eq. ( 42) in LBM units with ν lbm = 0.015 and ν B,lbm = 0.07563 and without thermal effects, we obtain α lbm = 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 x-axis, which is the axis of the source.In case 1 (Fig. 7A), 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. 7B), 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 (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 = Ent ( d 2 λ − 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 λ 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 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 α ph can be obtained as α ph =α lbm / x.With a space step x= 1.698 μm obtained in our case with ν ph of air and ν lbm = 0.015 (Tables 1 and 2), we get α ph ≈ 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 s 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 ν 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 x and time step t will decrease, in the same way, according to Eqs. ( 25) and (26).For example, the choice of ν lbm = 0.033, the value used by Haydock and Yeomans [30], gives x= 0.772 μm and 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 ≈ 22 MHz.Using ν lbm = 0.1 will still increase these effects leading, for example, to f ph ≈ 66.7 MHz and x = 0.254 μ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 ν lbm = 0.001 gives a relatively larger space step ( x= 25.47 μ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 x= 0.1 m is chosen, the time step is t= 1.7 10 −4 s, and, for a period T lbm = 40, the wave frequency is then f ph = 147.2Hz.Such simulation is performed with ν lbm = 0.06 (limitations due to the SRT model used) and corresponds then to ν 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 mediumfrequency waves in realistic fluids.
The present calculations with ν lbm = 0.015 could also be interpreted as an experiment in water.The choice of a bulk kinematic viscosity ν B,lbm = 7.563 10 −2 is a little too high, as previous studies rather take this bulk viscosity as μ B = 3μ [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, ν ph = 10 −6 m 2 /s for water, we get x≈ 26 nm, t≈10 −11 s, and, for a wavelength λ 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.

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: 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.
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 y /2 and the transverse plane at x = L x /2 in Fig. 9, whereas velocity vector plots and streamlines in the plane at y = L y /2 are shown in Fig. 10.We see that, for u a = 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.

(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.Figure 13A shows that the maximum velocity along the central x-axis is reached quite quickly, at only a quarter of the cavity.Figure 13B points out the weakness of the return flow compared to the central flow.

Vibration amplitude u
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 Fig. 13 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 a = 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 and x = 9 L 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 60,000 to 100,000.
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.

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.

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 ≈ 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 ≈ 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.

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 Depending on the velocity vibration amplitude u a , the toroidal roll is rather developed close to the transducer (small u ) 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 ν lbm cannot be chosen very small leads to small values of the space and time steps x and 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.

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

Fig. 3 Fig. 4 2 Fig. 5
Fig. 3 Longitudinal acoustic velocity field u produced after 1400 iterations (4.037 μs in the air) by a circular acoustic source vibrating with a period T lbm = 34.64(10.011MHz 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 −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 (colour figure online)

Fig. 6
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.011MHz 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 × 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 y = L y /2

Fig. 7
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.011MHz 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 λ) (orange vertical line) is compared to its position taking into account attenuation and obtained from (40) (black vertical line) (colour figure online)

Fig. 8
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 a = 0.01: x = L x /2 (black) and x = 9L x /10 (blue) (colour figure online)

Fig. 9 2 Fig. 10
Fig. 9 2D views of the streaming velocity U generated by a circular acoustic source vibrating with an amplitude u a = 0.01: A longitudinal section at y = L y /2 and B transverse section at x = L x /2

Fig. 11 Fig. 12
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 × 10 −8 .The isovalues given in the scale are also shown in the central plane at y = L y /2 (colour figure online)

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

Fig. 16 Fig. 17
Fig. 16 Velocity vectors (A) 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 Fig. 173D 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 × 10 −7 .The isovalues given in the scale are also shown in the central plane at y = L y /2 (colour figure online)

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

Fig. 19
Fig. 18 Time variation of the streaming velocity U obtained by the code with an imposed acoustic streaming force calculated for a source velocity amplitude u a = 0.1.The force is depicted in Fig. 17.The streaming velocity is given at two positions along the central x-axis: x = L x /2 (black) and x = 9L x /10 (blue) (colour figure online) 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= 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 ) (

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

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

Table 3
Factors a used to convert acoustic quantities from the LBM scale to the physical scale (φ ph = a φ lbm )