Acoustic lattice Boltzmann model for immiscible binary fluids with a species-dependent impedance

In this paper we consider the propagation of an acoustic wave or pulse in an immiscible binary fluid in which the speed of sound is different for each component. Aside from purely theoretical interests this phenomenon has a large number of applications, ranging from medical to industrial. Examples include foetal imagining, diagnosing and repairing muscle damage, removing gas from a fluid by cavitation, reducing damage to propeller heads caused by high-speed bubbles, and focusing of acoustic beams in, for example, lithotriptors. The simulations are performed using the lattice Boltzmann model LBM 1 . The LBM has been developed as a technique for simulating fluid flow and has proved to be particularly advantageous for simulating binary-fluid mixtures 2–6 . It has also been shown that the LBM can be applied to simulate acoustic waves 7–9 and the interaction between acoustic and velocity fields 10–14 . The LBM simulates the acoustic and flow phenomenon within the same model and so does not require the coupling of two distinct models. The LBM for acoustic waves applied in this paper is capable of simulation nonlinear waves with Mach number up to around 0.02. Higher Mach numbers can be achieved using a finite-difference lattice Boltzmann scheme see, for example, 15,16 . The approach presented here could be implemented in such a model.


I. INTRODUCTION
In this paper we consider the propagation of an acoustic wave or pulse in an immiscible binary fluid in which the speed of sound is different for each component.Aside from purely theoretical interests this phenomenon has a large number of applications, ranging from medical to industrial.Examples include foetal imagining, diagnosing and repairing muscle damage, removing gas from a fluid by cavitation, reducing damage to propeller heads caused by high-speed bubbles, and focusing of acoustic beams in, for example, lithotriptors.The simulations are performed using the lattice Boltzmann model ͑LBM͒ ͓1͔.The LBM has been developed as a technique for simulating fluid flow and has proved to be particularly advantageous for simulating binary-fluid mixtures ͓2-6͔.It has also been shown that the LBM can be applied to simulate acoustic waves ͓7-9͔ and the interaction between acoustic and velocity fields ͓10-14͔.The LBM simulates the acoustic and flow phenomenon within the same model and so does not require the coupling of two distinct models.The LBM for acoustic waves applied in this paper is capable of simulation nonlinear waves with Mach number up to around 0.02.Higher Mach numbers can be achieved using a finite-difference lattice Boltzmann scheme ͑see, for example, ͓15,16͔.The approach presented here could be implemented in such a model.

II. NUMERICAL MODEL
The binary fluid model applied here was originally developed by ͓2,3͔.The fluid system evolves on a regular twodimensional, hexagonal grid with six links i = 1 , . . ., 6, joining each site to its six nearest neighbors.A rest link i =0 is also defined.The fluid is described in terms of two distribution functions: f i representing the combined fluid and ⌬ i representing the order parameter of the binary system.The combined fluid density , the fluid velocity u, and the order parameter ⌬ are found from these distribution functions as The evolution of the distribution functions is determined by two lattice Boltzmann equations and where and ⌬ are independent relaxation parameters.Following ͓2,3͔, the form for the equilibrium distribution functions f i and ⌬ i is selected to describe a binary fluid consisting of two ideal gases with a repulsive interaction energy: and where Here ⌫ is the mobility, ⌬ is the chemical-potential difference between the two components, T is the temperature, is the interfacial energy, and ⌳ measures the strength of the interaction.This binary-fluid model can be shown to satisfy the Navier-Stokes equation ͓3͔ and the convection-diffusion equation

͑9͒
and

͑10͒
For T Ͻ⌳/ 2 the mixture separates into two immiscible phases ͓3͔.
In the LBM the density is not constrained to be constant and so acoustic waves can be simulated by introducing a density variation which will propagate at speed c s , where the pressure term in the Navier-Stokes equation can be expressed as p = c s 2 .This has been investigated for a single-phase fluid ͓7,8͔ for Mach numbers up to 0.02 including nonlinear acoustic phenomena.The limit on the Mach number arises from the assumption that variations in the density are negligibly small.This assumption is used in deriving the Navier-Stokes equation from the lattice Boltzmann equation.
An acoustic wave can be introduced in the same manner to the binary-fluid model described above.This gives a speed of sound c s = ͱ p 0 / .Additionally, for the binary-fluid model, introducing a density variation will perturb the terms in the equilibrium distribution functions.This will occur through the nonconstant density ͑as it does for a single-phase fluid͒ and also through the density gradient terms in Eq. ͑6͒.Consider a density ͑pressure͒ variation of amplitude ␦ and spatial extent ␦r.To maintain the simulation in an approximately incompressible limit, previous work has applied the limit ␦ / Ͻ 0.02.Also, in determining the Navier-Stokes equation from a Chapman-Enskog expansion it is assumed that the spatial resolution of the simulation is significantly larger than the lattice spacing-say, ␦r Ͼ 10.Thus, within the normal limits of the LBM, ٌ 2 ␦ ϳ ␦ / ͑␦r͒ 2 Ͻ 2 ϫ 10 −4 and ٌ͑͒ 2 ϳ͑␦͒ / ͑␦r͒ 2 Ͻ 4 ϫ 10 −6 for Ӎ 1.These terms are negligible compared to the dominant terms in Eq. ͑6͒, in particular the gradients of ⌬ at the interface.Also, from Eq. ͑10͒, p 0 Ӎ T away from the interface, giving c s Ӎ ͱ T.
Further, the speed of sound can be varied by introducing a forcing term dependent on the density gradient ͓17͔.This is introduced by modifying the lattice Boltzmann equation to which gives a modified speed of sound, c e = ͱ ͑c s 2 − ␤͒.Good agreement has been found with theory for the D2Q9 model ͑which has c s 2 =1/3͒ for −2 / 3 Ͻ ␤ Ͻ 1 / 3. We note that the difference between the form of the force term used in Eq. ͑11͒ and ͓17͔ is due to the form of the underlying lattice ͑see ͓18͔͒.
Here this model is extended such that the speed of sound is different in the two components of the binary fluid.This is done by replacing Eq. ͑2͒ with Eq. ͑11͒, where ␤ is now a function of the order parameter, ␤ = ␤(⌬͑x , t͒), such that the speed of sound is different in the two phases.In principle any function of ⌬͑x , t͒ may be applied.Here ␤͑x,t͒ = ͭ 0, ⌬ Յ 0, ␣, ⌬ Ͼ 0, ͑12͒ is used.Assuming constant density ͑see Sec.III͒, this gives an impedance difference of and a ratio of impedances where ⑀ ⌬ + and ⑀ ⌬ − are the acoustic impedances of the fluids represented by positive and negative ⌬, respectively.The acoustic impedance for each fluid is defined explicitly in the Appendix.

III. RESULTS
In this section the propagation of acoustic signals through an immiscible binary fluid mixture, with an impedance mismatch at the fluid interface, will be considered.Acoustic pulses and waves were generated at the left boundary ͑x =0͒ following ͓17͔.The distribution functions f i were set equal to their equilibrium values, f ¯i( = 0 + ␦͑t͒ , u =0), where 0 is the ambient density and ␦͑t͒ describes the pulse or acoustic wave being implemented.Satisfactory results were obtained by neglecting the ٌ 2 and ٌ 2 ⌬ terms.Thus f 0 ͑x =0,t͒ = ͑1−2T͓͒ 0 + ␦͑t͔͒ and f i ͑x =0,t͒ = T͓ 0 + ␦͑t͔͒ / 3 for i =1,2, ... ,6.The distribution functions g i were set to f i ⌬ ϱ , where ⌬ ϱ is the ͑measured͒ value of ⌬ far from the interface.At the right boundary the unknown distribution functions were determined ͓19͔ using a backward linear interpolation.Periodic boundary conditions were applied at the top and bottom boundaries.
In the following simulations the temperature, the interaction strength, and the mobility were fixed at T = 0.4, ⌳ = 1.1, and ⌫ = 0.1, respectively, to simulate two immiscible fluids.Two values for the interfacial energy were used: = 0.001 and 0.0001.The form of a plane interface is shown in Fig. 1.In both cases ⌬ varies over one lattice unit and there is no distinguishable difference between the two values of .There is also a dip in the density over the interface of approximately 0.3% for = 0.001 and 0.03% for = 0.0001.A density difference is also present across the interface of a circular bubble.The magnitude of this difference can be determined from Laplace's equation where r is the radius of the bubble, ⌬p is the pressure difference across the interface, and is the surface tension which is a function of the binary fluid parameters ͓4,20͔.For the parameters used here this gives a negligible density difference, significantly smaller than the density dip at the interface.
The effect of the small density dip at the interface was investigated by simulating an acoustic pulse propagating perpendicular to the planar interface between two binary fluids, with no impedance mismatch ͑␣ =0͒.The amplitude of the acoustic waves was selected to be comparable to the density change at the interface for = 0.001.This is shown in Fig. 2 where the pulse is shown propagating through the interface at selected times.There was no evidence of either a reflection at the interface or a change in the amplitude of the pulse.When there is an impedance mismatch at the interface we would expect an acoustic signal to be partially transmitted and partially reflected.This is shown for the same pulse in Figs. 3 and 4  in density at the interface is significantly smaller than the acoustic pulse.Figure 3 shows a plane interface between the two fluids, and Fig. 4 shows a circular bubble of fluid 2 ͑⌬ Ͼ 0͒ where the density is measured through the center of the bubble.In both cases there is clear evidence of a reflected wave at each interface and a change in the amplitude of the transmitted wave.In Fig. 4 a change in the width of the pulse is evident inside the bubble.This is expected due to the change in the speed of sound and is also evident in Fig. 3.It is also clear from Fig. 4 that any density difference between the inside and outside of the bubble, due to surface tension, is negligible for the interfacial energy applied here.
The simulation depicted in Fig. 3 was repeated for a range of ␣ and the amplitude of the incident, transmitted, and reflected components measured.From this the transmission coefficient C T and the reflection coefficient C R were calculated.These are shown in Fig. 5 and compared to the theoretical expressions derived in the Appendix.Good agreement between theory and simulation is observed over the range of ␣.
The model was then used to investigate an acoustic wave propagating through a bubble with a higher speed of sound.The bubble was simulated with a radius of 50 lattice units and ␣ was set to −0.4, giving a speed of sound of 0.89 inside the bubble and 0.63 outside.Plane acoustic waves were simulated moving from left to right through the fluid with a frequency of 0.1124.This gives a wavelength inside the bubble approximately equal to the radius.The resulting density field is depicted in Fig. 6 where the position of the bubble, determined by ⌬ = 0, is also shown.In front of the bubble the acoustic signal is dominated by the incident wave with a small component due to the reflected signal.Inside the bubble the transmitted component of the incident wave is also dominant; however, a reflected component from the second boundary can be observed.This can be seen in more detail in Fig. 7 which shows the density variation inside the bubble along a line through its center ͑y = 250͒.The position of the wave crest and trough were recorded at every time step over a wave period and are depicted by the cross in Fig. 7 to indicate the envelope of the acoustic signal.Also shown is a sample wave which sits inside this envelope.If no reflected component was present, the envelope would show a constant value for the density ͑other than a small decrease with increasing x due to the viscosity of the fluid͒.The envelope depicted in Fig. 7 represents the steady state of the system and is independent of the period selected and the initial phase of the wave.A similar effect was observed in front of the bubble but with a smaller standing-wave component.Behind the bubble there is an interference pattern.This is due to the imaging effect of the circular bubble.Directly behind the bubble an approximately plane wave was ob-  served to propagate away from the bubble with no standingwave component, but with a reduced amplitude relative to the incident wave.In Fig. 7 there are a number of points where the shape of the envelope is not completely smoothfor example, the top envelope at x = 0.These small deviations are due to the initialization of the bubble varying slightly from the final state.In particular, the system was initialized without the drop in density across the interface.

IV. CONCLUSIONS
A numerical technique based on an LBM for a binary fluid with species-dependent impedance was been investigated.Two immiscible fluids were simulated, and acoustic waves and pulses were simulated propagating across the interface.The speed of sound was varied between 0.1 and 1, and the reflection and transmission coefficients were measured.Excellent agreement was observed between the simulated coefficients and theory.The model was then applied to investigate a plane-wave propagating through a bubble with a higher speed of sound than the surrounding fluid.A partial standing wave inside the bubble and the interference pattern produced from the bubble lens were both observed.

APPENDIX: REFLECTION AND TRANSMISSION COEFFICIENTS
Consider an acoustic wave propagating from fluid 1 with speed of sound c 1 and acoustic impedance ⑀ 1 to fluid 2 with speed of sound c 2 and acoustic impedance ⑀ 2 .At the boundary i + r = t ͑A1͒ and where the subscripts i, r, and t refer to the incident, reflected, and transmitted components, respectfully.Equations ͑A1͒ and ͑A2͒ can be expressed in terms of the pressure using p = c s 2 and p = ±⑀u, where c s is the speed of sound in the particular fluid, ⑀ is the acoustic impedance, and ± refers to the direction of wave propagation.This gives FIG.1.Variation of the density and the order parameter ⌬ across a plane interface at x = 50.5 for = 0.001 and = 0.0001.