{'fraction_non_alphanumeric': 0.03276077806847913, 'fraction_numerical': 0.012598111705152997, 'mean_word_length': 4.646274887604368, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 2, 'zh_pornsignals.json': 0}}
Three-Dimensional Multi-Relaxation Time (MRT) Lattice-Boltzmann Models for Multiphase Flow

In this paper, three-dimensional (3D) multi-relaxation time (MRT) lattice-Boltzmann (LB) models for multiphase flow are presented. In contrast to the Bhatnagar-Gross-Krook (BGK) model, a widely employed kinetic model, in MRT models the rates of relaxation processes owing to collisions of particle populations may be independently adjusted. As a result, the MRT models offer a significant improvement in numerical stability of the LB method for simulating fluids with lower viscosities. We show through the Chapman-Enskog multiscale analysis that the continuum limit behavior of 3D MRT LB models corresponds to that of the macroscopic dynamical equations for multiphase flow. We extend the 3D MRT LB models developed to represent multiphase flow with reduced compressibility effects. The multiphase models are evaluated by verifying the Laplace-Young relation for static drops and the frequency of oscillations of drops. The results show satisfactory agreement with available data and significant gains in numerical stability.


Introduction
In recent years, computational methods based on the lattice-Boltzmann equation (LBE) have attracted much attention. They are based on the paradigm of simulating complex emergent physical phenomena by employing minimal discrete kinetic models that represent the interactions and spatial and temporal evolution of quasi-particles on a lattice [1,2]. Originally developed to overcome certain drawbacks such as the presence of statistical noise and the lack of Galilean invariance of lattice-gas automaton (LGA) [3], the lattice-Boltzmann equation (LBE) [4] has undergone a number of further refinements.
They include enhanced representation of collisions of particle populations through a relaxation process [5] which was further simplified by employing the Bhatnagar-Gross-Krook (BGK) approximation [6] later [7,8]. Also, its formal connection to the Boltzmann equation, developed in the framework of non-equilibrium statistical mechanics, was established [9].
The LBE has the beneficial feature that by incorporating physics at scales smaller than macroscopic scales, complex fluid flows such as multiphase flows can be simulated [10]. In particular, phase segregation and interfacial fluid dynamics can be simulated by incorporating inter-particle potentials [11], concepts based on free energy [12] or the kinetic theory of dense fluids [13,14,15].
The inter-particle potential approach in Ref. [11], is an earlier LBE approach for multiphase flows and is based on non-local pseudo-potentials between particle populations. On the other hand, free-energy based methods are derived from thermodynamic principles that naturally provide interfacial multiphase flow physics [12]. The LBE approaches based on kinetic theory of dense fluids represent physics based on mean-field interaction forces and exclusion volume effects [13,14,15]. More recently, LBE models have also been developed to handle high density ratio problems, which also overcome some of the other limitations of earlier approaches [16,17]. In general, the approaches above employed the BGK approximation to represent the collision process in the LBE and some of their applications in 3D are presented in Ref. [18].
It is well known that the BGK model, a single-relaxation time model, often results in numerical instability when fluids with relatively low viscosities are simulated [19]. The instability problems may be compounded in threedimensional (3D) flows when physics may not be adequately resolved owing to computational constraints. To address this limitation with the standard LB models, several approaches have been proposed. In one approach, known as the entropic lattice Boltzmann method(ELBE), the equilibrium distribution function is defined in such a way that it minimizes a certain convex function, a Lyapunov functional designated as the H−function, under the constraint of local conservation laws [20,21]. As a result, it ensures the positivity of the distribution function of particle populations and thus improves the numerical stability. While this approach is endowed with elegant and desirable physical features, its numerical accuracy is not established and it would incur relatively heavy computational overhead [22].
Alternatives to the BGK model have been proposed to improve numerical stability. In the multi-relaxation time (MRT) method [23], by choosing dif-ferent and carefully separated time scales to represent changes in the various physical processes due to collisions, the stability of the LBE can be significantly improved [19]. Another interesting approach is the fractional time step method [25]. In this approach, stability is improved by considering fractional propagation of particle populations by reducing time step, which increases computational time. While this approach was employed for the BGK model, the underlying idea is not limited to a particular collision model. Indeed such an approach could be employed in conjunction with the MRT model to further improve stability. Yet another possibility for improving stability is based on the so-called two-relaxation time (TRT) models (e.g. Ref. [26]). These are variants of the MRT models and can be considered as limiting cases of the more general MRT models.
In addition to the computational advantage of significantly improved stability, the MRT models also have physical advantages in that they are flexible enough to incorporate additional physics that cannot be naturally represented by the models based on the BGK approximation. In contrast to the BGK models, MRT models deal with the moments of the distribution functions, such as momentum and viscous stresses directly. This moment representation provides a natural and convenient way to express various relaxation processes due to collisions, which often occur at different time scales. Since the BGK model can represent only a single relaxation time for all processes, it cannot naturally represent physics in certain complex fluids. For example, as discussed in Ref. [27], BGK models cannot represent all the essential physics in viscoelastic flows in 3D. Moreover, MRT models make it possible to incorporate appropriate acoustic and thermal properties through adjustable Prandtl numbers for simulating thermo-hydrodynamics by employing a hybrid LBE/finite-difference model [28]. BGK-based models are not able to do so naturally. Also, in a more recent work, MRT models have been shown to be capable of handling more general forms of diffusional transport than possible with the BGK model [30]. In addition to the capability of dealing with genuinely anisotropic diffusion problems, their MRT formulation is able to suppress directional artifacts arising due to discrete lattice effects even in isotropic problems. Such features are possible with MRT models because they lend themselves readily to constructing and controlling various models to evolve at rates consistent with the dynamics of physical phenomena of interest. Indeed, approaches based on the MRT representation are inspired by the moment method developed in the seminal work of Maxwell [31] which was further developed by Grad [32].
While prior works have demonstrated the computational advantages of the MRT models for single-phase flows in 2D by Lallemand and Luo [19] and in 3D by d'Humières et al. [24] and for turbulence modeling [29], it has not been shown for multiphase flows in 3D in previous works. Multiphase flows involve additional physical complexity as a result of interfacial physics involved -i.e., phase segregation and surface tension effects. In this case, the accuracy of the numerical discretization of the source terms representing interfacial physics becomes an important consideration and needs to be incorporated carefully into the framework of the MRT model. These terms should be modeled in a way that, when we establish their relation with the various moments, the dynamical equations in the asymptotic limit should correspond to the desired macroscopic behavior for multiphase flow. These elements were not considered in single-phase MRT models. We introduce a second-order discretization of these source terms and avoid implicitness through a transformation which, to our best knowledge, is applied to MRT models for multiphase flows in 3D for the first time in this work. Recently, a MRT model for 2D multiphase problems was developed [34]. There are significant differences in the development and implementation of 2D and 3D MRT models. These will be discussed below.
Furthermore, realistic multiphase flows are 3D in nature and in this work, we develop a 3D MRT model for such problems by employing a model for interfacial physics based on the kinetic theory of dense fluids [13,14,15].
The specific contributions of this work will now be stated. We develop a 3D MRT model incorporating phase segregation and surface tension effects through source terms discretized by using trapezoidal rule integration and simplified through a transformation to achieve an effectively explicit 3D MRT model for multiphase flows. Since the underlying lattice structure for 2D and 3D models are different, the moment basis for the corresponding MRT models are different, and is more complicated in 3D. We present the theoretical developments based on the 3D moment basis for multiphase flow which, to our best knowledge, has not been considered in prior work. We derive the continuum equations for multiphase flow from the 3D MRT model through a Chapman-Enskog analysis [33], provide the dynamical relations between the various moments and forcing terms representing interfacial physics including surface tension forces, develop the relationships between the gradients of momentum fields in terms of the non-equilibrium parts of certain moments and their corresponding relaxation times, and explicitly derive the relationship between the transport coefficients for the fluid flow and appropriate relaxation times. To improve the stability of the approach, the model is then transformed in such a way that compressibility effects are reduced. Furthermore, we evaluate the accuracy and gains in stability of the MRT model for some canonical multiphase problems in 3D.
As the MRT models are endowed with greater computational stability and potential to incorporate additional physics, there is considerable interest in their applications to multiphase problems: The 2D MRT model developed by McCracken and Abraham [34] extended for axisymmetric problems using the axisymmetric LB model [35] has been employed to study the physics of break up of liquid jets [36]. This axisymmetric MRT model and the 3D MRT model developed in this work have also been employed to study the physics of head-on as well as off-center binary drop collisions, respectively [37]. The physically inspired MRT approach developed in this paper also provides a natural framework for incorporating additional physics such as viscoelastic or thermal effects in multiphase flows through an LB model, which are subjects of future work. The rest of the paper is organized as follows. In Section 2, the 3D MRT LBE multiphase models are developed. In Section 3, the macroscopic dynamical equations of these models are derived by using a Chapman-Enskog multiscale analysis. Section 4 transforms the model to simulate multiphase flow with reduced compressibility effects. In Section 5, the model is applied to benchmark problems to evaluate its accuracy and gains in stability. Finally, the paper closes with summary and conclusions in Section 6.

3D MRT LBE Model for Multiphase Flow
We develop a MRT model for multiphase flows in which the underlying interfacial physics is based on the kinetic theory of dense fluids. In particular, we consider that the particle populations representing the dense fluids experience mean-field interaction forces and respect Enskog effects for dense non-ideal fluids [13,14]. The effect of collisions of particle populations is represented though a generalized relaxation process in which the distribution functions for discrete velocity directions approach their corresponding local equilibrium values at characteristic time scales given in terms of a generalized collision or scattering matrix.
In the following, we consider subscripts with Greek symbols for particle velocity directions and Latin symbols for Cartesian components of spatial directions. We assume summation convention for repeated indices for the components of spatial directions. In addition, unless otherwise stated, we follow the convention that vectors corresponding to three-dimensional position space are represented by non-capitalized symbols with arrowheads; the vectors with a particle velocity basis are denoted by non-capitalized boldface symbols; the square matrices constructed from the particle velocity basis are represented by non-boldface capitalized symbols. We consider the following MRT LBE with a source term that gives rise to phase segregation and surface tension effects: Here, f α is the discrete single-particle distribution function, corresponding to the particle velocity, − → e α , where α is the velocity direction. The Cartesian component of the particle velocity c, is given by On the other hand, the right hand side (RHS) represents the effect of particle collisions and force interactions.
The form of the collision term Ω α that incorporates MRT collision processes will be discussed below. The source term S α which models the interfacial physics in multiphase flow needs to be accurately represented. In Eq. (1), we considered a second-order trapezoidal rule discretization of this term. It may be written as [13] where F I j is the interaction force term that models the interfacial physics, while F ext,j corresponds to external or imposed forces such as gravity. F I j is modeled as a function of density following the work of van der Waals [38].
The exclusion-volume effect of Enskog [33] is also incorporated to account for increase in collision probability due to the increase in the density of non-ideal fluids. These features account for the phase segregation and surface tension effects. For more details, the reader is referred to Ref. [13]. In Eq.
(2), f eq,M α is the local discrete Maxwellian and its functional expression will be presented below. Effectively, F I j may be written as In Eq. (3), F s,j represents the surface tension force and is related to the density ρ and its gradients by where κ is a surface tension parameter. It is related to the surface tension σ of the fluid by the equation [40] σ = κ ∂ρ ∂n where n is the direction normal to the interface. Thus, the surface tension is a function of both the parameter κ and the density profile across the interface.
The term ψ in Eq. (3) refers to the non-ideal part of the equation of state In this work, the Carnahan-Starling-van der Waals EOS [39], is employed where γ = bρ/4. The parameter a is related to the inter-particle pair-wise potential and b is related to the effective diameter d of the particle, and the mass m of a single particle, by b = 2πd 3 /3m. ψ assumes an important role in determining phase segregation [13]. The Carnahan-Starling-van der Waals EOS has a P − 1/ρ − T curve, in which dP/dρ < 0 for certain range of values of ρ, when the fluid temperature is below its critical value. This part of the curve represents an unstable physical situation and is the driving mechanism responsible for keeping the phases segregated and in maintaining a self-generated interface that is diffuse with a thickness of about 3-4 lattice grid points.
The local discrete Maxwellian f eq,M α (ρ, − → u ) in Eq. (2) is obtained from a truncated expansion in terms of fluid velocity − → u of its continuous version, after approximating it for discrete particle velocities − → e α [9]. As a result, it is a function of the fluid densities and velocities and is given by and for the three-dimensional, nineteen-velocity (D3Q19) model [7], shown in The corresponding particle velocity directions are and respectively.
We express the effect of collisions in Eq. (1) as a relaxation process through a multi-relaxation time (MRT) model for the collision term as shown below [19,23,24]: where Λ αβ is the component of the collision or the scattering matrix Λ. Here, the equilibrium distribution f eq β is an appropriate function of the conserved moments such as the density and momentum and is, in general, not necessarily the local discrete Maxwellian given by Eq. (8). The hydrodynamic field variables are obtained by taking the kinetic moments of the distribution functions Equation (1) As a result, Eq. (1) is replaced by the following effectively explicit scheme: where I αβ is the component of the identity matrix I.
Next, we construct appropriate sets of linearly independent moments from the distribution functions in velocity space. Since moments of the distribution function, such as the momentum and viscous stresses, directly represent physical quantities, the moment representation offers a natural and convenient way to express the relaxation processes due to collisions. One particular advantage of this representation is that the time scales of the various processes represented in terms of moments can be controlled independently. The moments are constructed from the distribution function through a transformation matrix T comprising a linearly independent set of vectors, i.e. wherē is the vector representing the particle distribution functions for the D3Q15 model and for the D3Q19 model.
In Eqs. (18) - (20), the superscript ′ T ′ represents transpose of a matrix and the vector f represents the moments and f = ρ, e, e 2 , j x , q x , j y , q y , j z , q z , 3p xx , 3π xx , p ww , for the D3Q15 and D3Q19 models, respectively. For the D3Q15 model, e and e 2 represent kinetic energy that is independent of density and square of energy, respectively, j x , j y and j z are the components of the momentum or mass 3π xx , 3π ww , m x , m y and m z . The first two of these moments have the same symmetry as the diagonal part of the traceless viscous tensor p ij , while the last three vectors are parts of a third rank tensor, with the symmetry of j k p mn [24].
The underlying principle for the construction of the transformation matrix is based on the observation that the collision matrix reduces to a diagonal form in an appropriate orthonormal basis, which is obtained by combinations of monomials of Cartesian components of the particle velocity directions − → e α [23].
Following d'Humieres et al. [24], for the D3Q15 model we introduce a trans-formation matrix T , which represents components of the 15 orthogonal basis column vectors v β that form the moment basis, i.e.
where the components of each column vector may be written as For the D3Q19 model, we have based on 19 orthogonal basis column vectors. The components of each of these These matrices are formed by a set of linearly independent orthogonal basis vectors (i.e. v α · v β = l α δ αβ , where l α is a normalizing constant) which are constructed by a Gram-Schmidt procedure such that they diagonalize the collision matrix Λ i.e. Λ = T ΛT −1 [24], where the matrix Λ in moment space is a diagonal matrix. This orthogonal set is built in increasing order of moments and then arranged in increasing order of complexity of the tensorial representation of the moments. An example of such a construction with some details for a 2D MRT model is given by Bouzidi et al. [41]. As an example, for unit lattice spacing and time steps, i.e. c = 1, the simplified form of the transformation matrices is given in the Appendix.
The equilibrium distribution functions f eq in moment space are related to those in velocity space, f eq , through the transformation matrix f eq = T f eq = ρ, e eq , e 2,eq , j x , q eq x , j y , q eq y , j z , q eq z , 3p eq xx , p eq ww , p eq xy , p eq yz , p eq zx , m eq xyz T (25) for the D3Q15 model and f eq = T f eq = ρ, e eq , e 2,eq , j x , q eq x , j y , q eq y , j z , q eq z , 3p eq xx , 3π eq ww , p eq ww , π eq ww , p eq xy , p eq yz , p eq zx , m eq x , m eq y , m eq given by [24] e eq = −ρ + (j 2 ; p eq xy = jxjy ρ ; p eq yz = jyjz ρ ; p eq xz = jxjz ρ ; m eq xyz = 0 for the D3Q15 model and ; π eq ww = − 1 2 p eq ww ; p eq xy = jxjy ρ ; p eq yz = jyjz ρ ; p eq xz = jxjz ρ ; m eq x = 0; m eq y = 0; m eq z = 0 for the D3Q19 model. The variables s 1 , s 4 , s 6 and s 8 are the relaxation parameters corresponding to the collision invariants ρ, j x , j y and j z , respectively. Since the collision process conserves these particular moments, the choice of their corresponding relaxation times is immaterial. This is because they are each specified to relax during collisions to their local equilibrium which are actually defined to be the corresponding pre-collision value of the respective quantities in the equilibrium distribution f eq (see Eqs. (25) and (26)). So, whatever be the relaxation parameters for these quantities, the collision process does not change their values. The relaxation parameters for these collision invariants could be set to zero when there is no forcing term S α in the LBE. However, care needs to be exercised in selecting the relaxation parameters of MRT LBE with forcing term. This has been pointed out in Ref. [42]. This is because the collision matrix also influences the forcing term in the effectively explicit MRT LBE (see Eq. (17)). To correctly obtain the continuum asymptotic limit of this LBE, i.e. To obtain the LBE model in moment space, the source terms are also transformed as follows: where ς is the column vector corresponding to the components of the source term S α ; ς corresponds to the components of S α . Now, pre-multiplying Eq.
(17) by the transformation matrix T , we obtain the 3D MRT model The hydrodynamic field variables can be obtained from the distribution functions in velocity space as or, more obviously, directly from the components of the moment space vector, f , i.e., f 1 and f β , where β = 4, 6, 8 for density and momentum, respectively.

Macroscopic Dynamical Equations for Multiphase Flow
In this section, the macroscopic dynamical equations for the 3D MRT LBE multiphase flow model developed above will be derived by employing the Chapman-Enskog multiscale analysis [33] for the D3Q15 velocity model. The macroscopic equations for other 3D MRT models, such as the D3Q19 model discussed in the previous section, can be found in an analogous way. Introducing the expansions [43] where ǫ = δ t in Eq. (1), the following equations are obtained as consecutive orders of the parameter ǫ: In equivalent moment space, obtained by multiplying Eqs.(37)-(39) by the transformation matrix T , they become where E i = T e αi T −1 .
First let us simplify the source term S α in Eq.(38) as by neglecting terms of the order of O(Ma 2 ) or higher in its definition, i.e.
Eq.(2) and Eq. (8). In Eq. (43), F x , F y and F z are the Cartesian components of the net force experienced by the particle populations including those that lead to phase segregation and surface tension effects and imposed forces (e.g., Upon substituting Eq. (43), after pre-multiplying by T , in Eq. (41), we get the components of the first-order equations in moment space, i.e. and The evolution equations express variations at the t 0 time scale level, the dynamical relationships between various moments constructed from the velocity space, their relaxation parameters and the net forces that include those that lead to phase segregation and surface tension effects and any external forces acting on the particle populations.
Similarly, the components of the second-order equations in moment space, i.e.
Eq. (42), can be obtained. Our interest is in the dynamical equations for the conserved moments. For brevity, here we express below only the second order equations of the conserved moments representing variations at the t 1 time scale.
and (51) and (62), by using ∂ t = ∂ t 0 + ǫ∂ t 1 , we get and respectively. In these equations, e (1) and p (1) xx , p (1) ww , p (1) xy , p (1) yz and p (1) zx are unknowns to be determined. Writing expressions for these terms from Eqs. (45) and Eqs. and where ζ and ν are the kinematic bulk and shear viscosities repectively. These are related to the relaxation parameters of the energy and the stress tensor moments by 11, 12, 13, 14. (77) Notice that some of the relaxation parameters (s 10 = s 11 = · · · = s 14 ) should be equal to one another to maintain the isotropy of the viscous stress tensor.
From Eqs. (76) and (77), we see that the bulk and shear viscosities can be chosen independently. In particular, this allows us to choose a bulk viscosity larger than the shear viscosity, so that acoustic modes are attenuated quickly. This is helpful in certain physical situations and could also aid in improving numerical stability. Finally, substituting for the forcing terms, i.e. F i = −∂ i ψ+ F s,i + F ext,i , and recognizing that the pressure is given by P = ψ + 1/3ρ, the macroscopic dynamical equations for the conserved moments, i.e. density and momentum, are given by where σ v ij is the deviatoric part of the viscous stress tensor,

3D MRT LBE Model for Multiphase Flow with Reduced Compressibility Effects
Multiphase LBE models based on inter-particle interactions [13] face difficulties for fluids far from the critical point and/or in the presence of external forces [44]. This difficulty is related to the calculation of the inter-particle force involving the term ∂ j ψ, in the model which becomes quite large across interfaces. As a remedy, He and co-workers [14] developed a suitable transformation of the distribution function f α , by invoking the incompressibility condition of the fluid, to g α that determines the hydrodynamic fields, i.e. pressure and ve-locity fields. A separate distribution function is also constructed to capture the interface through an order parameter, which is different from density. The introduction of these features helps in reducing compressibility effects in multiphase problems. In this section, we apply this idea to the 3D MRT model discussed in Section 2. Following He et al. [14], we replace the distribution function f α by another distribution function g α through the transformation By considering the fluid to be incompressible, i.e.
and using Eqs. (83) and (16), Eq. (1) becomes where g eq α = f eq α RT + ψ(ρ) The corresponding source terms are Equations (85)-(87) are numerically more stable than the original formulation since the term ∂ j ψ is multiplied by a factor proportional to the Mach number O(Ma), instead of being O(1) when the "incompressible" transformation is not employed. As a result, this term becomes smaller in the incompressible limit. Now, applying the transformation matrix to the above system, we get In this new framework, we still need to introduce an order parameter to capture interfaces. Here, we employ a function φ referred to henceforth as the index function, in place of the density as the order parameter. The evolution equation of the distribution function, whose emergent dynamics govern the index function, must have the term responsible for maintaining phase segregation and mass conservation. In this regard, we use Eqs. (1) and (2), by keeping the term involving ∂ j ψ and dropping the rest as they play no role in mass conservation. In addition, the density is replaced by the index function in these equations. Hence, the evolution of the distribution function for the index function is given bȳ where In moment space, this becomes It follows that hydrodynamical variables, such as the pressure and fluid velocity, can be obtained by taking appropriate kinetic moments of the distribution function g α , i.e.
The index function is obtained from the distribution functionf α by taking the zeroth kinetic moment, i.e.
The density is obtained from the index function through linear interpolation, i.e. As noted in Ref. [24], suitable code optimization techniques should be applied to avoid substantial increase in computational time due to the use of the MRT model. First, direct matrix computations, which can lead to increased computational overhead, should not be carried out in the transformation between velocity and moment spaces. Instead, the transformation between the components of the distribution functions in velocity space f α and the moments f β should be carried out explicitly using the expressions obtained from the transformation matrix, which is an integer matrix with several common elements and many zeroes. In these calculations, all the common sub-expressions should be computed just once. Second, the computational algorithm should be carried out as discussed in the previous paragraph. With such optimization strategies, the MRT model only incurs a moderate increase in computational time, about 35% more than the BGK model for the cases discussed in the next section, but with significantly improved stability.

Results and Discussion
The 3D MRT model developed in the previous section will now be evaluated for some benchmark multiphase problems. Unless otherwise specified, we express the results in lattice units, i.e., the velocities are scaled by the particle velocity c and the distance by the lattice spacing δ x . We first consider a static multiphase problem, namely, the verification of the well-known Laplace-Young  [13,14] is stable for only low density ratios. The MRT model has relatively little influence on this limitation. However, as will be shown below, it can significantly increase the stability at lower viscosities.
It may be noted that when the BGK model is employed to simulate multiphase flow with the viscosity used above, the computations become unstable. As  sults that have negligible difference between them. Henceforth, for comparison purposes we will only consider results obtained with the D3Q15 lattice. The density profile across the drop interface is found to be characterized well with no anisotropic effects. As with other LB models, and with methods based on the direct solution of Navier-stokes equations for multiphase flows, velocity currents around the interfaces are observed. These are generally small and is found to be proportional to the surface tension parameter κ.
Next, we consider a dynamical problem, namely, the oscillation of a liquid drop immersed in a gas. We employ the analytical solution of Miller and Scriven (1968) [46] for comparison with the computed time periods. According to Ref. [46], the frequency of the n th mode of oscillation of a drop is given by where ω n is the angular response frequency, and ω * n is Lamb's natural resonance frequency expressed as [47] (ω * n ) 2 = n(n + 1)(n − 1)(n + 2) R d is the equilibrium radius of the drop, σ is the interfacial surface tension, and ρ l and ρ g are the densities of the two fluids. The parameter χ is given by where µ l and µ g are the dynamic viscosities of the two fluids. Here, the subscripts g and l refer to the gas and liquid phases, respectively. We consider the second mode of oscillation, i.e. n = 2, and analytical expression for the time period T anal is obtained from Eq. (96) through T anal = 2π/ω 2 .
The initial computational setup consists of a prolate spheroid with a minimum and maximum radii of 11 and 15, respectively, and placed in the center of a domain discretized by 41 × 41 × 41 lattices. We consider a density ratio of 4 for different values of shear viscosity and surface tension parameters. Figure 4 shows the drop configurations at different times for ν = 1.6667 × 10 −2 and κ = 0.10. At t = 100, the drop has a prolate shape which temporarily assumes an approximate spherical shape at t = 500. At t = 700, notice that the    Figure 7 shows the interface locations when κ = 0.12. The 3D MRT model remains as stable as in the previous cases even for more complex dynamical problems such as when drops undergo shearing and other types of motions. In a recent study, this model was applied to simulate collision of a pair of drops under different conditions [37]. In particular, when employed to study off-center collisions it is able to reproduce experimentally observed [48] complex interfacial shape changes due to such dynamical events as shearing, rotation, stretching deformation and breakup, while maintaining an improved stability compared to the BGK model. Indeed, in some cases reported in Ref. [37] as much as an order of magnitude improvement in stability was observed by the use of MRT model in lieu of the BGK model.

Summary and Conclusions
In this paper, we develop 3D MRT LB models for multiphase flows. They employ the LBE with a generalized collision term together with forcing terms representing interfacial physics. By considering 3D lattice velocity models and a second-order discretization of the forcing terms, the LBE written in terms of the distribution function is transformed into an equivalent system represented by a set of conserved and non-conserved moments. The conserved or hydrodynamic moments include the fluid density and momentum, while the non-conserved or kinetic moments include heat fluxes and viscous stresses.
The models are developed such that when the collision term in the LBE is written in moment space, the moments relax to their equilibrium values at rates that can be adjusted independently.
By applying the Chapman-Enskog multiscale analysis to the MRT models, we show that they correctly recover the 3D hydrodynamical equations for multiphase flows in the continuum limit. The dynamical relationships between various moments and the forcing terms, including surface tension forces, are systematically derived. Some of the relaxation parameters in the collision term are shown to be related to the shear and bulk kinematic viscosities of the fluid.
The ability to independently adjust the relaxation parameters enhances the stability of the MRT models. The models are evaluated for accuracy by solving some multiphase test problems. It is shown that the 3D MRT models verify the Laplace-Young relation for static drops to within 8%. Computations with MRT models of drop oscillations show that shear viscosities can be lowered by a factor of 5 when compared to the BGK model. The computed time period of oscillations agrees with the analytical solution to within 5%. The MRT model, when applied to more complex multiphase problems such as binary drop collisions, is found to be as stable as that for these simple canonical problems. The MRT approach developed in this work can be extended to LBE multiphase flow models such as those developed in Refs. [16,17] to handle high density ratio problems. Moreover, it also provides a more natural framework to extend the LBE for more general situations such as viscoelastic or thermal effects in multiphase flows in future investigations.

A Appendix. Simplified Transformation Matrices for Unit Lattice
Spacing For unit lattice spacing and time steps, i.e. c = 1, the transformation matrix, Eq. (23), for the D3Q15 model simplifies to and that for the D3Q19 model to