Lattice Boltzmann modeling and simulation of liquid jet breakup

A three-dimensional color-fluid lattice Boltzmann model for immiscible two-phase flows is developed in the framework of a three-dimensional 27-velocity (D3Q27) lattice. The collision operator comprises the D3Q27 versions of three sub-operators: a multiple-relaxation-time (MRT) collision operator, a generalized Liu--Valocchi--Kang perturbation operator, and a Latva-Kokko--Rothman recoloring operator. A D3Q27 version of an enhanced equilibrium distribution function is also incorporated into this model to improve the Galilean invariance. Three types of numerical tests, namely, a static droplet, an oscillating droplet, and the Rayleigh--Taylor instability, show a good agreement with analytical solutions and numerical simulations. Following these numerical tests, this model is applied to liquid-jet-breakup simulations. The simulation conditions are matched to the conditions of the previous experiments. In this case, numerical stability is maintained throughout the simulation, although the kinematic viscosity for the continuous phase is set as low as $1.8\times10^{-4}$, in which case the corresponding Reynolds number is $3.4\times10^{3}$; the developed lattice Boltzmann model based on the D3Q27 lattice enables us to perform the simulation with parameters directly matched to the experiments. The jet's liquid column transitions from an asymmetrical to an axisymmetrical shape, and entrainment occurs from the side of the jet. The measured time history of the jet's leading-edge position shows a good agreement with the experiments. Finally, the reproducibility of the regime map for liquid-liquid systems is assessed. The present lattice Boltzmann simulations well reproduce the characteristics of predicted regimes, including varicose breakup, sinuous breakup, and atomization.


I. INTRODUCTION
Multiphase and multicomponent flows appear in many natural and industrial processes.A liquid jet injected into another fluid is an interesting example of such a flow.Significant efforts have been put into understanding the breakup of a liquid jet for more than a century.Since the pioneering works of Plateau [1] and Rayleigh [2], extensive studies on this subject have been performed theoretically, experimentally, and numerically [3][4][5][6].
Drops form directly from the nozzle at low injection velocities, and a liquid jet issues from the nozzle and then breaks into droplets in various patterns at higher injection velocities.The occurrence of such a regime is of interest in the study of liquid-jet breakup.Ohnesorge [7] classified his results into four types of breakup regimes: (0) dripping, (I) varicose, (II) sinuous, and (III) atomization [8,9].He also provided a regime map of liquid jets in a gas using the Ohnesorge and Reynolds numbers.After Ohnesorge's work, much research on this subject has been performed (e.g., [10][11][12]).The majority of investigations have focused upon liquid-gas systems (liquid jet into gaseous atmosphere).Breakup of jets in liquidliquid systems (liquid jets into another liquid) has not been investigated as extensively.Our focus in this paper is therefore on the breakup of liquid jets in immiscible liquid-liquid systems.
The liquid-liquid-jet systems can also be found in several fields, e.g., chemical processing [13][14][15] and CO 2 storage in oceans [16,17].In the field of nuclear engineering, interactions between melt and coolant must be well understood for safety design of nuclear reactors and have therefore been extensively investigated in the literature [18][19][20].In experiments, a high-temperature melt is often used to simulate the core melt materials.Abe et al. [20] discussed the relationships between the fragment size and the Rayleigh-Taylor and Kelvin-Helmholtz interfacial instabilities [21].
To better understand the fundamental interactions between melt-jet and coolant interactions, experiments using appropriate test fluids under isothermal conditions are also effective approaches as a separate effect of such interactions [19,22].Saito et al. [22] developed a breakup-regime map for a jet in immiscible liquidliquid systems based on experiments and phenomenological considerations.Figure 1 shows the dimensionless map and the corresponding visual images.The Ohnesorge classification [7] was extended to liquid-liquid systems.As can be seen, various flow regimes occur during liquidjet-breakup processes.The breakup transitions from the dropping regime to the atomization regime, and the generated droplet size drastically changes depending on the arXiv:1705.03141v2[physics.flu-dyn]28 Jun 2017 conditions.This implies that the breakups of liquid jets are essentially three-dimensional flows and possess multiscale phenomena such as droplet pinch-off and atomization.
Numerical simulations of liquid-liquid jets involve the solution of the Navier-Stokes equations for two fluids with specified boundary and interface conditions.Several approaches to solving these types of free-surface problems are available in the literature.As a first attempt, Richards et al. [23] investigated the axisymmetric steady-state laminar jet based on the volume-of-fluid (VOF) method [24].Thakre et al. [25] also used the VOF method provided by a commercial code, FLUENT, to simulate a melt jet into water.They successfully reproduced a variation in the breakup length [26]; later, a similar variation was confirmed by experiments [22,27].Homma et al. [28] numerically investigated liquid-liquidjet breakup using a front-tracking method [29,30].They mapped different breakup modes on a plot of Weber number vs. viscosity ratio.The drawback of their fronttracking simulations was that they neglected the coalescence of generated drops.The aforementioned numerical simulations [23,25,28] were limited to two-dimensional cases.
A completely different approach is the use of a lattice Boltzmann method.Several authors have investigated liquid-liquid-jet flows using multiphase lattice Boltzmann models [31][32][33].In recent years, this method has been recognized as a powerful tool for analysis of complex fluid dynamics, including multicomponent and multiphase flows [34].Compared with other macroscopic CFD methods based on the Navier-Stokes equations, the lattice Boltzmann method, which is constructed using mesoscopic kinetic equations, has several advantages.For instance, it is easy to incorporate mesoscale physics such as interfacial breakup or coalescence.Moreover, the computational cost for simulating realistic fluid flows is reasonable compared with particle-based methods (e.g., molecular dynamics).The relations of the scale properties in fluid flows are schematically illustrated in Fig. 2.
Two-phase or multiphase lattice Boltzmann models can be classified into four categories: • Color-fluid model [35,36], • Pseudo-potential model [37,38], • Free-energy model [39,40], • Mean-field model [41].This classification may not be exhaustive, for instance, the latter two models are sometimes identified as phasefield models [42] since the Cahn-Hilliard or similar interface tracking equations can be derived from them.For details about the multiphase lattice Boltzmann models, interested readers can refer to the comprehensive review papers [34,[42][43][44][45][46] and references therein.In this paper, we focus upon the color-fluid model.This model possesses many strengths in simulations of multiphase/multicomponent flows, including strict mass conservation for each fluid and flexibility in adjusting the interfacial tension [47].A static drop test is no longer needed to determine the interfacial tension; it can be directly obtained without any analysis or assumptions.Moreover, the color-fluid model shows a very small dissolution property for tiny droplets or bubbles [46].
Color-fluid models, which are often referred to as R-K or color-gradient models, were first developed by Gunstensen et al. [35], who extended the two-component lattice gas automata model of Rothman and Keller [48].Later, Grunau et al. [36] enabled the introduction of density and viscosity ratios by modifying the forms of the distribution functions.Latva-Kokko and Rothman [49] replaced Gunstensen's maximization-recoloring step with a formulaic segregation algorithm.Instead of widening the interface width, Latva-Kokko-Rothman's recoloring algorithm solves some issues with the previous color-fluidtype model, namely, the lattice-pinning problem and the spurious velocities.Reis and Phillips [50] extended the model to common a two-dimensional nine-velocity (D2Q9) lattice.They modified the perturbation operator to recover the Navier-Stokes equations correctly.
Leclaire et al. [51] demonstrated that integrating Latva-Kokko-Rothman's recoloring operator [49] into Reis-Phillips' perturbation operator [50] greatly improves the numerical stability and accuracy of solutions over a wide range of parameters.Using an isotropic gradient operator also enhanced numerical stability and accuracy [52].Liu et al. [53] derived a generalized perturbation operator using the phase-field (or order parameter) instead of a color-gradient and formulated the colorfluid model in three dimensions.Very recently, Leclaire et al. [54] generalized the color-fluid-type lattice Boltzmann model in two and three dimensions.
Galilean invariance is one of the issues to be improved in the color-fluid family.Following Holdych et al. [55], a source term to improve the Galilean invariance was derived by Leclaire et al. [56] and incorporated into an equilibrium distribution function.The enhanced equilibrium distribution function showed an improvement of the momentum-discontinuity problem through numerical tests on a layered Couette flow.Recently, Ba et al. [47] have modified an equilibrium distribution function based on the third-order Hermite expansion of the Maxwellian distribution.They also showed that discontinuous velocities were improved by this modification.
It is known that the LB method suffers from numerical instability in low-viscosity conditions.Modification of the collision operator is one method for overcoming this issue [57].A multiple-relaxation-time (MRT) collision operator or generalized lattice Boltzmann equation [58][59][60] has been widely used, even for multiphase flows, to enhance numerical stability and accuracy and to reduce spurious current near the interface.
We return to the issue of lattice Boltzmann simulations of liquid-jet breakup.McCracken and Abraham [61] successfully introduced an MRT operator to the multiphase lattice Boltzmann model and performed liquidjet breakup simulations [31].They assumed that the flow was axisymmetric in two dimensions.They investigated the influence of interfacial tension, injection velocity, and liquid viscosity under a density ratio of 5.However, three-dimensional simulations are required to further understand breakup characteristics, since liquidjet breakup is an essentially three-dimensional flow, as shown in Fig. 1.
The authors carried out lattice Boltzmann simulations of liquid-jet breakup in three dimensions [32,33].Matsuo et al. [32] used the three-dimensional two-phase lattice Boltzmann model, which was developed by Ebihara and Watanabe [62] based on the model of He et al. [41].They compared their simulation results with experiments and investigated the effect of the Froude number upon the jet-breakup length.In their simulations, however, the Reynolds number was limited to O(10 2 ).This was an order of magnitude smaller than that in the target experiments.In addition, the model of He et al. [41] suffered from the dissolution of tiny droplets [46]; thus, it would not be appropriate for liquid-jet-breakup simula-  [53] and applied this model to liquid-jet-breakup simulations.Although they could simulate liquid-jet breakup with the Reynolds number up to O(10 3 ), the kinematic-viscosity ratio was set to unity to avoid numerical instability.This meant that the kinematic viscosity of the surrounding liquid in their simulation was more viscous than that in the target experiment.Further improvement is required to compare the numerical results with experiments; this is the motivation of the present study.

Navier
In this paper, we present the three-dimensional twophase lattice Boltzmann model for immiscible two-phase flows and its application to liquid-jet breakup.In Sec.II, we formulate the three-dimensional two-phase lattice Boltzmann model for immiscible two-phase flows in the framework of a three-dimensional 27-velocity (D3Q27) lattice.The collision operator consists of D3Q27 versions of three sub-operators: an MRT-collision operator, a generalized-perturbation operator [53], and a formulaticrecoloring operator [49].A D3Q27 version of the enhanced equilibrium distribution functions [56] is also incorporated into this model to improve its Galilean invariance.In Sec.III, numerical tests, including those of a static droplet, an oscillating droplet, and the Rayleigh-Taylor instability, are used to validate the developed model.In Sec.IV, this model is applied to liquid-jetbreakup simulations.A simulation in which the parameters are exactly matched to the target experiment is performed and compared with experimental data.Finally, we assess the reproducibility of the breakup regimes expected by the dimensionless-regime map [22].Sec.V concludes this paper.

II. METHODOLOGY
The present model is formulated on a D3Q27 lattice.The key to the formulation is the combination of previous work and their extension to the D3Q27 framework.The main points of this process can be briefly summarized as follows: • Introducing a D3Q27 MRT collision operator [63] and relaxation parameters [64] • Extending an enhanced equilibrium distribution function [56] to the D3Q27 lattice • Extending a generalized perturbation operator [53] to the D3Q27 lattice For the present three-dimensional lattice Boltzmann model, distribution functions move on a D3Q27 lattice with the lattice velocity c i defined as follows: where c = δ x /δ t , δ x is the lattice spacing and δ t is the time step.The schematic structure of the D3Q27 lattice is shown in Fig. 3.The D3Q27 model is a straightforward extension of the D2Q9 model [65].
In this model, two immiscible fluids are represented as pseudo red and blue fluids, respectively.The distribution function, f k i , is introduced to represent the fluid k, where k = r and b denote the colors "red" and "blue," respectively, and i is the lattice-velocity direction.The total distribution function is defined as The evolution of the distribution function is expressed by the following lattice Boltzmann equation: where x and t are the position and time, respectively.The collision operator Ω k i is made up of three suboperators [66]: where (Ω k i ) (1) is the single-phase collision operator, (Ω k i ) (2) is the perturbation operator, and (Ω k i ) (3) is the recoloring operator.Using the MRT operator, the singlephase collision operator can be written as The density of the fluid k is given by The total fluid density is given by ρ = k ρ k , and the total momentum is defined as where F is the body force.Note that Eq. ( 6) indicates that the local velocity is modified to incorporate the spatially varying body force [67].In Eq. ( 2), M and K are, respectively, the 27 × 27 transformation and relaxation matrices.In this paper, Dirac's bra-ket notation is employed.Here, the "bra" operator f | denotes a row vector along each lattice-velocity direction, i.e., (f 1 , f 2 , . . ., f 27 ), and the "ket" operator |f denotes a column vector, i.e., (f 1 , f 2 , . . ., f 27 ) T , where the superscript "T" is the transpose operator.
The MRT collision operator in three dimensions is usually implemented with the D3Q15 or D3Q19 lattices.As an early attempt for the D3Q27 MRT collision operator, Dubois and Lallemand [68] and Premnath and Banerjee [69] presented the formulation independently in 2011.Dubois and Lallemand [68] arranged the D3Q27 orthogonal basis vectors for the moments according to their character (scalars, vectors, and tensors, etc.) and then used raw moments to formulate an MRT lattice Boltzmann model.On the other hand, Premnath and Banerjee [69] arranged the D3Q27 orthogonal basis vectors based on the increasing order of moments and used central moments to formulate an MRT lattice Boltzmann model.
Later, Geier et al. [63] provided the following orthogonal moment set for the D3Q27 lattice: where The transformation matrix transfers the distribution functions from velocity space to moment space.Using the MRT collision operator instead of a traditional single-relaxation-time (or BGK) collision operator contributes to enhancement of numerical stability and accuracy, even with additional computational costs.The practical forms of the transformation matrix and its inverse are given in Appendix A. The relaxation matrix K is the diagonal matrix given by [64] where the elements 0 < s i < 2 represent both the hydrodynamic and non-hydrodynamic relaxation parameters.The hydrodynamic parameters satisfy the following relations: where ν and ζ are the kinematic viscosity and the bulk viscosity, respectively.In this paper, we use the optimized parameters proposed by Suga et al. [64]: s 1 = 0, s 5 = 1.5, s 11 = 1.5, s 14 = 1.83, s 17 = 1.4,s 18 = 1.61, s 19 = s 21 = 1.98, and s 24 = s 27 = 1.74.We confirmed that these parameters significantly enhanced the numerical stability even at extremely low kinematic viscosity with the order of O(10 −4 ).
For the single-phase collision operator, an enhanced equilibrium distribution function proposed by Leclaire et al. [56] is used in this paper: In the case of Φ k i = 0, Eq. ( 11) recovers the common form of an equilibrium distribution function.Using the form of Eq. ( 11), the Galilean invariance is improved for variable density and viscosity ratios under the hypothesis of a small pressure gradient [56,70,71].However, it should be noted that this can only partly restore Galilean invariance, as still the third order diagonal equilibrium moments are not independently supported, and are related to the corresponding first moments.Consequently, there will be cubic velocity errors in Galilean invariance even for the D3Q27 lattice, which could become perceptible for flows under high shear.This issue can be solved by making some additional corrections to the collision operator, as suggested recently by Dellar [72].
The weights, w i , are those of a standard D3Q27 lattice [65]: Moreover, in the D3Q27 lattice, one can derive and where ⊗ is the tensor product and the symbol ":" indicates tensor contraction; ν is the kinematic viscosity interpolated by [53,66] Here, φ is the order parameter to distinguish the two components in a multicomponent flow, defined as [53] φ The values of the order parameter φ = 1, −1, and 0 correspond to a purely red fluid, a purely blue fluid, and the interface, respectively [66].In the framework of D3Q27 lattice, the tensor G k in Eq. ( 14) is defined as As established in [36], the density ratio between the fluids, γ, must be taken into account as follows to obtain a stable interface: where the superscript "0" indicates the initial value of the density at the beginning of the simulation [56].In each homogeneous phase region, the pressure of the fluid k is given by for the D3Q27 lattice.This corresponds to an isothermal equation of state.In this paper, we choose α b = 8/27, in which c b s = 1/ √ 3 [33,70].The term |F in Eq. ( 4) represents the discrete forcing term accounting for the body force F. In the MRT framework, the forcing term reads as [73] where Eqs. (20) and ( 21) reduce to Guo et al.'s original forcing scheme [67] when using a single-relaxation time [73] To model the interfacial tension, Liu et al. [53] derived a generalized perturbation operator based on the CSF [74] concept, and the work of Reis and Phillips [50] is employed to obtain the interfacial tension: Eq. ( 22) satisfies the correct form of the interfacialtension force in the Navier-Stokes equations when the lattice-specific variables B i are chosen correctly.We derived the values of B i in the framework of the D3Q27 lattice as follows: In this model, the interfacial tension can be directly given by where we assumed that A = A r = A b , τ is the relaxation time.Parameter A controls the strength of interfacial tension, σ.
Although the perturbation operator, (Ω k i ) (2) , generates interfacial tension, it does not guarantee the immiscibility of both fluids.To promote phase segregation and maintain the interface, the following recoloring operator is applied [49,51,75] ( ( where (e) i is evaluated using Eq. ( 11), a zero velocity, and In the present model, the following continuity and Navier-Stokes equations can be derived via Chapman-Enskog analysis [53,67,76] where is the viscous stress tensor with the shear viscosity ν given by Eq. ( 9) and the bulk viscosity ζ given by Eq. ( 10); p = p r + p b is the pressure.In Eq. ( 30), the term ∇ • S arises from the perturbation operator given by Eq. ( 22) and is equivalent to the interfacial force based on the CSF concept [53]; the capillary stress tensor, S, is given by For the computation of the gradient operator for an arbitrary function χ, the following second-order isotropic central scheme [53,[77][78][79] is adopted.
All the simulations in this paper are carried out in lattice units.In this paper, δ x and δ t are set to 1 as in the usual lattice Boltzmann simulations.The aforementioned formulation focuses on the two-component systems.It should be straightforward to implement the present model in three or more component systems according to the work of Leclaire et al. [80].
One can consider the key to the developed lattice Boltzmann model in this paper to be based on a combination of previous works and their extension to the framework of the D3Q27 lattice.A brief procedure to determine the related lattice-specific coefficients [Eqs.(13), (14), and (23)] is described in Appendix B.

III. NUMERICAL TESTS A. Static droplet
Static-droplet tests are performed to test the validity of interfacial tension predicted by Eq. ( 24).The computational domain is discretized into an 85 × 85 × 85 lattice.A steady red droplet with radius R is immersed in a blue fluid.The density field of each phase is initialized as follows: where 2 with the central position of the computational domain (x c , y c , z c ).We set the density and kinematic-viscosity ratios as 1.5 and 1, respectively.The kinematic viscosity for each phase is set to be 1.0 × 10 −3 , and gravity is neglected throughout the simulations.The parameters are A in Eq. ( 24) and initial droplet radius R (see Table I).Periodic boundary conditions are imposed on all sides of the computational domain.
The Laplace equation in three dimensions is given by where ∆p is the pressure difference across the droplet interface.The pressure for each phase is evaluated by Eq. ( 19) and is measured after 80,000 iterations using the same procedure as Leclaire et al. [51].Figure 4 shows the measured pressure differences.For both the higher [Fig.4(a)] and lower [Fig.4(b)] interfacial tensions, the results are in proportion to the droplet curvature 1/R, with which the Laplace law was satisfied.The theoretical prediction by Eq. ( 24), shown in Fig. 4 as a solid line, also agrees well with the measured pressure differences.Table I summarizes the simulation parameters and the evaluated error.The error E is calculated as [53] where σ th and σ Lap are the interfaceial tension predicted by Eq. ( 24) and that measured by the Laplace equation Eq. ( 36), respectively.We confirm that the lattice Boltzmann model developed in Sec.II can predict the interfacial tension for a static case within a maximum error of 1.6%.We should mention the influence of lattice isotropy on the so-called spurious velocity.Fig. 5 compares the droplet shape and velocity field at the equilibrium state for D3Q19 [Fig.5(a)] and D3Q27 lattices [Fig.5(b)].The numerical test using D3Q19 is based on [33].The maximum spurious velocities |u| max are 1.2 × 10 −2 for D3Q19 and 5.8 × 10 −3 for D3Q27, respectively.Although the conditions are same except for the employed lattice geometry, the simulation using the D3Q27 lattice shows better result.Enhancing the lattice isotropy (from D3Q19 to D3Q27) contributes to reducing the spurious velocity.24).The present simulations show a good agreement with the theoretical prediction.

B. Oscillating droplet
Oscillatory-droplet tests are performed to test the model validity in an unsteady case.The simulation setup and parameters for the present numerical tests follow Premnath and Abraham [81], except for the interfacialtension coefficient.A droplet with a prolate-spheroid shape is placed at the center of a computational domain discretized by a 41 × 41 × 41 lattice.The droplet's minimum and maximum radii are 11 and 15, respectively.We set a density ratio of 4. The kinematic viscosity ratio is set to be unity (ν = ν r = ν b ), and gravity is neglected throughout the simulations.The parameters are ν and A (see Table II).Parameter A is the same as the ones provided in Sec.III A.
The analytical solution of Miller and Scriven [82] is used for comparison with the computed time periods.The frequency of the n-th mode is given by where ω n is Lamb's natural resonance frequency [83] ( Note that R e in Eq. ( 39) is the equivalent radius of a droplet.Subscripts 1 and 2 refer to the dispersed and continuous phases, respectively, so they can be replaced by r and b in this paper.Parameter χ is given by We are only interested in the second mode (n = 2) here.The analytical expression for the time period is obtained by T th = 2π/ω 2 .
Figure 6 shows the transient change of the oscillating droplet shapes with ν = 3.333 × 10 −3 and A = 1.0 × 10 −2 .After assuming a spherical shape at t = 700, the droplet becomes an oblate spheroid at t = 1, 400.The configuration turns into a prolate spheroid at t = 2, 500.Such a series of oscillations can be also seen in Ref. [81] The interfacial location is measured, and the results are shown in Fig. 7 as a function of time.The interfacial locations are recorded per 10 time steps.Under higher-interfacial tension [Fig.7(a)], all cases attenuate with oscillation regardless of the kinematic viscosity.The higher the kinematic viscosity, the earlier attenuation occurs.For the lower-interfacial-tension case [Fig.7(b)], the timescale of oscillation becomes qualitatively longer; this tendency agrees with the theoretical expectation by Eqs.(38) and (39).Only in the case of ν = 1.667 × 10 −2 does the interfacial location reach an equilibrium spherical shape without a series of oscillations, as shown in Fig. 6.The effect of viscous damping distinguishes rather than the effect of interfacial tension in this case.
The simulation parameters and evaluated errors in the oscillating period are summarized in Table II.The error E is calculated as where T sim is the measured oscillation period.Note that, in the case of ν = 1.667 × 10

C. Rayleigh-Taylor instability
To assess the validity of the body force implementation [Eqs.(20) and ( 21)], the Rayleigh-Taylor instability is selected as the next numerical test.The Rayleigh-Taylor instability is a fundamental interfacial instability that is induced when a heavy fluid is placed over a light one subjected to a slightly disturbed interface in a gravitational field [21].The Rayleigh-Taylor instability has received considerable attention owing to its extensive applications, e.g., in the fundamental process of melt jet breakup [20,22].To our knowledge, this is the first time a color-fluid model has been applied to the threedimensional Rayleigh-Taylor instability.
We refer to the computational setup adopted in He et al. [84].A schematic diagram of the computational setup is illustrated in Fig. 8.The top and bottom boundaries are no-slip walls; the lateral boundaries are periodic.
As in the manner of [84], the single-mode initial perturbation is initially imposed as at the mid-plane, where W is the width of the computational domain.The body force in Eq. ( 21) for this problem is incorporated as with g = (0, 0, −g).The gravitational acceleration g is chosen to satisfy (W g) 1/2 = 0.04 [41].The computational domain is set to be W × W × 4W with W = 128.The domain's size resulted in a 128 × 128 × 512 lattice.The Atwood number is fixed at 0.5 throughout the simulations.The interface tension is neglected; thus, the perturbed interface is expected to always be unstable in the inviscid case.The kinematic-viscosity ratio is set to unity.Another dimensionless group is the Reynolds number, defined as We use three patterns of Reynolds numbers: Re = 512, 1, 024, and 5, 120.
Figure 9 shows the interfacial evolution of the Rayleigh-Taylor instability.The dimensionless time is given by In the initial stages by t * = 2, the Reynolds number dependence on the interfacial configuration is small.We can see the edge of the spike rolled up at t * = 3.At later stages (t * = 3, 4), the higher the Reynolds number is, the more unstable the interface becomes.For the Re = 5, 120 case [Fig.9(c)], the interface becomes especially complicated.The Kelvin-Helmholtz instability appears conspicuously as the Reynolds number increases.(44)] is fixed at 0.5 in all simulations.At later stages, the higher the Reynolds number is, the more unstable the interface becomes.
The time history of the positions of the bubble front, spike tip, and saddle point are calculated and plotted in Fig. 10(a) (see Fig. 9 for the locations of the bubble front, spike tip, and saddle point).As can be seen in Fig. 10(a), the differences in the positions are very small in our simulations, unlike the interface structure shown in Fig. 9.He et al. [84] pointed out that, when Re > 512, the Reynolds number dependence is negligible.The present simulation, using the forcing scheme in Eqs. ( 20) and ( 21), shows a similar trend, supporting their suggestion.
Using the same computational setup and the parameters of He et al. [84], three-dimensional Rayleigh-Taylor simulations were carried out using the lattice Boltzmann [85] and phase-field methods [86].Wang et al.
[85] used the phase-field-based MRT lattice Boltzmann model; Lee and Kim [86] directly solved the Navier-Stokes-Cahn-Hilliard equations.Here, we compare our results with previous works [84][85][86].Figure 10(b) shows the time histories of the positions at Re = 1, 024 and At = 0.5.Comparing the results, the time changes of the bubble front and the saddle point are almost the same, irrespective of the method used.For the spiketip evolution, no significant difference is observed among the lattice Boltzmann simulations.At later stages, these simulations are found to penetrate more deeply than the phase-field simulation.The difference between the lattice Boltzmann and phase-field methods is considered to arise from the different wall-boundary conditions implemented.

A. Setup
Figure 11 illustrates a schematic diagram of the boundary conditions for liquid-jet-breakup simulations.This computational setup is the same as that of Saito et al. [33], except for the outflow boundary.In the initial state, the computational domain is filled with blue-particledistribution functions, f b i , with zero velocity.The boundaries consist of an inflow boundary, a wall boundary, and an outflow boundary.A circular inflow boundary is implemented at the top within (x − x c ) 2 + (y − y c ) 2 < (D j0 /2) 2 , where x c and y c are central locations on the x-y plane.The uniform velocity u j0 is implemented, with the corresponding equilibrium functions given at this site.No artificial disturbances are considered at the inflow boundary.A wall boundary is implemented on the top (except for an inflow-boundary site) and on the lateral sites.A free-slip condition [87] is implemented as a wall boundary condition.At the outflow, a convective boundary condition [88] is used, in which the following convective equation for the distribution functions is solved, where N is the node on the outflow boundary.Following Lou et al. [88], we set two ghost-nodes N + 1 and N + 2. The discretized form can be given by the first-order implicit scheme, where λ = U c (t + δ t )δ t /δ x .For the convective velocity normal to the outflow boundary U c , several choices can be considered, e.g., the local velocity, the average velocity, and maximum velocity [89].Through some numerical tests, we determine that the local velocity is suitable to the present system, that is, where u z (x, t) = u z (x, y, z, t) is the z-direction component of fluid velocity u.
The body force in Eq. ( 21) for the liquid-jet simulations is set as with g = (0, 0, g).Eq. ( 50) means that the gravitational force acts only in the dispersed phase [45].

B. Comparison with experiments
Using the lattice Boltzmann model developed in Sec.II and the boundary conditions provided in Sec.IV A, we perform here numerical simulations of liquid-jet breakup.The parameters for this simulation are set to be exactly the same as in the experiments [33], and the results are compared.In the target experiments, a glycerin-watermixture jet was injected into a silicon-oil pool.These test fluids were immiscible with each other.In [33], the lattice Boltzmann simulation was also carried out using the MRT color-fluid model based on the D3Q19 lattice.
In the framework of linear theory, one can choose the following dimensionless groups to describe the problem [4] The conditions of the experiments and simulations are summarized in Table III and in the dimensionless groups in Eqs. ( 51)- (55).The parameters for the target experiments and the present simulation match exactly.In Table III, only the kinematic-viscosity ratio η in the previous simulation differs from the others.This was because and previous works with the lattice Boltzmann [84,85] and phase-field methods [86], where all simulations are performed with Re = 512.For the spike-tip evolution, no significant difference is observed among the lattice Boltzmann simulations.At later stages, these simulations are found to penetrate more deeply than the phase-field simulation.
筑波大学 University of Tsukuba 30 the previous D3Q19 MRT color-fluid model could not maintain the numerical stability when η = 4.2.Note that the present lattice Boltzmann model remains stable under the condition in Table III.

Boundary conditions
According to the grid refinement study in [33], we determine the inlet diameter D j0 = 25 and the computational domain 8D j0 × 8D j0 × 20D j0 .In this paper, the density of the dispersed phase ρ j (= ρ 0 r ) and the inlet velocity u j0 are set to be 1.0 and 0.1, respectively.Other parameters can be determined using the relations of dimensionless groups [Eqs.( 51)-( 55)] in lattice units: Figure 12 shows the interface evolution of the present simulation.Time is non-dimensionalized by the inlet diameter D j0 and the inlet velocity u j0 as in the liquid-jet-breakup simulations.At t * = 3, a mushroom-like head shape is created.At t * = 9 to 15, the interface of the jet becomes unstable.Observation shows that this interfacial instability is triggered by the return flow of a mushroom head generated in an initial stage of jet injection.Such a characteristic flow pattern leads to the onset of interfacial instability, although no artificial spatial or temporal perturbation has been assumed in the initial or boundary conditions.At later stages (t * = 30 and 72), the jet's leading-edge collapses, and the entrainment behavior, namely, the droplet formation from the side of the jet, is also observed.Finally, the liquid core becomes asymmetric, as can be seen at a The D3Q19 MRT color-fluid model was used.In the simulation, the kinematic-viscosity ratio η was set to be unity owing to the limitation of numerical stability.
t * = 72.Such a series of processes is similar to an experimental observation [33].
The history of jet penetration can be evaluated by both experiments and simulations.A comparison between experiments and simulations in terms of the position of the jet's leading-edge is provided in Fig. 13, where the position is normalized by the inlet diameter D j0 .Focusing on the results of lattice Boltzmann simulations, the time histories are almost the same until t * = 15.Afterward, the differences in positions gradually increase and the present simulation results tend to be close to the experimental data.This is thought to be due to the influence of ambient viscosity.In the present simulation, the lower viscosity of the surrounding fluid facilitates the penetration of the jet.The simulation using the developed lattice Boltzmann model shows a better agreement with the experimental data than that using the previous D3Q19 model [33].The enhanced numerical stability of the model enables us to simulate more realistic conditions.

C. Transition of breakup regimes
Recently, Saito et al. [22] have classified jet breakup regimes in liquid-liquid systems based on observations.This was an extension of Ohnesorge's for liquid-gas [7][8][9].The classified breakup regimes are follows: 0-dripping, I-varicose breakup, II-sinuous breakup, and III-atomization.Regime II was further divided into two sub-regimes: IIa-sinuous without entrainment, IIb-sinuous with entrainment.On the basis of the observations and phenomenological considerations, the flow-transition criteria were derived in [22]: for Regimes I and II, and for Regimes II and III, where the Reynolds number Re is defined as in Eq. ( 53) and the Ohnesorge number Oh can be given by Eqs. ( 53) and ( 54) as Oh = We 1/2 /Re.By using the above transition criteria, one can predict the breakup regimes of an immiscible liquid-liquid jet from initial parameters.With the present lattice Boltzmann model, we assess here the potential for the reproducibility of the breakup regimes.The simulation conditions are summarized in Table IV.The corresponding values on the dimensionless diagram are shown in Fig. 14.The description "Ref.case" represents the same condition as mentioned in Sec.IV B. As can be seen, the simulation results provided in the previous section are in Regime II, where sinuous breakup with or without entrainment is expected.The other parameters, namely, density ratio γ, viscosity ratio η, and Froude number Fr, are set to be the same as in Table III.In Case 1, Re is decreased by increasing the jet viscosity ν j and decreasing jet velocity u j0 while Oh is fixed at the same value as in the Ref. case.As can be seen in Fig. 14, Case 1 is in Regime 0 or I.It is expected that the dripping-or varicose-breakup regimes will appear.In Case 2, Re is fixed at the same value as in the Ref. case while Oh is increased by decreasing the interfacial tension σ.Since Case 2 is in Regime III in Fig. 14, it is expected that atomization will appear.Note that it is impossible to change only one parameter in the experimental procedures; only the effect of a focused parameter can be taken into account in the numerical simulations.
Figure 15 shows the simulation results for Case 1 in Table IV.We set the inlet diameter and velocity to D j0 = 15 and u j0 = 0.05, respectively.The other parameters are as follows: σ = 9.1 × 10 −3 , ν j = ν r = 1.6 × 10 −3 , ν c = ν b = 3.9 × 10 −4 , and g = 2.0 × 10 −5 .At t * = 17, the swollen part is generated around the inlet.The mushroom-like head does not appear.The swollen part moves downward with the growth of the neck part at t * = 26, and the corresponding part breaks up into a single droplet at t * = 30.At this time, the following swollen part is generated on a liquid column.The formation of the swollen part, the growth of the neck part, and breakup into a single droplet are observed through simulation.In this case, the so-called satellite-droplet formation just after the primary-droplet formation is not observed.According to the experimental data in liquidliquid systems [22], the average size of satellite droplet is about 0.3 times smaller than the nozzle diameter D j0 .To reproduce the satellite-droplet formation, the higher resolution would be required.The order of the droplet sizes has the same extent as the inlet diameter or is larger than it.This series of processes is a characteristic of the socalled pinch-off behavior, which is similar to the scenario considered in Rayleigh's breakup for a liquid column [2].This corresponds to the varicose breakup or Regime I in Fig. (14).The characteristics of the varicose breakup regime are reproduced by the present lattice Boltzmann simulation.
Figure 16 shows the simulation results for Case 2 in Table IV.We set the inlet diameter and velocity to D j0 = 30 and u j0 = 0.1, respectively.The other parameters are as follows: σ = 3.0 × 10 −5 , ν j = ν r = 8.8 × 10 −4 , ν c = ν b = 2.1×10 −4 , and g = 3.9×10 −5 .As in the Ref. case, a mushroom-like head appears at t * = 8.Through t * = 18-33, the jet continues penetration with active entrainment.The sizes of the generated droplets are much smaller than the inlet diameter, and the number of droplets is much   [33], and experimental result (square symbol) [33].The present simulation shows a better agreement with the experimental data than that using the previous D3Q19 model.
greater than in the Ref. case and Case 1. Owing to the entrainment from a jet interface, the liquid column is fully covered, and it is difficult to identify this column in the snapshots shown in Fig. 16.The series of processes corresponds to the atomization or Regime III in Fig. (14).The characteristics of the atomization regime are also reproduced by the present lattice Boltzmann simulations.
Before summarizing the simulation results, let us mention the jet-breakup regime shown in Fig. 12 again.During the penetration, the jet maintains an axisymmetric configuration.This is reasonable because the computa-筑波大学 University of Tsukuba tional setup described in Sec.IV A is exactly symmetric on the x-y plane.However, the jet column finally winds and becomes an asymmetric shape (see t * = 72 in Fig. 12).In addition, the entrainment also can be seen in this case.This type of breakup regime corresponds to Regime IIb in Fig. (14).This implies that the physical balance among hydrodynamic forces, including inertia, viscous force, and interfacial-tension force can be naturally reproduced via the lattice Boltzmann simulation.This fact numerically supports the validity of the dimen- sionless diagram proposed in Ref. [22].

V. CONCLUSIONS
In this paper, a three-dimensional lattice Boltzmann model for immiscible two-phase flows was developed in the framework of a D3Q27 lattice.An MRT collision operator for the D3Q27 lattice [63,68] was introduced.The choice of relaxation parameters optimized by Suga et al. [64] greatly improved the numerical stability of the model.A generalized perturbation operator [53] and an enhanced equilibrium distribution function [56] were also successfully incorporated into the present D3Q27 model.
Using the formulated lattice Boltzmann model, three types of numerical tests were carried out: a static droplet, an oscillating droplet, and the Rayleigh-Taylor instability.The static-droplet test shows that the measured pressure difference satisfied the Laplace law, and the interfacial tension agreed well with that predicted by Eq. ( 24) within a maximum error of 1.6%.The oscillating-droplet test was performed to compare the oscillation period with the analytical solution.For the high-interfacial-tension case, the error with the analytical solution was within 2.4%.Under low interfacial tension, a droplet reached equilibrium without oscillations in the most viscous case; the maximum error for the available data was 10.8%.The Rayleigh-Taylor simulations were performed with a computational setup and parameters that were strictly similar to those of He et al. [84].The positions of the bubble front, spike tip, and saddle point were measured for comparison with previous works using the lattice Boltzmann [84,85] and phase-field methods [86].The time changes of the bubble front and the saddle point were almost the same, irrespective of the method used.At later stages, a little difference between the lattice Boltzmann and phase-field simulations was observed for the spiketip evolution owing to the difference in implementation of the wall-boundary conditions.
The developed model was applied to liquid-jet-breakup simulations.First, we chose the parameters to match the experimental conditions [33].The five dimensionless groups were employed to determine the physical properties in lattice units.The developed D3Q27 model maintained numerical stability throughout the simulations; the previous work using a D3Q19 color-fluid lattice Boltzmann model [33] failed to simulate stably under the same conditions.The present results showed that the charac-teristic interfacial evolution captured the experimental results.A mushroom-like head was formed at the early stages; later, the configuration of the liquid core transitioned from asymmetric to axisymmetric, and interface entrainment also naturally occurred.The time history of the jet's leading-edge was compared with that obtained by experiment.Quantitative comparisons agreed well with the experimental data.
By choosing the parameters based on the regime map for jet breakup in liquid-liquid systems [22], we performed simulations to evaluate the reproducibility of the regime map.In the varicose regime, pinch-off-type breakup, i.e., the droplet formation in the tip of the liquid column, occurred.Satellite-drop formation was not confirmed in the present simulation.In the atomization regime, entrainment from the liquid column was distinguished.The liquid column was fully covered with entrainment droplets, which were much smaller than the inlet diameter.In conclusion, the breakup regimes appearing in the simulations successfully reproduced the predicted regimes, including the varicose, sinuous, and atomization regimes.
The authors also tested the stable ranges of the developed color-fluid model through simple test cases as in Sec.III A. For such a static case, at least, it is confirmed that the available maximum density ratio γ (= ρ r /ρ b ) was up to 1,000 at the unit kinematic-viscosity ratio (ν r = ν b = 1.0×10 −3 ); the maximum kinematic-viscosity ratio η (= ν r /ν b ) was up to 1,000 at γ = 1.5.From such a point of view, the authors believe that the threedimensional color-fluid lattice Boltzmann model developed in this paper is generally suitable for various other applications within the aforementioned stability ranges.
First, we focus on ϕ k i in the equilibrium function Eq. (11).The equilibrium function f k(e) i can be chosen arbitrarily to satisfy the conservations of mass and momentum: The form of target equation here is To satisfy Eqs.(B1) and (B2), the required condition for As in [50], we assume (B7) Next, we derive coefficients of the additional term Φ k i in Eq. (11).Let us consider the conservation of colorblinded variables Φ i = k Φ k i for simplicity.From the conservation of the 0-th to third order moments in the equilibrium functions, f , the relations to be satisfied are as follows [44,56]  Finally, we move on B i in Eq. (23).Following Ref. [53], the conditions for B i to be satisfied are to satisfy Eqs.(B17)-(B19).Note that the choice of coefficients in Eq. (B21) is somewhat arbitrary as in [53].

FIG. 2 .
FIG.2.Relationship between multiscale properties of a fluid flow and a simulation method.The lattice Boltzmann method is the so-called mesoscopic simulation method between the microscopic particle-based (e.g., molecular dynamics) and macroscopic Navier-Stokes-based methods.

FIG. 8 .
FIG. 8. Boundary conditions for the Rayleigh-Taylor instability.(a) The computational domain is discretized into an 128 × 128 × 512 lattice and the single-mode initial perturbation is initially imposed.(b) The top and bottom are no-slip boundaries and the lateral boundaries are periodic.

FIG. 9 .
FIG. 9. Interface evolution of the single-mode Rayleigh-Taylor instability with (a) Re = 512, (b) Re = 1, 024, and (c) Re = 5, 120.The Atwood number [Eq.(44)] is fixed at 0.5 in all simulations.At later stages, the higher the Reynolds number is, the more unstable the interface becomes.

FIG. 10 .
FIG. 10.Time history of the positions of the bubble front, spike tip, and saddle point.(a) Reynolds number dependence of the positions with Re = 512 (dash-dotted line), Re = 1, 024 (solid line), and Re = 5, 120 (dotted line).As pointed out by He et [84], when Re > Reynolds number dependence is negligible.(b) Comparison of the present study (solid line) and previous works with the lattice Boltzmann[84,85] and phase-field methods[86], where all simulations are performed with Re = 512.For the spike-tip evolution, no significant difference is observed among the lattice Boltzmann simulations.At later stages, these simulations are found to penetrate more deeply than the phase-field simulation.

FIG. 11 .
FIG. 11.Boundary conditions for a liquid-jet simulation.(a)The computational domain is discretized into a 8Dj0 ×8Dj0 × 20Dj0 lattice.The boundaries consist of an inflow boundary, a wall boundary, and an outflow boundary.(b) A circular inflow boundary with a uniform velocity uj0 is implemented at the top within Dj0.Free-slip[87] and convective boundary conditions[88] are implemented as the wall and outflow boundaries, respectively.

筑波大学 University of Tsukuba 10 t 72 FIG. 12 .
FIG. 12. Interface evolution of liquid-jet breakup under the same conditions as the target experiment: γ = 1.4,η = 4.2, Re = 3.4 × 10 3 , We = 2.2 × 10 2 , and Fr = 8.5.The computational domain is set to be 200 × 200 × 500.The interfacial instability is triggered by the return flow of a mushroom head generated in an initial stage of jet injection.The liquid core finally becomes asymmetric with entrainment behavior.Numerical stability is maintained throughout the simulation, although the kinematic viscosity for the continuous phase is set as low as 1.8 × 10 −4 .

FIG. 14 .
FIG.14.Location of simulation parameters on the regime map proposed in[22].The description "Ref."represents the same condition as mentioned in Sec.IV B. In Ref. case, sinuous-type regime is expected and will appear; in Case 1 and Case 2, it is expected that dripping/varicose and atomization, respectively, will appear.

Finally, we summarize
the simulation results in the dimensionless diagram of Fig. 14.Typical snapshots of liquid-jet breakup from Figs. 12, 15, and 16 are summarized in Fig 17.The characteristics of varicose breakup (Regime I), sinuous breakup (Regime II), and atomization (Regime III) are successfully simulated.On the breakup regimes, we can conclude that the lattice Boltzmann model developed in Sec.II under the boundary conditions described in Sec.IV A well reproduce the breakup characteristics expected by the regime map proposed in Ref. [22].

FIG. 17 .
FIG. 17. Summary of the present results of liquid-jet-breakup simulations in the dimensionless diagram for liquid-liquid systems[22].The lattice Boltzmann model developed in this paper well reproduce the breakup characteristics expected by the regime map.

TABLE I .
Parameters and evaluated errors of static droplet tests.

TABLE III .
Conditions of the experiments and simulations.

TABLE IV .
Simulation conditions for dimensionless numbers to be investigated.The density ratio, viscosity ratio, and Froude number are set to be the same as in TableIII