2.1. LBM model for flow field

The LBM model has three ingredients. The first ingredient is a discrete phase space defined by a regular lattice together with a set of symmetric discrete velocities { e _α | α = 0,1, …, N } and density distribution functions { f _α | α = 0,1,…, N }. The second ingredient includes a collision matrix S and equilibrium distribution functions . The third ingredient is the evolution equation

where S is the collision matrix. In the velocity space, the full matrix S results in the difficulty to implement the collision step. The full collision matrix in the moment space can be transformed into a diagonal matrix, which supports the convenience in performing the collision process. ^{[ 16 – 18 ]} Equation ( 1 ) is transformed to

where m ( x ,t) and m ^eq ( x ,t) are moments and their corresponding equilibriums, and S̃ is a diagonal matrix whose components represent the relaxation times of the moments. The transformation between the velocity space and the moment space is achieved by using matrix M , which serves as a transformation matrix and maps the distribution functions f ( x , t ) to their moments as

The transformation matrix M is constructed via the Gram–Schmidt orthogonalization procedure from some polynomials of the discrete velocities component.

In this work, the D3Q19 model is adopted as shown in Fig. 1 , and the particle velocity directions are

where c = δ _x / δ _t is the lattice velocity.

The equilibrium distribution moments m ^eq are related to the conserved moments as

The nineteen equilibrium distribution moments of this model are

where j = ( j _x ,j _y ,j _z ) is the momentum, are the components of the energy flux, , and are the components of the symmetric traceless viscous stress tensor, and are parts of a third-rank tensor. The collision matrix S̃ in the moment space is

The parameters s ₀ , s ₃ , s ₅ , and s ₇ are the relaxation times corresponding to the collision invariants. The parameters s ₉ , s ₁₁ , s ₁₃ , and s ₁₅ are related to the kinematic viscosity ν by , where .

The fluid macroscopic fluctuation density δρ and velocity u are computed from the moments of the distribution functions as

2.2. Front-tracking method

When the evolution equation ( 2 ) is solved numerically on an Eulerian grid, the fluid properties, i.e., density and viscosity, on this grid are required. In terms of the two-phase flow, each phase is assumed as an incompressible flow, and the fluid properties are taken as constant in each phase. Thus the properties of the fluid are physically discontinuous across the interface, and this discontinuity at grid points adjacent to the interface results in many difficulties in simulating numerically the two-phase flow. In order to simulate accurately the interface movement, Tryggvason et al . ^{[ 1 ]} adopted an Eulerian grid as a background grid to solve the governing equations of flow, and applied an additional Lagrangian grid to track explicitly the interface motion. The discontinuities across the interface are distributed from the interface to the background grid, and continuous distributions of the fluid properties can be reconstructed on the Eulerian grid. In this study, the entire flow domain is solved by only one MRT LBM evolution equation ( 2 ) and the different fluids are treated as one fluid with varying properties across the interface. The mixed Eulerian and Lagrangian grids are exhibited in Fig. 2 .

The reconstruction of the fluid properties b ( x , t ) at time t in the whole background grid is achieved through constructing an indicator function I ( x , t ). Let this indicator function be zero in one phase and one in another phase, then we can define

where b ( x , t ) is either fluid density or viscosity, and the subscripts 0 and 1 refer to fluid 1 and fluid 2, respectively. According to the interface location, the indicator function I ( x , t ) is obtained by solving the following Poisson equation:

The Poisson equation is solved by a standard second-order centered difference method, and the linear equation is solved with a traditional successive-over-relaxation method. Here G ( x , t ) is related to the interface location and normal

where n _f is the unit normal vector of the interface, Γ represents the interface between the two phases, x _f is the coordinates of the interface points, and D ( x − x _f ) is a distribution function for the interface. In this study, we adopt the traditional Peskin distribution function to approximate the D ( x − x _f ) function

The partial derivatives in Eq. ( 11 ) are computed as

where w _α is the weight coefficient

This discretization is based on the lattice-based stencil, instead of the standard stencil based on the coordinate directions. A smooth indicator function can be constructed by solving the Poisson equation ( 11 ). This indicator function is zero in one phase, varies continuously from zero to one in a certain thickness region, and one in another phase.

The governing equation ( 2 ) is solved on the stationary background grid. And then the velocity of the interface makers is interpolated from the flow field on the background grid. Subsequently, the interface makers move in a Lagrangian manner from a position at the n time step to a position at the n + 1 time step

where the velocity interpolation is carried out using the same distribution function as Eq. ( 13 ).

As the interface makers move, the interface grid resolution may change, which has a strong effect on the information exchange between Lagrangian grid and Eulerian grid. Therefore, an adaptive Lagrangian grid is adopted to maintain the resolution during the simulation. In this study, restructuring the Lagrangian grid for the interface consists of two operations, edge split and edge deletion. For a short edge, the edge is removed, and for a long edge, an additional maker is placed in the middle of the edge. Generally, we take the maximum edge length as the Eulerian grid size, and the minimum edge length as 0.3 times the Eulerian grid size.

2.3. Surface tension

Calculating accurately the surface tension is a great challenge in simulating the two-phase flow. In this study, the surface tension is computed from the tensile forces on the three edges of all interface elements m (see Fig. 3 )

where t _k is the tangent vector of the edge shared by element m and neighboring element k , and n _k is its unit normal vector. The summation in Eq. ( 19 ) is carried out over all three edges of the element m . The three tangent vectors t _k of the element m can be computed directly from the known positions of the three corner points of the element. Once the tangent vectors are obtained, the unit normal vectors can be found from the cross-product of two different tangent vectors.

3.1. Laplace law of a stationary drop

Initially, a three-dimensional spherical drop locates at the center of the domain, and a homogenous pressure is set initially. The periodic boundary conditions are imposed at all boundaries. According to the Laplace law, the pressure difference across the drop interface at equilibrium is related to the surface tension

where R is the drop radius, and p _in and p _out are the pressures inside and outside of the drop, respectively. The pressure inside the drop is measured by averaging the pressures at all lattice points inside the spherical drop at 0.7 R , and the outside pressure is measured by averaging the pressures at all lattice points at 1.3 R away from the drop center.

In this numerical computation, we consider four different radii 10, 12, 14, and 16. The corresponding three-dimensional domain is discretized by 6 R × 6 R × 6 R lattice units. In all the cases, the surface tension σ is set to be 0.01, the kinematic viscosity ν is equal to 0.1, and two different density ratios are adopted. Table 1 demonstrates the comparison between the simulating pressure difference Δ p _simulation across the drop interface and the theoretical pressure difference Δ p _theory computing from Eq. ( 20 ). The error for the pressure difference across the interface is within 3.5%. In addition, we take the two-dimensional diffusive model ^{[ 19 ]} as the reference. The maximum error compared to the theoretical solution is 12% for R = 40 lattice units in Ref. [ 19 ], and the maximum error compared to the theoretical solution is less than 4% in the present model. In fact, with increasing drop radius, the error compared to the theoretical value decreases apparently. Therefore, a relatively large drop is required to achieve a reasonable precision in the diffusive model, ^{[ 19 ]} because the interface transition zone in the diffusive model is thicker than that in the present interface tracking method. Finally, the variations of the drop volume are also investigated, the results indicate that the decrement of drop volume in the equilibrium is less than 0.5% compared to the initial volume.

3.2. Oscillation of a three-dimensional ellipsoidal drop

The oscillation of a three-dimensional ellipsoidal drop, which immerses in another fluid, is considered. The frequency of the oscillatory drop is given in Ref. [ 20 ]

where ω _n is the angular frequency, and is Lamb’s natural frequency expressed as

with R _E being the equilibrium radius of the drop. The parameter χ is given by

where μ _o and μ _i are the dynamic viscosities of the two fluids, respectively. The theoretic expression for the period T _theory is obtained from Eq. ( 21 ) through T _theory = 2 π / ω ₂ with n = 2.

Initially, the minimum and the maximum radii of the ellipsoidal drop are 11 and 13, respectively. The ellipsoidal drop is placed in the domain center, and the domain is 72 × 72 × 72 lattice units. Table 2 indicates that the simulated periods coincide well with the theoretic periods T _theory for different viscosities and different surface tensions. The maximum relative error for σ = 0.01 is 1.26% and the maximum relative error for σ = 0.005 is 2.73%. Figure 5 exhibits the interface locations of the oscillating drop as a function of time for three different viscosities. Due to the surface tension and the viscous force, the ellipsoidal drop will reach the equilibrium state with a spherical shape, the spherical radius at equilibrium is R = 11.88. In addition, it takes more time to reach the equilibrium state with a smaller viscosity.

3.3. Drop deformation in a shear flow

In this section, a three-dimensional numerical simulation is performed on the classical Taylor experiment, the drop deformation within a shear flow as illustrated in Fig. 6 . The drop deformation depends on two dimensionless parameters, the Reynolds and the capillary numbers, which are defined as

where γ = U ₀ / H is the shear rate and H is the half height of the parallel plates. For small Re and Ca , the drop eventually tends to be in a steady state, and a parameter D _f is adopted to characterize the drop deformation

where L and B are the largest and the smallest distances from the drop surface to its center (major and minor axes). Taylor ^{[ 21 ]} presented the theoretic result for small D _f

The parameters are adopted as ρ _i = ρ _o = 1, σ = 0.01, Re = 0.1, and the viscosity of the drop is equal to the surrounding fluid. The capillary number varies from 0.05 to 0.30. The computational domain is 12 R × 6 R × 6 R lattice units. Figure 7 indicates that the deformation parameters D _f in this simulation are in good agreement with the theoretic Taylor relation at small capillary numbers. Figure 8 shows that the deformation of the drop is more obvious with a larger capillary number. For Ca = 0.1, the deformation parameter obtained from the simulation is D _f = 0.109, and for Ca = 0.26, the deformation parameter is D _f = 0.309.