A computational framework based on immersed boundary-lattice Boltzmann coupling for ultrasound-excited encapsulated microbubbles

1 Introduction

Encapsulated microbubbles (EMBs) have attracted considerable attention in biomedical ultrasound because of their important roles in diagnostic imaging, targeted drug delivery, sonoporation, mechanotransduction, and therapeutic ultrasound applications (de Jong et al., 1991; Lindner, 2004; Hernot and Klibanov, 2008; Ferrara et al., 2007). When EMBs are used to enhance ultrasound imaging contrast, they are commonly referred to as ultrasound contrast agents. Under acoustic excitation, EMBs exhibit complex oscillatory behavior governed by the interaction between the surrounding fluid, the encapsulating viscoelastic shell, and the imposed ultrasound field. Accurate numerical modeling of EMB dynamics is therefore essential for understanding ultrasound-mediated biophysical phenomena and for improving the design and optimization of ultrasound-based biomedical technologies.

A substantial portion of the existing theoretical and numerical studies on ultrasound-driven EMBs has been based on Rayleigh-Plesset-type (RP-type) formulations. Classical RP-type models, including those proposed by the researchers (de Jong et al., 1992; Church, 1995; Hoff et al., 2000; Chatterjee and Sarkar, 2003; Marmottant et al., 2005; Sarkar et al., 2005; Tsiglifis and Pelekasis, 2008; Paul et al., 2010; Chabouh et al., 2021), extend the original RP equation (RPE) by incorporating elasticity, viscous damping, buckling and rupture, and compressibility and anisotropy of the shell. These models have significantly advanced the understanding of coated bubble acoustics and remain widely used because of their computational efficiency and analytical tractability. However, RP-type approaches fundamentally rely on the assumption of spherical symmetry and therefore are primarily restricted to predicting the radial oscillations of EMBs. Consequently, they are unable to fully resolve nonspherical deformations, interfacial instabilities, jet formation, or complex fluid-structure interactions that may arise under strong acoustic excitation or in the presence of nearby boundaries.

To address these limitations, several studies have investigated nonspherical EMB dynamics using boundary element or boundary integral formulations. Hsiao and Chahine (Hsiao and Chahine, 2013) examined ultrasound-induced breakup of viscous shell microbubbles using a simplified shell representation, while Tsigklifis and Pelekasis (Tsigklifis and Pelekasis, 2013) studied transient breakup and saturation phenomena of insonated contrast agents. Furthermore, Wang et al. (Wang et al., 2015) developed a three-dimensional (3-D) boundary integral model for nonspherical contrast agent dynamics, Mobadersany and Sarkar (Mobadersany and Sarkar, 2018) investigated EMB collapse and jet formation near a membrane during sonoporation, Wu et al. (Wu et al., 2024) developed a 3-D nonlinear model for describing the nonspherical deformation of an EMB under acoustic excitation, and Furukawa et al. (Furukawa et al., 2025) numerically investigated the nonspherical dynamics of thin-shelled micron-sized bubbles exposed to ultrasound. These investigations demonstrate the importance of nonspherical interfacial dynamics in acoustically-driven EMBs. Nevertheless, boundary element and boundary integral approaches commonly assume potential flow in the surrounding liquid and therefore neglect liquid viscosity. Although this assumption can be reasonable in some inertial cavitation regimes, viscous effects become increasingly important in many biomedical ultrasound applications. In particular, viscous damping of the surrounding medium may substantially influence EMB oscillations under ultrasound excitation (Lauterborn and Kurz, 2010; Nazari-Mahroo et al., 2018; Qin et al., 2023). Moreover, when a bubble oscillates or collapses near a stationary rigid wall, such as during sonoporation and ultrasound-mediated drug delivery, solving the surrounding viscous flow with a no-slip wall condition naturally incorporates near-wall viscous effects, including velocity gradients and wall shear stresses (Popinet and Zaleski, 2002; Minsier et al., 2009). Consequently, neglecting surrounding liquid viscosity may limit the physical realism of potential-flow-based formulations in biologically relevant environments.

Fully resolved computational fluid dynamics (CFD) approaches offer a promising alternative for overcoming these limitations because they can directly solve the surrounding multiphase flow while resolving the bubble-fluid interaction without imposing spherical symmetry or potential flow assumptions. Among available CFD techniques, the lattice Boltzmann method (LBM) has emerged as an attractive framework for simulating multiphase and interfacial flows, owing to its mesoscopic formulation, parallel efficiency, and suitability for handling complex moving interfaces (Chen and Doolen, 1998; Krüger et al., 2017). In addition, the immersed boundary (IB) method provides an effective approach for representing deformable interfaces and fluid-structure interactions without requiring body-fitted meshes or interface reconstruction procedures (Mittal and Iaccarino, 2005).

In the present study, a computational framework based on IB-LB coupling is developed for simulating ultrasound-excited EMBs. The proposed framework combines an axisymmetric multicomponent multiphase (MCMP) LB flow solver with an IB representation of the microbubble’s viscoelastic shell, thereby enabling fully resolved two-way coupling between the surrounding fluid domains and the encapsulating shell within an acoustically-driven environment. Unlike RP-type reduced-order models, the present methodology directly resolves the surrounding flow field and interfacial dynamics while naturally incorporating viscous effects in both the liquid and shell regions. Furthermore, unlike boundary element or boundary integral formulations, the proposed framework does not rely on inviscid flow assumptions and is therefore more suitable for investigating EMB dynamics in confined and biologically relevant environments where viscous effects become significant.

The shell viscous effects are incorporated via the Boussinesq-Scriven constitutive law and the shell elasticity is integrated using an exponential elasticity model (EEM), allowing systematic investigation of shell property effects on acoustic bubble dynamics. Acoustic forcing is introduced through a time-dependent far-field pressure condition representative of diagnostic and therapeutic ultrasound regimes. The developed framework is validated against modified RP solutions for spherical EMB oscillations and is subsequently used to investigate nonlinear acoustic regimes and the influence of shell viscoelastic properties on the EMB response.

Overall, the proposed solver aims to bridge the gap between classical reduced-order bubble models and fully resolved multiphase CFD methodologies. By combining physical fidelity, numerical robustness, and computational flexibility, the present method provides a general computational platform for future investigations of nonspherical EMB dynamics, bubble-boundary interactions, ultrasound-mediated drug delivery, mechanotransduction, and other biomedical ultrasound phenomena involving complex fluid-structure coupling.

2 Methods

2.1 Multicomponent multiphase lattice Boltzmann method

Within the LB framework, fluid motion is described at a mesoscopic level by representing the fluid as fictitious particles that propagate along discrete lattice directions and undergo collisions at lattice sites (Mohamad, 2011). The number of allowed propagation directions is determined by the selected lattice stencil. In the present study, the D2Q9 stencil is adopted, meaning that the flow is treated as two-dimensional (2-D) and particles can move along nine discrete velocity directions (see Figure 1). This stencil is chosen because it provides a suitable compromise between computational efficiency, isotropy, and the accurate recovery of the Navier-Stokes and continuity equations (Succi, 2001). The discrete lattice velocity vectors c k are defined as:where c = h / Δ t, with h being the lattice spacing and Δ t the lattice time step size.

c k = 0 , 0 , k = 0 , cos π k − 1 2 , sin π k − 1 2 c , k = 1 t o 4 , cos π k − 9 / 2 2 , sin π k − 9 / 2 2 2 c , k = 5 t o 8 , (1)

The fundamental variable in the LBM is the particle distribution function f k x , t, which represents the probability of finding a particle at position x and time t moving with velocity c k (Krüger et al., 2017). At each time step, the distributions undergo two successive processes: collision and streaming. During the collision step, the populations are locally redistributed, while in the streaming step they propagate to neighboring lattice nodes (Qian et al., 1992).

In this work, the Shan-Chen (Shan and Chen, 1993; Shan, 2010) MCMP LBM is employed to model the flow. This approach introduces interaction forces between fluid components to enable phase separation. For each component i, a separate distribution function f k ( i ) is defined, and its evolution is governed by the LB equation:where the standard Bhatnagar-Gross-Krook collision operator is given by:Here, in Equation 3, τ ( i ) is the relaxation time associated with fluid component i, and S mass k ( i ) is a source term arising from the axisymmetric formulation, which will be defined later. The equilibrium distribution functions f k eq ( i ) are expressed as a low-Mach-number approximation of the Maxwell-Boltzmann distribution:where u eq ( i ) is the equilibrium velocity of component i, c s is the LBM’s sound speed, and w k are direction-dependent weights, with values w 0 = 4 / 9 , w 1 − 4 = 1 / 9, and w 5 − 8 = 1 / 36 (He and Luo, 1997). The Chapman-Enskog analysis relates the kinematic viscosity ν ( i ) to the relaxation time τ ( i ) of the i-th component as:with c s = c / 3. Several velocity definitions are required. The velocity u ( i ) of each component is:with ρ ( i ) being the i-th component density. The barycentric velocity u b of the mixture is:This is the physical velocity satisfying the Navier-Stokes equations (Shan and Doolen, 1995), and ρ denotes the mixture density. Also, the Shan-Chen forcing scheme (Martys and Chen, 1996) defines the equilibrium velocity as:where the common velocity is defined by the weighted average:Since the total momentum of particles of all components should be conserved by the collision operator at each lattice site (Krüger et al., 2017; Shan and Doolen, 1995), the relaxation time τ ( i ) is used in the common velocity definition given by Equation 9. The total force acting on each component is:where F SC ( i ) represents Shan-Chen surface tension forces, F mom ( i ) accounts for axisymmetric corrections, and F IB arises from the IB, i.e., the EMB surface. Because external, i.e., non-surface-tension, forces are distributed to the components according to their concentration, the density fraction ρ ( i ) / ρ is multiplied by such forces in Equation 10 (Krüger et al., 2017). An axisymmetric flow is considered in this study and the cylindrical coordinates ( r , z , θ ), wherein the flow in the r-z plane is axisymmetric with respect to the z-axis, i.e., there are no variations with respect to the azimuthal angle θ, are used. Furthermore, two immiscible fluid components are considered in the present investigation, with i = 1 and i = 2 reflecting the liquid and gas components, respectively.

f k i x + c k Δ t , t + Δ t = f k i x , t + Ω k i x , t + Δ t S mass k i x , t , (2)

Ω k i x , t = − Δ t τ i f k i x , t − f k eq i x , t . (3)

f k eq i = w k ρ i 1 + u eq i ⋅ c k c s 2 + u eq i ⋅ c k 2 2 c s 4 − u eq i ⋅ u eq i 2 c s 2 , (4)

ν i = c s 2 τ i − Δ t 2 , (5)

u i = 1 ρ i ∑ k f k i c k ; ρ i = ∑ k f k i , (6)

u b = 1 ρ ∑ i ∑ k f k i c k + F i Δ t 2 ; ρ = ∑ i ρ i . (7)

u eq i = u ′ + τ i F i ρ i , (8)

u ′ = ∑ i ρ i u i τ i ∑ i ρ i τ i . (9)

F i = F SC i + ρ i ρ F mom i + F IB , (10)

Surface tension is modeled through the intra- and inter-molecular forces given by Equations 11a,b (Garousi, 2025): where F 11 and F 22 are the intra-molecular forces, and F 12 and F 21 are the inter-molecular forces of a Cartesian 2-D flow. Their continuum forms are (Benzi et al., 2006; Lycett-Brown and Luo, 2015; Krüger et al., 2017): The parameter G denotes the interaction strength, ψ and φ represent intra- and inter-molecular pseudopotentials, respectively, and the gradient ∇ and Laplacian ∇ 2 operators are given by ∇ = ( ∂ r , ∂ z ) and ∇ 2 = ∂ r r + ∂ z z. The additional terms ensuring proper 3-D surface tension representation are (Garousi, 2025): Inter-molecular pseudopotentials given by Equations 14a,b are defined as (Yu et al., 2007): where the constants in the lattice system are set as a 0 = 0.005, ρ 10 = − 0.0008 / l o g ( a 0 ), and ρ 20 = 0.0003 (Yu, 2009). In this work, the gas component is considered to be an ideal fluid, thus G 22 = 0. Additionally, the inter-molecular interaction strengths are assumed to be G 12 = G 21 = 0.0001 (Moslemi et al., 2023). The liquid component uses a non-ideal equation of state (EOS) and its intra-molecular pseudopotential is computed as (Krüger et al., 2017):while the Carnahan-Starling EOS is used to give (Huang et al., 2015):In Equation 16, T is the temperature in the lattice system, and the lattice values of the parameters in the Carnahan-Starling EOS are a = 1, b = 4, and R = 1 (Yuan and Schaefer, 2006). Also, G 11 is set to ± 1 to ensure that the radicand in Equation 15 always remains positive.

F SC 1 = F 11 + F 12 + F 1 1 axi + F 1 2 axi , (11a)

F SC 2 = F 22 + F 21 + F 2 2 axi + F 2 1 axi , (11b)

F 11 = − ψ 1 G 11 c s 2 Δ t 2 ∇ ψ 1 + c s 4 Δ t 4 2 ∇ ∇ 2 ψ 1 , (12a)

F 12 = − φ 1 G 12 c s 2 Δ t 2 ∇ φ 2 + c s 4 Δ t 4 2 ∇ ∇ 2 φ 2 , (12b)

F 22 = − ψ 2 G 22 c s 2 Δ t 2 ∇ ψ 2 + c s 4 Δ t 4 2 ∇ ∇ 2 ψ 2 , (12c)

F 21 = − φ 2 G 21 c s 2 Δ t 2 ∇ φ 1 + c s 4 Δ t 4 2 ∇ ∇ 2 φ 1 . (12d)

F 1 1 axi = − ψ 1 G 11 c s 4 Δ t 4 2 ∇ ∂ r ψ 1 r , (13a)

F 1 2 axi = − φ 1 G 12 c s 4 Δ t 4 2 ∇ ∂ r φ 2 r , (13b)

F 2 2 axi = − ψ 2 G 22 c s 4 Δ t 4 2 ∇ ∂ r ψ 2 r , (13c)

F 2 1 axi = − φ 2 G 21 c s 4 Δ t 4 2 ∇ ∂ r φ 1 r . (13d)

φ 1 = 1 − e x p − ρ 2 / ρ 20 , (14a)

φ 2 = a 0 − e x p − ρ 1 / ρ 10 , (14b)

ψ 1 = 2 c s 2 Δ t 2 G 11 P EOS − ρ 1 c s 2 , (15)

P EOS = ρ 1 R T 1 + b ρ 1 4 + b ρ 1 4 2 − b ρ 1 4 3 1 − b ρ 1 4 3 − a ρ 1 2 , (16)

The definitions of the mass source term S mass k ( i ) in Equation 2 and of the additional momentum force F mom ( i ) in Equation 10 for an axisymmetric flow are according to the Srivastava et al.’s (Srivastava et al., 2013) work. Assuming zero azimuthal velocity u θ ( i ), i.e., no swirl, the final forms of S mass k ( i ) (as given by Equation 17a), and of the r- and z-components of the vector F mom ( i ) (as defined by Equations 17b,c) are: Furthermore, the derivatives of a scalar field χ, appearing in Equation 12a through Equations 12d, 13a through Equation 13d, and in Equations 17b,c, are evaluated using the isotropic fifth-order finite differences expressed by Equation 18a through Equation 18e as (Srivastava et al., 2013): where c k r and c k z are respectively the r- and z-components of the lattice velocity vector given by Equation 1.

S mass k i = − w k ρ i u r i r , (17a)

F mom r i = 1 r 2 ρ i ν i ∂ r u r i − u r i r − ρ i u r i 2 , (17b)

F mom z i = 1 r ρ i ν i ∂ r u z i + ∂ z u r i − ρ i u r i u z i . (17c)

∂ r χ = 1 36 ∑ k = 0 8 8 χ x + c k Δ t − χ x + 2 c k Δ t c k r , (18a)

∂ r r χ = 1 36 ∑ k = 0 8 8 ∂ r χ x + c k Δ t − ∂ r χ x + 2 c k Δ t c k r , (18b)

∂ z χ = 1 36 ∑ k = 0 8 8 χ x + c k Δ t − χ x + 2 c k Δ t c k z , (18c)

∂ z z χ = 1 36 ∑ k = 0 8 8 ∂ z χ x + c k Δ t − ∂ z χ x + 2 c k Δ t c k z , (18d)

∂ r z χ = 1 36 ∑ k = 0 8 8 ∂ z χ x + c k Δ t − ∂ z χ x + 2 c k Δ t c k r , (18e)

2.2 Immersed boundary method

In the IB method, the interface is represented by discrete Lagrangian markers r j, while the fluid field is solved on an Eulerian grid. Figure 2 shows the combination of these two grid systems. The coupling between the interface and fluid occurs through two key steps: velocity interpolation and force spreading. During velocity interpolation (Equation 19a), the velocity of each marker j is obtained from nearby fluid velocities, whereas force spreading (Equation 19b) distributes the forces from the markers back to the fluid as a body force term F IB. The governing discrete relations are (Peskin, 2002): where g j represents the EMB shell force acting on each marker, determined from an appropriate viscoelastic constitutive formulation. The discrete Dirac delta function D is factorized as D = ϕ ( r ) ϕ ( z ) / h 2, with ϕ defining the interpolation-spreading stencil.

r ̇ j = ∑ x h 2 u b x , t D r j − x , (19a)

F IB x , t = ∑ j g j t D r j − x , (19b)

2.3 Derivation of the EMB shell force

The derivation begins with the interfacial linear momentum transport equation, expressed as (Edwards et al., 1991):where ρ s denotes the shell density, D / D t is the material derivative, and V s represents the shell velocity vector. The tensor P s is the shell pressure dyadic, while F s corresponds to the excess shell force density vector (e.g., gravitational or electromagnetic forces). The vector n is the unit normal to the interface, and P ̄ denotes the mean pressure dyadic of the bulk flow. For a generic tensor field G, the notation G is defined as G = G I − G II, where G I and G II are the values of the tensor on opposite sides of the interface. Owing to the extremely small thickness of the interfacial transition layer, which is especially relevant for lipid-coated microbubbles (Klibanov, 2006), the shell density can reasonably be neglected, i.e., ρ s ≈ 0. Consequently, Equation 20 simplifies to:The shell pressure dyadic P s can generally be decomposed into isotropic and deviatoric parts according to:where σ is the effective surface tension, T s denotes the shell stress dyadic, and I s is the surface identity dyadic. Following the Boussinesq-Scriven constitutive relation for Newtonian interfaces, the shell stress dyadic takes the form (Scriven, 1960):where μ s and κ s represent the shell shear viscosity and shell dilatational viscosity, respectively. The Boussinesq-Scriven constitutive model is a classical continuum mechanics framework for modeling viscous interfacial rheology. The model introduces interfacial shear and dilatational viscosities and therefore enables the shell to resist tangential deformation and surface area variation in a physically meaningful manner. Since EMB shells behave as rheologically complex interfaces coated with lipids, proteins, polymers, or surfactants, incorporation of interfacial viscous stresses is important for accurately capturing dissipative effects during oscillation. The Boussinesq-Scriven formulation has widely been used in interfacial fluid mechanics studies involving surfactant-laden interfaces, deformable droplets, liquid films, and viscoelastic interfaces (Yu and Zhou, 2011; Schwalbe et al., 2011; Martínez-Calvo and Sevilla, 2018; Cunha et al., 2021; Jaensson et al., 2021; Cunha et al., 2023), because it provides a continuum description of interfacial viscous stresses while remaining mathematically tractable for numerical implementation. The shell rate-of-deformation dyadic D s is defined as:In Equation 23, the colon operator ‘:’ indicates the scalar contraction between the dyadics I s and D s, while in Equation 24, the superscript ‘Tr’ denotes transpose of the dyadic ∇ V s. Substituting Equation 24 into Equation 23, then substituting the resulting T s into Equation 22, and finally inserting the resulting P s into Equation 21, yield the stress boundary condition (BC) for a Newtonian interface:where B s = − ∇ n is the surface curvature dyadic and H s = − 0.5 ∇ ⋅ n is the mean surface curvature. Equation 25 may be decomposed into normal and tangential stress components relative to the EMB surface:The decomposition given by Equation 26 consists of the following normal and tangential components: Equations 27a,b correspond to the normal and tangential stress BCs for a Newtonian interface. To account for shell elasticity, the effective surface tension is defined as (Paul et al., 2010):where σ 0 is the reference surface tension and the function E depends on the shell elasticity E s, a dimensionless shell parameter α s, the instantaneous EMB surface area A, and the unstrained equilibrium surface area A E. Using the EEM proposed by Paul et al. (Paul et al., 2010), the function E is written as E = E s ⁡ exp ( − α s δ A A E ) δ A A E, where the area change is δ A = A − A E. Also, the unstrained equilibrium surface area A E is defined as (Paul et al., 2010):with R E and R 0 in Equation 29 being the unstrained equilibrium and equilibrium radii of the EMB, respectively. The EEM is chosen because of its capability to capture the experimentally observed strain-softening behavior of ultrasound contrast agent shells under finite-area deformation. Assuming zero excess shell force density, i.e., F s = 0, and taking both the reference surface tension and shell elasticity to be constant, the effective surface tension given by Equation 28 is substituted into Equations 27a,b. After multiplying both sides of these equations by the Voronoi area A j associated with Lagrangian marker j, and retaining the terms involving the shell parameters μ s, κ s, and E s, the normal g n j and tangential g t j forces exerted by the EMB’s viscoelastic shell at marker j position become (Garousi, 2025): The total shell force acting at marker j is therefore:

ρ s D V s D t − ∇ ⋅ P s − F s = n ⋅ P ̄ , (20)

− n ⋅ P ̄ = F s + ∇ ⋅ P s . (21)

P s = I s σ + T s , (22)

T s = κ s − μ s I s : D s I s + 2 μ s D s , (23)

D s = 1 2 ∇ V s ⋅ I s + I s ⋅ ∇ V s Tr . (24)

− n ⋅ P ̄ = F s + 2 H s σ n + ∇ σ + κ s + μ s ∇ ∇ ⋅ V s + 2 μ s n B s − 2 H s I s : ∇ V s + 2 H s n κ s + μ s ∇ ⋅ V s + μ s n × ∇ ∇ × V s ⋅ n − 2 μ s B s − 2 H s I s ⋅ ∇ V s ⋅ n , (25)

n ⋅ P ̄ = P ̄ n + P ̄ t . (26)

P ̄ n = n ⋅ P ̄ ⋅ n n = − F s ⋅ n n − 2 H s σ n − 2 μ s n B s − 2 H s I s : ∇ V s − 2 H s n κ s + μ s ∇ ⋅ V s , (27a)

P ̄ t = n ⋅ P ̄ ⋅ I s = − F s ⋅ I s − ∇ σ − κ s + μ s ∇ ∇ ⋅ V s − μ s n × ∇ ∇ × V s ⋅ n + 2 μ s B s − 2 H s I s ⋅ ∇ V s ⋅ n . (27b)

σ = σ 0 + E E s , α s , A , A E , σ > 0 0 , σ ≤ 0 , (28)

A E = 4 π R E 2 = 4 π R 0 1 + 1 − 1 + 4 σ 0 α s / E s 2 α s − 1 / 2 2 , (29)

g n j = − 2 H j E − 2 μ s B j − 2 H j I j : ∇ r ̇ j − 2 H j κ s + μ s ∇ ⋅ r ̇ j A j n j , σ > 0 − 2 μ s B j − 2 H j I j : ∇ r ̇ j − 2 H j κ s + μ s ∇ ⋅ r ̇ j A j n j , σ ≤ 0 , (30a)

g t j = − A j κ s + μ s ∇ ∇ ⋅ r ̇ j − A j μ s n j × ∇ ∇ × r ̇ j ⋅ n j + 2 μ s A j B j − 2 H j I j ⋅ ∇ r ̇ j ⋅ n j . (30b)

g j = g n j + g t j . (31)

To evaluate the unit normal vector at marker j position, denoted by n j, the Lagrangian marker position vector r j is parameterized in terms of the arc length ξ according to Equation 32:As given by Equations 33a,b, the functions R j and Z j are represented using cubic spline interpolants over each subinterval ξ j , ξ j + 1 for each j = 0,1 , … , N − 1: The first marker ( j = 0 ) is located at the north pole of the EMB, while j = N represents the last marker. The tangent vector at marker j position is obtained from the derivative of the parametric representation according to Equation 34:As expressed by Equation 35, a normal vector is then generated by rotating the tangent vector by 9 0 0, giving:and the corresponding unit normal vector is (Xie et al., 2012):In the present work, the arc length ξ j at each marker j position is evaluated at every time step using the Kucera’s iterative procedure (Kucera and Blake, 1990; Best, 1991; Pearson, 2003).

r j = R j ξ = ξ j e r + Z j ξ = ξ j e z , j = 0,1 , … , N . (32)

R j ξ = r 1 j + r 2 j ξ − ξ j + r 3 j ξ − ξ j 2 + r 4 j ξ − ξ j 3 , (33a)

Z j ξ = z 1 j + z 2 j ξ − ξ j + z 3 j ξ − ξ j 2 + z 4 j ξ − ξ j 3 . (33b)

t j = d r j d ξ | ξ = ξ j = d R j d ξ | ξ = ξ j e r + d Z j d ξ | ξ = ξ j e z . (34)

n o j = − d Z j d ξ | ξ = ξ j e r + d R j d ξ | ξ = ξ j e z , (35)

n j = n o j n o j . (36)

2.4 Multicomponent multiphase IB-LBM algorithm

To simulate the dynamics of EMBs, the coupling between the MCMP LBM and the IB method is implemented through the following iterative procedure:

• Starting from the current IB configuration r j:
• Apply the Kucera’s iterative method to determine the local unit normal vector n j on the EMB surface according to Equation 36.
• Evaluate the total EMB shell force g j using Equation 31, where the normal and tangential contributions are computed from the pair of Equations 30a,b.
• Spread the total EMB shell force g j onto the surrounding fluid domain through Equation 19b in order to compute the Eulerian body force associated with the IB, namely, F IB.
• Determine the total body force acting on each fluid component i using Equation 10.
• Compute the bare component velocity u ( i ) from Equation 6, the common velocity u ′ from Equation 9, and the equilibrium velocity u eq ( i ) via Equation 8.
• Evaluate the equilibrium distribution functions f k eq ( i ) using Equation 4, and then solve the LB Equation 2, to advance the populations f k ( i ).
• Calculate the fluid mixture’s barycentric velocity u b using Equation 7.
• Interpolate the Lagrangian grid velocity r ̇ j from the barycentric velocity u b through Equation 19a.
• Employ a desired advection scheme, e.g., the explicit forward Euler method, to advect the Lagrangian markers and obtain the updated IB configuration.
• Proceed to the next time step by returning to Step 1.

3 Results

3.1 Young-Laplace test for a 3-D bubble

The equilibrium configuration of a gas-filled, free (uncoated) bubble suspended in a liquid serves as a suitable benchmark problem for validating an MCMP flow solver. Similar Young-Laplace (Y-L) validation studies have also been conducted in previous investigations (Lallemand et al., 2007; Xie et al., 2012; Srivastava et al., 2013; Xie et al., 2015; Su et al., 2017). The EMB shell force is set to zero in this part of the study, i.e., g j = 0. The center of a spherical bubble with an initially non-equilibrium radius R init is positioned at ( r , z ) = ( 0 , L z / 2 ) inside an axisymmetric computational domain, as illustrated in Figure 2. To track the IB representing the bubble surface, N mark Lagrangian markers are initially distributed along the interface with spacing h between neighboring markers, following the recommendation of Krüger et al. (Krüger et al., 2017). Given the initial bubble radius R init, the initial arc length of the axisymmetric bubble is calculated as L 0 = π R init. The total number of Lagrangian markers then equals N mark = ( L 0 / / h ) + 1, where the double slash ‘//’ signifies the integer division. The first marker is positioned at the north pole, and the remaining markers are uniformly distributed toward the south pole in a clockwise manner. The coordinates of the marker located at the north pole are used to compute the instantaneous bubble radius as R = ( r j = 0 ) 2 + ( z j = 0 − L z / 2 ) 2, where r j = 0 and z j = 0 are the r- and z-components of the north-pole marker position vector r j = 0. The Lagrangian markers are advanced in time using the first-order forward Euler method according to Equation 37:For the evaluation of the discretized Dirac delta function, the stencil ϕ 2 defined by Equation 38 as (Peskin, 2002):is employed, where x denotes either of the coordinate axes r or z, and the index ‘2’ specifies the count of Eulerian nodes used along each axis for interpolation and spreading.

r j t + Δ t = r j t + r ̇ j t Δ t . (37)

ϕ 2 x = 1 − | x | , 0 ≤ | x | ≤ h , 0 , | x | > h , (38)

Any length L in the lattice system can be related to its physical counterpart L ph as L / h = L ph / h ph, where h ph is the physical grid size. Likewise, lattice time t is related to physical time t ph through t / Δ t = t ph / Δ t ph, where Δ t ph is the physical time step size. In the present investigation, the lattice time step size is set to Δ t = 1 (Krüger et al., 2017). Within the computational domain shown in Figure 2, the left boundary corresponds to the symmetry axis, periodic BCs are imposed on the top and bottom boundaries, and a free-slip BC is prescribed on the right boundary. The dimensions of the domain are chosen as L r = 100 and L z = 200 in lattice units. The ratio of the relaxation times of the two fluid components is selected as τ ( 1 ) / τ ( 2 ) = 0.8436 to reproduce the water-to-perfluorobutane kinematic viscosity ratio. This value is obtained from Equation 5 together with the kinematic viscosities of the two fluids. Perfluorobutane ( C 4 F 10 ) is a typical gas agent used to fill commercial microbubbles, e.g., Sonazoid (Kudo et al., 2008; Moriyasu and Itoh, 2009). It should be emphasized that in MCMP LBM-based solvers, both fluid components are present at every Eulerian node in the computational domain. Densities in the lattice system, i.e., density of component 1 in the gas region ρ g ( 1 ), density of component 2 in the gas region ρ g ( 2 ), density of component 1 in the liquid region ρ l ( 1 ), and density of component 2 in the liquid region ρ l ( 2 ), are initialized as ρ g ( 1 ) = 0.000108, ρ g ( 2 ) = 0.00023, ρ l ( 1 ) = 0.41556, and ρ l ( 2 ) = 4.58 × 1 0 − 8 (Yu, 2009). To initialize the populations f k ( i ), the Mei et al.’s (Mei et al., 2006) approach is adopted, generating consistent non-equilibrium distribution functions, thereby producing stable initial populations. Following initialization of the densities and populations, the MCMP fluid system is allowed to evolve until equilibrium is reached. Once equilibrium is established, the bubble equilibrium radius R 0 is measured. Also, the following relation (Chen et al., 2015):is used to calculate the equilibrium pressures inside the bubble P g and outside it P l. In Equation 39, i and q refer to different fluid components. Using the measured values of R 0, P g, and P l, the reference surface tension is computed from the Y-L relation for a 3-D bubble as σ 0 = R 0 P g − P l / 2. The equilibrium mixture densities in the liquid and gas regions, i.e., ρ l and ρ g, are also measured, enabling calculation of the liquid-to-gas density ratio.

P = ρ c s 2 + 1 2 c s 2 ∑ G i i ψ i 2 + 1 2 c s 2 ∑ i ≠ q G i q φ i φ q , (39)

Table 1 shows the sensitivity of the equilibrium pressure jump across the interface, i.e., P g − P l, to the lattice spacing h. It is seen that the relative difference between the pressure jump values for h = 0.5 and h = 0.25 is 0.57%. Therefore, in the rest of the study, the lattice spacing is set to h = 0.25 to ensure a sufficiently fine grid. It should be noted that T r denotes reduced temperature, defined as the ratio of the lattice temperature T to the lattice critical temperature T crit, i.e., T r = T / T crit. The critical temperature is set to T crit = 0.0943 (Huang et al., 2015), below which the separation of liquid and gas phases is possible. Figure 3 presents the temporal evolution of the normalized bubble radius until equilibrium is attained at t = 31,025. For the same case, Figure 4 displays the normalized mixture density distribution at z = 100, corresponding to the axial location of the bubble center. The equilibrium analysis is repeated for bubbles with different initial radii. Once equilibrium is achieved in each case, the interfacial pressure jump, i.e., Δ P = P g − P l, and the inverse equilibrium radius ( 1 / R 0 ) are calculated. Figure 5 plots Δ P as a function of 1 / R 0. The resulting linear relationship confirms that the present MCMP IB-LBM correctly reproduces the Y-L law for spherical bubbles in static equilibrium. The influence of the reduced temperature T r on the achievable liquid-to-gas density ratio is also investigated. Several bubbles with identical initial radii are allowed to reach equilibrium at different reduced temperatures, after which the equilibrium density ratios are evaluated. Table 2 summarizes the results. Density ratios ranging from 70.7551 to 1387.0 are achieved for T r values between 0.7 and 0.55. Since realistic multiphase systems such as air bubbles in water typically exhibit density ratios on the order of 1000, these results demonstrate that the present MCMP flow solver is capable of handling physically relevant density contrasts. In addition, decreasing the reduced temperature is found to increase the achievable density ratio, which agrees with trends reported previously in the literature (He et al., 2020; Ezzatneshan and Vaseghnia, 2021).

TABLE 1. Sensitivity of the interfacial pressure jump, i.e., P g − P l, to the lattice spacing h when R init = 57 and T r = 0.6.

TABLE 2. Effect of reduced temperature T r on the equilibrium liquid-to-gas density ratio ( ρ l / ρ g ) for a spherical bubble when R init = 55.

3.2 Radial oscillations of an EMB exposed to ultrasound

CFD results of capturing the radial oscillations of an EMB submerged in a liquid and subject to ultrasound, i.e., acoustic forcing, using the developed MCMP IB-LBM are presented here. After conducting the Y-L test for a 3-D uncoated bubble and achieving an equilibrium MCMP flow field, the EMB shell force g j, defined by Equation 31, comes into play to enable the modeling of an EMB dynamics. The left boundary of the computational domain shown in Figure 2 remains the symmetry axis, while the Zou-He (Zou and He, 1997) pressure BC is applied at the top, bottom, and right boundaries. The Zou-He BC is based on the bounce-back of the non-equilibrium portion of the populations, which corresponds to the deviation between the actual and equilibrium populations and contains information related to viscous stresses, momentum transport, and higher-order effects (Aidun and Clausen, 2010; Krüger et al., 2017). At a boundary node, the populations whose discrete velocities point inward after streaming are unknown. The Zou-He method reconstructs these unknown populations by prescribing the pressure at the boundary and ensuring momentum consistency through reflection of the non-equilibrium parts of opposite populations. Further implementation details can be found in the literature (Zou and He, 1997; Feiger et al., 2019).

The EMB is initially in static equilibrium and then is exposed to a time-varying acoustic pressure P a defined as P a ( t ) = P 0 − Δ P a ⁡ sin ( 2 π F t ), where F is the frequency of oscillations. The pressure variation amplitude is given by Δ P a = ϵ P 0, with ϵ being the dimensionless amplitude and P 0 the hydrostatic liquid pressure. To validate the numerical predictions, the IB-LBM results are compared with solutions of a modified RPE. Derived from the continuity and radial momentum equations in the spherical coordinates assuming there is no mass transfer across the EMB surface, this second-order nonlinear ordinary differential equation (ODE) characterizes the dynamic behavior of a viscoelastic-shelled microbubble under ultrasound. Using the EEM, the modified RPE given by Equation 40 reads (Sarkar et al., 2005):Here, the overdot denotes differentiation with respect to time, R is the instantaneous EMB radius, and the internal bubble pressure is expressed using the polytropic relation as P b ( t ) = P g R 0 / R 3 κ, where P g is the equilibrium gas pressure inside the bubble evaluated from the Y-L law as P g = P 0 + 2 σ 0 / R 0. Furthermore, ν l is the liquid kinematic viscosity, and κ denotes the polytropic exponent. Using the initial conditions R ( t = t 0 ) = R 0 and R ̇ ( t = t 0 ) = 0, and assuming isothermal gas behavior inside the bubble, i.e., κ = 1, the modified RPE is solved through the Python ODE solver odeint. The isothermal assumption is consistent with the isothermal nature of the MCMP IB-LBM developed in this study. Water and C 4 F 10 at 25 °C and the reduced temperature of T r = 0.6623, which form the equilibrium density ratio of ρ l / ρ g = 100.3573, are taken as the liquid and gas phases, respectively. The physical parameters used in solving the modified RPE are: P 0 ph = 1 a t m, ρ l ph = 997.05 k g / m 3, ν l ph = 8.917 × 1 0 − 7 m 2 / s, R 0 ph = 1.6 μ m, and F ph = 3 M H z. Also, the interfacial rheological parameters are considered as: σ 0 ph = 0.019 N / m, E s ph = 0.55 N / m, κ s ph = 0.012 m s P = 0.012 × 1 0 − 6 N ⋅ s / m, and α s = 1.5, following Paul et al. (Paul et al., 2010). They determined the rheological values through matching linearized EEM dynamics with attenuation measurements. These rheological values are representative of lipid-coated microbubbles commercially available as Sonazoid. To ensure equilibrium similarity between the lattice and physical systems (Krüger et al., 2017), the equilibrium EMB radius of R 0 = 24.7241 in the lattice unit is chosen from the law of similarity for the Laplace number La, i.e., L a = σ 0 ρ l R 0 / ( μ l ) 2 = σ 0 ph ρ l ph R 0 ph / ( μ l ph ) 2, where symbols with the superscript ‘ph’ denote parameters in physical units, and those without it represent lattice values. The pressure variation amplitude Δ P a in the lattice unit is related to its physical counterpart Δ P a ph according to Equation 41 (Viggen, 2014):and the frequency F in the lattice unit is associated with the physical frequency F ph via the law of similarity for the Womersley number Wo, i.e., W o = R 0 2 π F / ν l = R 0 ph 2 π F ph / ν l ph (Krüger et al., 2017).

R R ̈ + 3 2 R ̇ 2 = − 4 ν l R R ̇ − 2 σ ρ l R − 4 k s R ̇ ρ l R 2 + P b t − P a t ρ l . (40)

Δ P a = Δ P a ph ρ l ρ l ph R 0 ph ν l R 0 ν l ph 2 , (41)

When the EMB is driven for 10 acoustic cycles, Figure 6 illustrates the influence of the stencil employed in the discretized Dirac delta function on the accuracy of the IB-LBM predictions for the normalized EMB radius ( R / R 0 ). The results obtained using the stencils ϕ 2, ϕ 3, and ϕ 4, are compared against the solution of the modified RPE. The stencils ϕ 3 and ϕ 4 are defined by Equations 42,43 as (Peskin, 2002):andrespectively, where indices ‘3’ and ‘4’ specify that three and four Eulerian nodes are used along each coordinate axis r and z for interpolation and spreading by the stencils ϕ 3 and ϕ 4, respectively. Looking at Figure 6, all three stencils overall are capable of reproducing the oscillatory response of the EMB reasonably well. However, noticeable differences in prediction accuracy are observed. The lower-order stencil ϕ 2 produces the largest deviation from the modified RPE solution, particularly near the oscillation peaks where the bubble expansion is overpredicted. The stencil ϕ 3 improves the agreement with the theoretical solution, whereas the stencil ϕ 4 yields the closest match throughout the simulation. As shown more clearly in the enlarged view in Figure 6B, increasing the stencil order reduces the discrepancy between the IB-LBM and modified RPE results, particularly in the oscillation amplitude, while also slightly improving the temporal alignment of the oscillations. The improved performance of the higher-order stencils can be attributed to their smoother interpolation and spreading characteristics between the Eulerian and Lagrangian grids. By distributing the IB forcing over a wider support region, the higher-order stencils reduce numerical errors associated with the fluid-structure coupling and provide a more accurate representation of the bubble interface dynamics. Consequently, among the tested stencils, ϕ 4 demonstrates the highest accuracy in capturing the spherical oscillations of the microbubble.

ϕ 3 x = 1 + 1 − 3 x 2 3 , 0 ≤ | x | ≤ h 2 , 5 − 3 | x | − − 2 + 6 | x | − 3 x 2 6 , h 2 < | x | ≤ 3 h 2 , 0 , | x | > 3 h 2 , (42)

ϕ 4 x = 3 − 2 | x | + 1 + 4 | x | − 4 x 2 8 , 0 ≤ | x | ≤ h , 5 − 2 | x | − − 7 + 12 | x | − 4 x 2 8 , h < | x | ≤ 2 h , 0 , | x | > 2 h , (43)

Figure 7 presents the effect of the advection scheme used for updating the Lagrangian marker positions on the accuracy of the IB-LBM predictions for the normalized EMB radius ( R / R 0 ). The results obtained using the first-order Euler, second-order Adams-Bashforth, and fourth-order Runge-Kutta schemes are compared against the solution of the modified RPE. The Adams-Bashforth scheme is expressed according to Equation 44 (Doddi and Bagchi, 2008):and the Runge-Kutta method advects the Lagrangian markers based on Equation 45 (Chapra and Canale, 2015):where the stage abscissas are defined by Equation 46 as:The unknown velocity vectors are evaluated as: with the asterisked position vectors for the Lagrangian markers defined by Equations 48a, b: Since the Eulerian grid velocity field u b at times t + Δ t / 2 and t + Δ t is not yet available, u b at time t is used in Equations 47a,b. This treatment resembles the lagged-variable strategy frequently used in deferred-correction-type CFD formulations (Versteeg and Malalasekera, 2007). As seen in Figure 7, all three advection schemes successfully reproduce the oscillatory behavior of the bubble. Nevertheless, differences in numerical accuracy are evident. Among the tested schemes, the Euler method exhibits the largest deviation from the modified RPE solution, particularly near the oscillation peaks where the bubble expansion amplitude is slightly overpredicted. The Adams-Bashforth method improves the agreement with the theoretical solution, while the Runge-Kutta scheme provides the closest match throughout the simulation. As illustrated more clearly in the enlarged view in Figure 7B, increasing the order of the advection scheme reduces the discrepancy between the IB-LBM and modified RPE results in both oscillation amplitude and temporal alignment. The higher-order schemes more accurately track the motion of the IB markers and therefore reduce numerical errors associated with interface advection. Consequently, the Runge-Kutta method yields the most accurate prediction of the spherical microbubble dynamics among the schemes considered in this study.

r j t + Δ t = r j t + 3 2 r ̇ j t − 1 2 r ̇ j t − Δ t Δ t , (44)

r j t + Δ t = r j t + Δ t 6 ζ 1 + 4 ζ 2 + ζ 3 , (45)

ζ 1 = r ̇ j t ; ζ 2 = r ̇ j t + Δ t 2 ; ζ 3 = r ̇ j t + Δ t . (46)

r ̇ j t + Δ t 2 = ∑ x h 2 u b x , t D r * j t + Δ t 2 − x , (47a)

r ̇ j t + Δ t = ∑ x h 2 u b x , t D r * j t + Δ t − x , (47b)

r * j t + Δ t 2 = r j t + r ̇ j t Δ t 2 , (48a)

r * j t + Δ t = r j t + r ̇ j t Δ t . (48b)

Figure 8 demonstrates the capability of the IB-LBM to accurately capture nonlinear acoustic oscillations of an EMB subjected to sinusoidal forcing, where substantial departures from purely sinusoidal behavior occur. Unlike the previously presented linear-regime cases, the present response exhibits strongly nonlinear oscillations characterized by asymmetric expansion and compression phases, variations in oscillation amplitude, and the appearance of higher-harmonic features in the time history of the normalized EMB radius ( R / R 0 ). The IB-LBM predictions obtained using the different discretized Dirac delta stencils remain in close agreement with the modified RPE solution throughout the simulation. Although differences among the stencils ϕ 2, ϕ 3, and ϕ 4 are still observable, particularly near the oscillation peaks, all stencils successfully reproduce the main nonlinear features of the bubble dynamics. The higher-order stencil ϕ 4 again yields the closest agreement with the modified RPE solution. Overall, the figure confirms that the proposed IB-LB framework is capable of accurately modeling nonlinear acoustic bubble dynamics in addition to the linear oscillatory regimes examined previously.

Figure 9 further validates the ability of the IB-LBM to reproduce nonlinear acoustically-driven EMB dynamics through comparison of different Lagrangian marker advection schemes. In contrast to the earlier linear-regime cases, the present oscillations are distinctly nonlinear and exhibit significant deviations from sinusoidal behavior, including unequal oscillation amplitudes, waveform distortion, and enhanced nonlinear expansion-compression asymmetry. The normalized EMB radius ( R / R 0 ) histories predicted using the Euler, Adams-Bashforth, and Runge-Kutta advection schemes are compared against the modified RPE solution. The results indicate that all three advection schemes are capable of capturing the nonlinear oscillatory behavior with good agreement relative to the modified RPE solution. As in the previous cases, the higher-order Runge-Kutta method provides the closest agreement, especially near the oscillation extrema where nonlinear effects become more pronounced. These results demonstrate that the proposed IB-LB framework remains accurate and stable even in nonlinear acoustic regimes where departures from harmonic oscillations become substantial.

When the EMB is driven for 10 acoustic cycles and the stencil ϕ 4 and the Runge-Kutta advection scheme are used, Figure 10 presents the IB-LBM contours of mixture density and barycentric velocity magnitude at four representative instants. The solid black line overlaying the barycentric velocity magnitude contours represents the EMB surface constructed by the Lagrangian markers. The mixture density contours clearly distinguish the gas core (blue region) from the surrounding liquid (red region) while capturing the finite-thickness diffuse interface characteristic of the MCMP LB formulation. As the microbubble undergoes successive expansion and compression, the diffuse interface maintains a continuous density transition between the gas and liquid phases, demonstrating the robustness of the proposed IB-LBM in resolving large interface motions. The corresponding barycentric velocity magnitude contours illustrate the fluid motion induced by the oscillating microbubble. During expansion phases, the outward radial motion of the shell displaces the surrounding liquid, producing a relatively broad region of elevated velocity adjacent to the interface. In contrast, during collapse, the contracting shell accelerates the surrounding liquid toward the bubble, generating a localized region of higher velocity concentrated near the interface. The largest velocity magnitudes consistently occur immediately adjacent to the shell, where the interfacial motion directly transfers momentum to the surrounding liquid, while the velocity rapidly decreases with increasing distance from the bubble, for the disturbance propagates outward and dissipates over a larger fluid volume. Comparison of the four snapshots further illustrates the strong coupling between the interfacial motion and the surrounding flow field. Larger bubble sizes are associated with more spatially extended velocity fields owing to the larger displaced liquid volume, whereas smaller bubble sizes produce a more compact but stronger near-interface flow. These results demonstrate that the proposed MCMP IB-LBM accurately resolves both the evolving gas-liquid interface and the accompanying transient viscous flow structures throughout the acoustic exposure, thereby providing a physically resolved description of ultrasound-driven microbubble dynamics.

Moreover, when the EMB is driven for 10 acoustic cycles and the stencil ϕ 4 and the Runge-Kutta advection scheme are employed, Figure 11 presents the IB-LBM time histories of the elastic ( g elas ), viscous ( g vis ), and total ( g elas + g vis ) shell forces together with the time history of the applied acoustic pressure ( P a ) during the last two oscillation cycles. The elastic part of the shell force is defined as g elas = − 2 H j = 0 E A j = 0 n j = 0, where j = 0 denotes the Lagrangian marker at the north pole of the EMB and the brackets ‘| |’ signify the magnitude of a vector. The viscous part is expressed by Equation 49 as follows:

g vis = − 2 μ s B j = 0 − 2 H j = 0 I j = 0 : ∇ r ̇ j = 0 − 2 H j = 0 κ s + μ s ∇ ⋅ r ̇ j = 0 A j = 0 n j = 0 + g t j = 0 . (49)

Based on Figure 11, it is seen that the shell forces respond in a distinctly nonlinear manner because they depend on the instantaneous shell deformation and deformation rate. The elastic component increases whenever the shell undergoes significant surface area variation, as the shell elasticity generates restoring stresses that oppose interface deformation. Accordingly, pronounced elastic force peaks occur near the extrema of the bubble oscillation, where the shell experiences its largest deformation. In contrast, the viscous component is governed by the rate of shell deformation and therefore becomes most significant during periods of rapid expansion or contraction. The total shell force reflects the combined action of these elastic restoring and viscous dissipative mechanisms. The largest total-force peaks occur during rapid, large-amplitude oscillations when both the shell deformation and deformation rate are simultaneously high. Conversely, the total shell force approaches zero as the shell passes through states of relatively small deformation and low deformation rate. These results demonstrate that the proposed IB-LBM successfully resolves the time-dependent viscoelastic response of the encapsulating shell under ultrasound excitation.

3.3 IB-LBM analysis of shell property effects on ultrasound-driven EMBs

Having previously established the accuracy of the proposed IB-LB framework through comparisons with the modified RPE, the effects of the shell mechanical properties on the nonlinear acoustic response of EMBs are investigated here. In particular, the influence of shell elasticity E s and dilatational viscosity κ s on the EMB oscillatory behavior under ultrasound excitation is explored. The obtained results indicate that the present IB-LBM is capable of reproducing rich nonlinear physics associated with viscoelastic shell dynamics and their interaction with the surrounding fluid.

Figure 12 illustrates the effect of shell elasticity on the normalized EMB radius ( R / R 0 ) when the shell dilatational viscosity is ignored. The results show that introducing shell elasticity substantially modifies the acoustic response of the EMB. Compared with the uncoated microbubble, the encapsulated cases exhibit changes in oscillation amplitude, waveform shape, and long-term oscillatory behavior. In particular, increasing the shell elasticity initially strengthens the elastic restoring stresses opposing surface area expansion, thereby suppressing the early radial growth of the EMB. This observation is consistent with the one reported by Mobadersany and Sarkar (Mobadersany and Sarkar, 2018). It is seen that the stiffer shell also enhances the mechanical stability of the interface through increasing the duration of sustained oscillations. This is because a stiffer shell stores a larger portion of the acoustic energy elastically, consequently reducing irreversible bubble collapse tendencies and thereby stabilizing the EMB against dissolution (Church, 1995). Furthermore, due to nonlinear elastic energy storage and the strain-softening behavior of the EEM employed in the present work, larger shell elasticities eventually promote larger subsequent expansions after the initial oscillation stage. These results confirm that the proposed IB-LBM is capable of resolving the nonlinear coupling between interfacial elasticity and multiphase fluid motion. It is also worth noting that the consistency of the trends shown in Figure 12 with those observed by Mobadersany and Sarkar (Mobadersany and Sarkar, 2018) and by Church (Church, 1995) provides independent support for the physical behavior predicted by the proposed IB-LBM.

Figure 13 presents the effect of shell dilatational viscosity κ s on the EMB acoustic response when the shell elasticity is E s ph = 0.55 N / m. Introducing shell dilatational viscosity noticeably alters the EMB dynamics by damping the oscillatory motion and reducing the oscillation amplitude. As the shell dilatational viscosity increases, the oscillations become progressively attenuated, and the EMB response becomes smoother and less nonlinear. The results further show that shell dilatational viscosity strongly affects the persistence and stability of the oscillations. The purely elastic shell produces large-amplitude oscillations with substantial nonlinear behavior, whereas increasing κ s suppresses extreme expansions and compressions of the bubble surface. Physically, this behavior originates from viscous dissipation within the shell, which converts part of the acoustic energy into dissipative losses and therefore weakens the EMB oscillatory response. It should be noted that the simulation case when κ s ph = 0.012 m s P was previously validated against the corresponding modified RPE solution, and the remaining cases differ only in the shell dilatational viscosity value while employing the same governing equations and numerical implementation.

Overall, Figures 12 and 13 demonstrate that the proposed IB-LB framework can successfully capture the complex influence of shell viscoelastic properties on nonlinear ultrasound-driven EMB dynamics. The solver accurately reproduces how shell elasticity modifies interface stiffness and nonlinear oscillatory behavior, while shell dilatational viscosity introduces damping and energy dissipation effects. These results highlight the capability of the present IB-LB methodology to provide detailed physical insight into the mechanics of EMBs subjected to acoustic forcing.

3.4 Nonspherical collapse of an EMB near a membrane induced by ultrasound

To further evaluate the capability of the MCMP IB-LBM developed in this study for simulating more complex EMB dynamics, the nonspherical collapse of an ultrasound-driven EMB in the vicinity of a membrane is modeled. This problem is of particular interest because ultrasound-induced EMB-membrane interactions play a central role in sonoporation, a technique widely employed to enhance intracellular drug and gene delivery (Tachibana et al., 1999; Pitt et al., 2004; Kooiman et al., 2011; Escoffre et al., 2013). In the present simulation, C 4 F 10 is used as the gas encapsulated within the EMB, while water serves as the surrounding liquid medium. The left boundary of the computational domain illustrated in Figure 2 is kept as the axis of symmetry, the bottom boundary is assigned a no-slip condition, the top boundary employs the Zou-He pressure BC, and the right boundary is modeled using a free-slip condition. Also, the computational domain dimensions are set to L r = 200 and L z = 400 in lattice units. The microbubble is initially at static equilibrium and at t > t 0, an oscillatory acoustic pressure, P a ( t ) = P 0 − Δ P a ⁡ sin ( 2 π F t ), is imposed at the top boundary to initiate the insonation process. To validate the IB-LBM predictions, the computed EMB dynamics are compared with the numerical results reported by Mobadersany and Sarkar (Mobadersany and Sarkar, 2018), who used the boundary integral method to solve the problem. To ensure consistency with their study, the following physical parameters are employed: the equilibrium bubble radius is specified as R 0 ph = 1.6 μ m, the excitation frequency is chosen as F ph = 2 M H z, and the reference surface tension and shell elasticity are prescribed as σ 0 ph = 0.019 N / m and E s ph = 0.55 N / m, respectively. Furthermore, the shell dilatational viscosity is neglected, the dimensionless shell parameter is set to α s = 1.5, the dimensionless acoustic pressure amplitude is taken as ϵ = 5, and the offset parameter, representing the ratio of the initial bubble-center-to-membrane distance to the equilibrium bubble radius, is specified as λ = 4. It is to be noted that the stencil ϕ 4 and the Runge-Kutta scheme are employed for respectively evaluating the discretized Dirac delta function and advancing the Lagrangian markers on the EMB surface within the IB-LBM.

Figure 14 compares the EMB shapes predicted by the present IB-LBM with those computed through the boundary integral method by Mobadersany and Sarkar (Mobadersany and Sarkar, 2018) at different times during the insonation process. The close agreement between the two sets of results demonstrates that the proposed MCMP IB-LBM accurately reproduces the nonspherical collapse dynamics of an EMB adjacent to a membrane. Figure 15 presents the IB-LBM predictions of the mixture density and barycentric velocity magnitude at different times during the insonation process of an EMB near a membrane. The left half of each panel illustrates the mixture density contours, which delineate the gas-liquid interface, while the right half shows the corresponding barycentric velocity magnitude of the fluid mixture. Together, these fields provide insight into the evolution of the surrounding flow and the mechanisms governing jet formation and EMB collapse.

Figure 15A corresponds to the early expansion stage shortly after the acoustic excitation is applied. The rarefaction phase of the ultrasound, during which the acoustic pressure drops below the hydrostatic liquid pressure, causes the gas inside the EMB to expand, producing an outward displacement of the surrounding liquid. The velocity field remains nearly radially symmetric, with relatively high velocity magnitudes seen around the bubble interface inside the gas region. The nearly concentric velocity contours indicate that inertial effects are still weak and that the oscillation remains predominantly spherical. The mixture density contours show a smooth, well-defined gas-liquid interface with only minor deformation. Figure 15B depicts the instant at which the EMB reaches its maximum volume while maintaining an almost spherical shape. During the expansion phase, the surrounding liquid acquires outward momentum. Although Figure 15B corresponds to the maximum EMB volume, the acoustic pressure imposed at the upper boundary has already entered its compressive phase before t = 1341. Because the pressure disturbance must propagate from the upper boundary toward the EMB, the bubble response exhibits a finite delay relative to the imposed pressure signal. Thus, at t = 1341, the maximum size of the EMB is achieved, while the surrounding liquid has already begun to respond to the developing compressive pressure field. The compressive pressure wave approaches from the upper side, while the bottom no-slip membrane restricts the lower-side liquid motion and introduces viscous resistance. The free-slip side permits stronger tangential motion and less resistance than the membrane side, causing a high-velocity patch observed near the east pole of the EMB. Consequently, the flow is not driven or constrained symmetrically around the EMB. Figure 15C illustrates the early stage of jet formation. As the collapse intensifies, the pressure gradient above the EMB becomes substantially larger than that below it, because of the nearby membrane and the asymmetric flow field. The liquid above the bubble is therefore accelerated downward, giving rise to a narrow, high-speed jet directed toward the membrane. This process is reflected by the localized region of elevated barycentric velocity immediately above the EMB. Simultaneously, the mixture density contours reveal a pronounced indentation of the upper bubble surface, indicating that the initially spherical interface has transitioned to a strongly nonspherical configuration. Lastly, Figure 15D shows the final stage of collapse, during which the liquid jet rapidly penetrates the collapsing EMB and impinges on the membrane. At this stage, the velocity magnitude reaches its maximum value owing to the strong conversion of the stored potential energy accumulated during the expansion phase into liquid kinetic energy. The concentrated high-velocity region surrounding the jet tip reflects the intense momentum transfer associated with jet impact. The mixture density contours demonstrate that the bubble has fragmented into a highly distorted shape, while the velocity field indicates heightened hydrodynamic loading on the membrane, i.e., the liquid velocity adjacent to the membrane at t = 1842 is approximately 50 times larger than when t = 100. The high-speed microjet generated during the final stage of collapse is recognized as the principal mechanism responsible for membrane permeabilization during sonoporation.

4 Discussion

Proposed RP-type formulations in the literature have historically served as the principal theoretical framework for investigating coated bubble acoustics because of their computational efficiency and analytical simplicity. Nevertheless, these reduced-order formulations inherently assume spherical symmetry and therefore cannot directly resolve surrounding flow structures or complex fluid-structure interactions, e.g., nonspherical deformations of EMBs near membranes induced by ultrasound. On the other hand, prior investigations of nonspherical EMB dynamics using boundary element or boundary integral formulations neglect liquid viscosity through the assumption of potential flow in the surrounding medium. This assumption may become increasingly restrictive in biomedical ultrasound applications where viscous dissipation plays a significant role, particularly during large-amplitude oscillations or near-wall bubble dynamics.

Motivated by these shortcomings, the present study introduces a computational framework based on IB-LB coupling for simulating ultrasound-excited EMBs. The proposed methodology combines an axisymmetric MCMP LB flow solver with an IB representation of the viscoelastic shell, thereby enabling fully resolved two-way coupling between the surrounding fluid phases and the EMB shell within an acoustically-driven environment. Through direct solution of the continuity and Navier-Stokes equations in LB form, the developed IB-LBM naturally incorporates surrounding liquid viscosity, and using the Boussinesq-Scriven constitutive law along with an EEM, the shell mechanics are integrated into the proposed framework.

The Y-L validation study confirms that the developed MCMP IB-LB formulation can accurately reproduce equilibrium pressure differences across curved interfaces over a wide range of bubble sizes. The predicted linear relationship between the interfacial pressure jump and inverse bubble radius in equilibrium demonstrates consistency with the Y-L relation, while the realistically high liquid-to-gas density ratios achieved by the solver verify its capability to model physically relevant multiphase conditions. These results are important because accurate representation of interfacial tension and MCMP equilibrium forms the foundation for reliable simulation of acoustically-driven EMB dynamics.

The simulations of ultrasound-driven EMB oscillations further demonstrate that the proposed IB-LBM accurately reproduces the predictions of a modified RPE in spherical oscillation regimes. The close agreement between the IB-LBM predictions and the modified RPE solutions therefore illustrates that the present fully resolved CFD framework remains consistent with established EMB acoustics theory while providing substantially greater modeling generality. The stencil-sensitivity and advection-scheme studies additionally demonstrate that the numerical accuracy of the proposed framework depends on the treatment of IB’s interpolation-spreading and interface advection. Higher-order interpolation-spreading stencils and higher-order advection schemes consistently improve agreement with the modified RPE solution by reducing numerical errors associated with velocity interpolation, force spreading, and Lagrangian marker motion. In particular, the fourth-order stencil and Runge-Kutta advection scheme yield the closest agreement with the reference solution. These observations are physically consistent with the role of IB’s interpolation-spreading and surface advection in accurately exchanging momentum and force between Eulerian and Lagrangian grids.

Importantly, it is seen that the proposed IB-LBM remains accurate in nonlinear acoustic regimes where the EMB response departs substantially from purely sinusoidal oscillations. The predicted oscillation histories exhibit strong nonlinear behavior, including asymmetric expansion-compression cycles and complex waveform modulation. Moreover, the shell property studies further demonstrate that the proposed framework can reproduce complex viscoelastic effects associated with EMB encapsulation. Increasing shell elasticity significantly alters the nonlinear oscillatory response by modifying the interface stiffness and changing the energy exchange between the acoustic field and the bubble interface. Likewise, increasing shell dilatational viscosity introduces substantial damping and reduces oscillation amplitudes through viscous dissipation within the shell. These findings are consistent with those reported in prior works that emphasize the importance of shell rheology in determining coated bubble dynamics.

The presented simulation of the nonspherical collapse of an ultrasound-driven EMB near a membrane demonstrates the broader applicability of the proposed IB-LB framework. In contrast to reduced-order RP-type models, which are restricted to spherical oscillations, and potential-flow-based formulations, which neglect viscous effects in the surrounding liquid, the present methodology directly resolves the coupled evolution of the gas-liquid interface, the surrounding viscous flow field, and the viscoelastic shell during strongly nonspherical collapse. The computed mixture density and barycentric velocity fields provide detailed insight into the physical mechanisms governing jet initiation, jet development, and membrane impact, while the close agreement with the established numerical results validates the capability of the proposed solver to accurately reproduce complex nonspherical EMB dynamics. These results highlight the potential of the IB-LBM as a predictive CFD tool for investigating ultrasound-mediated bioeffects, including sonoporation and targeted drug delivery, where fluid viscosity, wall interactions, and large interfacial deformations play essential roles.

When modeling the radial fluctuations of an EMB, the proposed MCMP IB-LBM is inherently more computationally demanding than solving a modified RPE, because it resolves the complete multiphase flow field, fluid-structure interaction, and shell mechanics throughout the computational domain. Consequently, simulations require substantially more computational resources than reduced-order models. For example, when the stencil ϕ 4 and the Runge-Kutta advection scheme are employed and using a serial implementation on a desktop computer equipped with an Intel® Core™ i7-3610QM processor (2.30 GHz) and 8 GB of RAM, a representative IB-LBM simulation of the radial oscillations of an EMB on a computational domain of 401 × 801 Eulerian nodes requires 6200 time steps and takes approximately 177 min, whereas the corresponding modified RPE solution requires a few seconds. Although the computational expense is considerably higher, the IB-LBM enables direct simulation of nonspherical deformation, viscous flow structures, and EMB-boundary interactions, that are inaccessible to reduced-order models. It is also worth noting that the reported IB-LBM runtime of a representative serial implementation could be substantially reduced through parallelization, which is well suited to LB-based frameworks.

Several future research directions naturally emerge from this study. The most immediate extension involves the simulation of EMB shape mode oscillations that may evolve into jet formation, fragmentation, and shell rupture under sufficiently strong acoustic excitation. Since the present framework inherently includes surrounding medium’s viscous effects and moving deformable boundaries, it may also provide new physical insights into physiologically more realistic near-membrane EMB dynamics and sonoporation mechanisms by resolving the bidirectional fluid-structure interaction between an oscillating EMB and a deformable membrane, together with the surrounding viscous flow. Furthermore, extension of the framework to fully 3-D simulations would allow investigation of asymmetric instability modes and complex bubble interactions that cannot be captured in axisymmetric formulations. Another promising direction involves coupling the present framework with models for heat transfer, mass transport, or biological tissue mechanics to enable multiscale investigation of ultrasound-mediated drug delivery, mechanotransduction, and bioeffects.

The present work raises several broader theoretical questions that could be explored in future studies. For example, the interaction between shell viscoelasticity and surrounding flow viscosity may influence the onset of nonspherical instabilities and cavitation thresholds in ways not captured by classical RP-type formulations. Likewise, the competition between shell elasticity, shell viscosity, and near-wall viscous stresses may govern transition pathways among stable oscillations, asymmetric collapse, and membrane poration. Fully resolved CFD methodologies such as the present IB-LB framework provide a promising platform for testing such hypotheses and for advancing the current understanding of ultrasound-driven EMB physics in biologically relevant environments.

5 Limitations

The present IB-LB framework was developed to provide a robust and physically resolved computational methodology for simulating ultrasound-driven EMBs in biologically relevant environments. Nevertheless, limitations of the current implementation should be acknowledged.

First, the present framework employs a single-relaxation-time (SRT) LB formulation for each fluid component. As a result, the range of liquid-to-gas dynamic viscosity ratios that can be achieved is limited. In the current implementation, the proposed solver is capable of handling liquid-to-gas dynamic viscosity ratios below approximately 250. However, this limitation does not significantly restrict the applicability of the framework to many biomedical ultrasound problems involving EMBs, since the effective viscosity ratios encountered in numerous biologically relevant applications fall within this range. To be specific, the viscosity ratio of the MCMP fluid system considered in this study is μ water / μ C 4 F 10 = 68.974 at 25 °C. As the viscosity ratio is increased beyond the above-mentioned limit, the SRT formulation becomes increasingly susceptible to numerical instability, manifested by amplified spurious currents near the gas-liquid interface, deterioration of interface sharpness, and, in severe cases, divergence of the numerical solution. These instabilities arise because different kinetic moments relax at a single rate, reducing the ability of the method to independently damp nonphysical modes. If simulations involving substantially larger viscosity contrasts become necessary, the present IB-LBM can be extended through incorporation of a multi-relaxation-time (MRT) LB formulation, which is known to provide improved numerical stability for high-viscosity-ratio multiphase flows. The MRT formulation requires additional moment-space transformations and therefore is expected to modestly increase the computational cost.

Second, the present formulation is not intended for strongly compressible flow regimes, such as supersonic flows involving strong shock wave propagation. This limitation arises from the weakly compressible nature of the LBM employed in the current framework. Nevertheless, this restriction is not expected to affect the intended biomedical ultrasound applications of the present work. In many EMB-related biomedical applications, the surrounding liquid can reasonably be treated as incompressible, while the gas inside the microbubble typically remains within weakly compressible acoustic regimes under diagnostic and many therapeutic ultrasound conditions. Consequently, the present IB-LB methodology remains well-suited for modeling the dominant fluid-structure interaction physics associated with ultrasound-driven EMB dynamics in biologically relevant environments.

On the whole, these limitations primarily reflect the scope of the current numerical implementation rather than fundamental restrictions of the IB-LB coupling strategy itself. The modular nature of the proposed framework provides a clear pathway for future extensions involving higher-viscosity-ratio multiphase flows, improved compressibility treatment, and more advanced LB-based formulations.

6 Conclusion

In the present study, a computational framework based on IB-LB coupling was developed for simulating ultrasound-excited EMBs. The proposed methodology combined an axisymmetric MCMP LB flow solver with an IB representation of the EMB viscoelastic shell, thereby enabling fully resolved two-way coupling between the surrounding fluid phases and the shell within an acoustically-driven environment.

The developed IB-LBM was validated through the Y-L test for a 3-D bubble and through comparison with solutions of a modified RPE for ultrasound-driven EMB oscillations. The results demonstrated that the proposed framework accurately reproduces interfacial pressure jumps, physically relevant liquid-to-gas density ratios, and spherical EMB oscillations under both linear and nonlinear acoustic regimes. In addition, higher-order interpolation-spreading stencils and higher-order advection schemes were found to improve the numerical accuracy of the IB-LBM predictions.

The proposed methodology also successfully captured the influence of shell viscoelastic properties on nonlinear EMB dynamics. It was shown that increasing shell elasticity alters the interface stiffness and modifies the nonlinear oscillatory response, while increasing shell dilatational viscosity was shown to introduce stronger damping and dissipative effects. Moreover, the successful simulation of jet formation and nonspherical collapse of an ultrasound-driven EMB near a membrane illustrated that the proposed IB-LBM can accurately resolve complex viscous flow structures and fluid-structure interactions beyond the capabilities of classical RP-type and potential-flow-based formulations, thereby extending its applicability to biologically relevant ultrasound-mediated phenomena. These simulations demonstrated the capability of the present IB-LB framework to reproduce complex fluid-structure interaction physics associated with ultrasound-driven EMBs.

Altogether, the present work demonstrated that the proposed IB-LBM provides a robust and physically resolved CFD framework for modeling EMB dynamics in biologically relevant acoustic environments. The methodology established in this study provides a foundation for future investigations of further nonspherical EMB dynamics (e.g., EMB shape mode oscillations), fully bidirectional bubble-boundary interactions, cavitation-based tissue ablation, and other biomedical ultrasound phenomena involving complex multiphase fluid-structure coupling.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: John Sharer Allen, University of Hawaii at Manoa, United States

Reviewed by: Phuong Nguyen, UPR9080 Laboratoire de Biochimie Théorique (LBT), France

Zhenxiang Ji, Beijing Institute of Technology, China

Qianxi Wang, University of Birmingham, United Kingdom