Lattice Boltzmann modeling and simulation of liquid jet breakup

TL;DR

This paper develops a 3D lattice Boltzmann model for simulating liquid jet breakup, accurately capturing experimental phenomena and regime behaviors in two-phase flows.

Contribution

It introduces a novel D3Q27 lattice Boltzmann model with enhanced stability and accuracy for simulating liquid jet breakup in two-phase flows.

Findings

Model accurately reproduces jet breakup regimes

Simulations match experimental conditions and results

Reproduces varicose, sinuous, and atomization regimes

Abstract

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.

Figures

Tables

Table 1. Table 1: Parameters and evaluated errors of static droplet tests.

$R$	$A$	$E$%
$30$	$1.0\times {10}^{-2}$	0.77
$25$	$1.0\times {10}^{-2}$	0.92
$20$	$1.0\times {10}^{-2}$	1.6
$30$	$4.0\times {10}^{-4}$	0.63
$25$	$4.0\times {10}^{-4}$	0.77
$20$	$4.0\times {10}^{-4}$	1.1

Table 2. Table 2: Parameters and evaluated errors of oscillating droplet tests.

$\nu$	$A$	$\sigma$	$E$%
$3.333\times {10}^{-3}$	$1.0\times {10}^{-2}$	$2.27\times {10}^{-3}$	0.60
$6.667\times {10}^{-3}$	$1.0\times {10}^{-2}$	$2.31\times {10}^{-3}$	0.83
$1.667\times {10}^{-2}$	$1.0\times {10}^{-2}$	$2.44\times {10}^{-3}$	2.4
$3.333\times {10}^{-3}$	$4.0\times {10}^{-4}$	$9.07\times {10}^{-5}$	2.7
$6.667\times {10}^{-3}$	$4.0\times {10}^{-4}$	$9.24\times {10}^{-5}$	10.8
$1.667\times {10}^{-2}$	$4.0\times {10}^{-4}$	$9.78\times {10}^{-5}$	–

Table 3. Table 3: Conditions of the experiments and simulations.

	Experiment and present simulation	Previous simulation111The D3Q19 MRT color-fluid model was used. In the simulation, the kinematic-viscosity ratio $\eta$ was set to be unity owing to the limitation of numerical stability.
$\gamma$	1.4	1.4
$\eta$	4.2	1.0
$\mathrm{Re}$	$3.4\times {10}^3$	$3.4\times {10}^3$
$\mathrm{We}$	$2.2\times {10}^2$	$2.2\times {10}^2$
$\mathrm{Fr}$	8.5	8.5

Table 4. Table 4: Simulation conditions for dimensionless numbers to be investigated. The density ratio, viscosity ratio, and Froude number are set to be the same as in Table  3

	$\mathrm{Re}$	$\mathrm{Oh}$	$\mathrm{We}={\left(\mathrm{Re}\cdot \mathrm{Oh}\right)}^2$
Ref. case	$3.4\times {10}^3$	$4.4\times {10}^{-3}$	$2.2\times {10}^2$
Case 1	$4.6\times {10}^2$	$4.4\times {10}^{-3}$	4.1
Case 2	$3.4\times {10}^3$	$3.0\times {10}^{-3}$	$1.0\times {10}^4$

Equations

$${\mathbf{c}}_i=(c_{ix},c_{iy},c_{iz})=\{\begin{array}{ll}(0,0,0)c,& i=1,\\ (\pm 1,0,0)c,(0,\pm 1,0)c,(0,0,\pm 1)c,& i=2,3,\dots ,7,\\ (\pm 1,\pm 1,0)c,(0,\pm 1,\pm 1)c,(\pm 1,0,\pm 1)c,& i=8,9,\dots ,19,\\ (\pm 1,\pm 1,\pm 1)c,& i=20,21,\dots ,27,\end{array}$$

$$f_i^k(\mathbf{x}+{\mathbf{c}}_i{\delta }_t,t+{\delta }_t)-f_i^k(\mathbf{x},t)={\Omega }_i^k(\mathbf{x},t),$$

$${\Omega }_i^k={\left({\Omega }_i^k\right)}^{(3)}\left[{\left({\Omega }_i^k\right)}^{(1)}+{\left({\Omega }_i^k\right)}^{(2)}\right],$$

$${\left(|{\Omega }^k⟩\right)}^{(1)}=-{\mathbf{M}}^{-1}\mathbf{K}\mathbf{M}\left(|f^k⟩-|f^{k(e)}⟩\right)+|F⟩.$$

$${\rho }_k=\sum _i{f}_i^k.$$

$$\rho \mathbf{u}=\sum _i\sum _k{f}_i^k{\mathbf{c}}_i+\frac{1}{2}\mathbf{F}{\delta }_t,$$

$${\bf M}=\left[\begin{array}[]{c}\bra{1}\\ \bra{c_{ix}}\\ \bra{c_{iy}}\\ \bra{c_{iz}}\\ -2\bra{1}+\bra{|{\bf c}_{i}|^{2}}\\ 2\bra{c_{ix}^{2}}-\bra{c_{iy}^{2}+c_{iz}^{2}}\\ \bra{c_{iy}^{2}-c_{iz}^{2}}\\ \bra{c_{ix}c_{iy}}\\ \bra{c_{iy}c_{iz}}\\ \bra{c_{iz}c_{ix}}\\ -4\bra{c_{ix}}+3\bra{c_{ix}c_{iy}^{2}+c_{iz}^{2}c_{ix}}\\ -4\bra{c_{iy}}+3\bra{c_{iy}c_{iz}^{2}+c_{ix}^{2}c_{iy}}\\ -4\bra{c_{iz}}+3\bra{c_{iz}c_{ix}^{2}+c_{iy}^{2}c_{iz}}\\ 4\bra{c_{ix}}-6\bra{c_{ix}c_{iy}^{2}+c_{iz}^{2}c_{ix}}+9\bra{c_{ix}c_{iy}^{2}c_{iz}^{2}}\\ 4\bra{c_{iy}}-6\bra{c_{iy}c_{iz}^{2}+c_{ix}^{2}c_{iy}}+9\bra{c_{ix}^{2}c_{iy}c_{iz}^{2}}\\ 4\bra{c_{iz}}-6\bra{c_{iz}c_{ix}^{2}+c_{iy}^{2}c_{iz}}+9\bra{c_{ix}^{2}c_{iy}^{2}c_{iz}}\\ 4\bra{1-|{\bf c}_{i}|^{2}}+3\bra{c_{ix}^{2}c_{iy}^{2}+c_{iy}^{2}c_{iz}^{2}+c_{iz}^{2}c_{ix}^{2}}\\ -8\bra{1}+12\bra{|{\bf c}_{i}|^{2}}-18\bra{c_{ix}^{2}c_{iy}^{2}+c_{iy}^{2}c_{iz}^{2}+c_{iz}^{2}c_{ix}^{2}}+27\bra{c_{ix}^{2}c_{iy}^{2}c_{iz}^{2}}\\ 2\bra{c_{iz}^{2}+c_{iy}^{2}}+3\bra{c_{ix}^{2}c_{iy}^{2}+c_{iz}^{2}c_{ix}^{2}}-4\bra{c_{ix}^{2}}-6\bra{c_{iy}^{2}c_{iz}^{2}}\\ 2\bra{c_{iz}^{2}-c_{iy}^{2}}+3\bra{c_{ix}^{2}c_{iy}^{2}-c_{iz}^{2}c_{ix}^{2}}\\ -2\bra{c_{ix}c_{iy}}+3\bra{c_{ix}c_{iy}c_{iz}^{2}}\\ -2\bra{c_{iy}c_{iz}}+3\bra{c_{ix}^{2}c_{iy}c_{iz}}\\ -2\bra{c_{iz}c_{ix}}+3\bra{c_{ix}c_{iy}^{2}c_{iz}}\\ \bra{c_{ix}c_{iy}^{2}-c_{iz}^{2}c_{ix}}\\ \bra{c_{iy}c_{iz}^{2}-c_{ix}^{2}c_{iy}}\\ \bra{c_{iz}c_{ix}^{2}-c_{iy}^{2}c_{iz}}\\ \bra{c_{ix}c_{iy}c_{iz}}\end{array}\right],$$

$$\begin{array}{rl}\mathbf{K}& =\mathrm{diag}\left[s_1,s_2,\dots ,s_{27}\right]\\ & =\mathrm{diag}\left[s_1,s_1,s_1,s_1,s_5,s_6,s_6,s_8,s_8,s_8,s_{11},s_{11},s_{11},s_{14},s_{14},s_{14},s_{17},s_{18},s_{19},s_{19},s_{21},s_{21},s_{21},s_{24},s_{24},s_{24},s_{27}\right],\end{array}$$

$$\nu$$

$$\zeta$$

$$f_i^{k(e)}({\rho }_k,\mathbf{u},{\alpha }_k)$$

$$={\rho }_k\left({\phi }_i^k+w_i\left[\frac{3}{c^2}({\mathbf{c}}_i\cdot \mathbf{u})+\frac{9}{2c^4}({\mathbf{c}}_i\cdot \mathbf{u}{)}^2-\frac{3}{2c^2}{\mathbf{u}}^2\right]\right)+{\Phi }_i^k.$$

$$w_i=\{\begin{array}{ll}8/27,& i=1,\\ 2/27,& i=2,3,\dots ,7,\\ 1/54,& i=8,9,\dots ,19,\\ 1/216,& i=20,21,\dots ,27.\end{array}$$

$${\phi }_i^k=\{\begin{array}{ll}{\alpha }_k,& i=1,\\ 2(1-{\alpha }_k)/19,& i=2,3,\dots ,7,\\ (1-{\alpha }_k)/38,& i=8,9,\dots ,19,\\ (1-{\alpha }_k)/152,& i=20,21,\dots ,27,\end{array}$$

$${\Phi }_i^k=\{\begin{array}{ll}-3\bar{\nu }(\mathbf{u}\cdot \nabla {\rho }_k)/c,& i=1,\\ +16\bar{\nu }({\mathbf{G}}_k:{\mathbf{c}}_i\otimes {\mathbf{c}}_i)/c^3,& i=2,3,\dots ,7,\\ +4\bar{\nu }({\mathbf{G}}_k:{\mathbf{c}}_i\otimes {\mathbf{c}}_i)/c^3,& i=8,9,\dots ,19,\\ +1\bar{\nu }({\mathbf{G}}_k:{\mathbf{c}}_i\otimes {\mathbf{c}}_i)/c^3,& i=20,21,\dots ,27,\end{array}$$

$$\bar{\nu }=\frac{1+\varphi }{2}{\nu }_r+\frac{1-\varphi }{2}{\nu }_b.$$

$$\varphi (\mathbf{x},t)=\frac{{\rho }_r(\mathbf{x},t)-{\rho }_b(\mathbf{x},t)}{{\rho }_r(\mathbf{x},t)+{\rho }_b(\mathbf{x},t)}.$$

$${\mathbf{G}}_k=\frac{1}{48}\left[\mathbf{u}\otimes \nabla {\rho }_k+(\mathbf{u}\otimes \nabla {\rho }_k{)}^{\mathrm{T}}\right].$$

$$\gamma =\frac{{\rho }_r^0}{{\rho }_b^0}=\frac{1-{\alpha }_b}{1-{\alpha }_r},$$

$$p_k={\rho }_k{\left(c_s^k\right)}^2={\rho }_k\frac{9(1-{\alpha }_k)}{19}{c}^2,$$

$$|F⟩={\mathbf{M}}^{-1}\left(\mathbf{I}-\frac{1}{2}\mathbf{K}\right)\mathbf{M}|F^{\prime }⟩,$$

$$|F^{\prime }⟩=w_i\left[3\frac{{\mathbf{c}}_i-\mathbf{u}}{c^2}+9\frac{({\mathbf{c}}_i\cdot \mathbf{u}){\mathbf{c}}_i}{c^4}\right]\cdot \mathbf{F}{\delta }_t.$$

$${\left({\Omega }_i^k\right)}^{(2)}=\frac{A_k}{2}|\nabla \varphi |\left[w_i\frac{({\mathbf{c}}_i\cdot \nabla \varphi )}{|\nabla \varphi {|}^2}-B_i\right].$$

$$B_i=\{\begin{array}{ll}-(10/27)c^2,& i=1,\\ +(2/27)c^2,& i=2,3,\dots ,7,\\ +(1/54)c^2,& i=8,9,\dots ,19,\\ +(1/216)c^2,& i=20,21,\dots ,27.\end{array}$$

$$\sigma =\frac{4}{9}A\tau c^4{\delta }_t,$$

$$({\Omega }_i^r{)}^{(3)}=\frac{{\rho }_r}{\rho }{f}_i+\beta \frac{{\rho }_r{\rho }_b}{{\rho }^2}\cos ({\theta }_i)f_i^{(e)}(\rho ,0,\bar{\alpha }),$$

$$({\Omega }_i^b{)}^{(3)}=\frac{{\rho }_b}{\rho }{f}_i-\beta \frac{{\rho }_r{\rho }_b}{{\rho }^2}\cos ({\theta }_i)f_i^{(e)}(\rho ,0,\bar{\alpha }),$$

$$\cos ({\theta }_i)=\frac{{\mathbf{c}}_i\cdot \nabla \varphi }{|{\mathbf{c}}_i||\nabla \varphi |},$$

$$\bar{\alpha }=\frac{1}{2}(1+\varphi ){\alpha }_r+\frac{1}{2}(1-\varphi ){\alpha }_b.$$

$$\frac{\partial \rho }{\partial t}$$

$$\frac{\partial (\rho \mathbf{u})}{\partial t}+\nabla \cdot (\rho \mathbf{u}\mathbf{u})$$

$$\mathbf{\Pi }=\rho \nu [\nabla \mathbf{u}+(\nabla \mathbf{u}{)}^{\mathrm{T}}]+\rho (\zeta -\nu )(\nabla \cdot \mathbf{u})\mathbf{I},$$

Full text

††preprint: Physical Review E

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 (Plateau, 1873) and Rayleigh (Rayleigh, 1878), extensive studies on this subject have been performed theoretically, experimentally, and numerically (McCarthy and Molloy, 1974; Lin and Reitz, 1998; Villermaux, 2007; Eggers and Villermaux, 2008).

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 (1936) classified his results into four types of breakup regimes: (0) dripping, (I) varicose, (II) sinuous, and (III) atomization (Kolev, 2005; McKinley and Renardy, 2011). 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.,  (Merrington and Richardson, 1947; Tanasawa and Toyoda, 1954; Grant and Middleman, 1966)). The majority of investigations have focused upon liquid-gas systems (liquid jet into gaseous atmosphere). Breakup of jets in liquid-liquid 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 (Meister and Scheele, 1969; Takahashi and Kitamura, 1971; Das, 1997) and CO2 storage in oceans (Riestenberg et al., 2004; Tsouris et al., 2007). 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 (Kondo et al., 1995; Dinh et al., 1999; Abe et al., 2006). In experiments, a high-temperature melt is often used to simulate the core melt materials. Abe et al. (2006) discussed the relationships between the fragment size and the Rayleigh–Taylor and Kelvin–Helmholtz interfacial instabilities (Chandrasekhar, 1961).

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 (Dinh et al., 1999; Saito et al., 2017). Saito et al. (2017) developed a breakup-regime map for a jet in immiscible liquid-liquid systems based on experiments and phenomenological considerations. Figure 1 shows the dimensionless map and the corresponding visual images. The Ohnesorge classification (Ohnesorge, 1936) was extended to liquid-liquid systems. As can be seen, various flow regimes occur during liquid-jet-breakup processes. The breakup transitions from the dropping regime to the atomization regime, and the generated droplet size drastically changes depending on the 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. (1993) investigated the axisymmetric steady-state laminar jet based on the volume-of-fluid (VOF) method (Hirt and Nichols, 1981). Thakre et al. (2015) 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 (Bürger et al., 1995); later, a similar variation was confirmed by experiments (Saito et al., 2017, 2015). Homma et al. (2006) numerically investigated liquid-liquid-jet breakup using a front-tracking method (Unverdi and Tryggvason, 1992; Tryggvason et al., 2001). They mapped different breakup modes on a plot of Weber number vs. viscosity ratio. The drawback of their front-tracking simulations was that they neglected the coalescence of generated drops. The aforementioned numerical simulations (Richards et al., 1993; Thakre et al., 2015; Homma et al., 2006) 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 (McCracken and Abraham, 2005a; Matsuo et al., 2015; Saito et al., 2016). In recent years, this method has been recognized as a powerful tool for analysis of complex fluid dynamics, including multicomponent and multiphase flows (Aidun and Clausen, 2010). 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 (Gunstensen et al., 1991; Grunau et al., 1993),

• •

Pseudo-potential model (Shan and Chen, 1993, 1994),

• •

Free-energy model (Swift et al., 1995, 1996),

• •

Mean-field model (He et al., 1999a).

This classification may not be exhaustive, for instance, the latter two models are sometimes identified as phase-field models (Li et al., 2016) 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 (Chen and Doolen, 1998; Nourgaliev et al., 2003; Aidun and Clausen, 2010; Chen et al., 2014; Liu et al., 2016; Li et al., 2016) 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 (Ba et al., 2016). 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 (Liu et al., 2016).

Color-fluid models, which are often referred to as R-K or color-gradient models, were first developed by Gunstensen et al. (1991), who extended the two-component lattice gas automata model of Rothman and Keller (1988). Later, Grunau et al. (1993) enabled the introduction of density and viscosity ratios by modifying the forms of the distribution functions. Latva-Kokko and Rothman (2005) 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-fluid-type model, namely, the lattice-pinning problem and the spurious velocities. Reis and Phillips (2007) 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. (2012) demonstrated that integrating Latva-Kokko–Rothman’s recoloring operator (Latva-Kokko and Rothman, 2005) into Reis–Phillips’ perturbation operator (Reis and Phillips, 2007) 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 (Leclaire et al., 2011).  Liu et al. (2012) derived a generalized perturbation operator using the phase-field (or order parameter) instead of a color-gradient and formulated the color-fluid model in three dimensions. Very recently,  Leclaire et al. (2017) 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. (1998), a source term to improve the Galilean invariance was derived by Leclaire et al. (2013a) 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. (2016) 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 (Luo et al., 2011). A multiple-relaxation-time (MRT) collision operator or generalized lattice Boltzmann equation (d’Humières, 1994; Lallemand and Luo, 2000; d’Humieres et al., 2002) 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 (2005b) successfully introduced an MRT operator to the multiphase lattice Boltzmann model and performed liquid-jet breakup simulations (McCracken and Abraham, 2005a). 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 liquid-jet 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 (Matsuo et al., 2015; Saito et al., 2016). Matsuo et al. (2015) used the three-dimensional two-phase lattice Boltzmann model, which was developed by Ebihara and Watanabe (2003) based on the model of He et al. (1999a). 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. (1999a) suffered from the dissolution of tiny droplets (Liu et al., 2016); thus, it would not be appropriate for liquid-jet-breakup simulations with tiny-droplet formation.

Saito et al. (2016) incorporated the MRT operator into the three-dimensional 19-velocity (D3Q19) color-fluid model proposed by Liu et al. (2012) 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.

In this paper, we present the three-dimensional two-phase 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 (Liu et al., 2012), and a formulatic-recoloring operator (Latva-Kokko and Rothman, 2005). A D3Q27 version of the enhanced equilibrium distribution functions (Leclaire et al., 2013a) 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-jet-breakup 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 (Saito et al., 2017). 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 (Geier et al., 2015) and relaxation parameters (Suga et al., 2015)

• •

Extending an enhanced equilibrium distribution function (Leclaire et al., 2013a) to the D3Q27 lattice

• •

Extending a generalized perturbation operator (Liu et al., 2012) to the D3Q27 lattice

For the present three-dimensional lattice Boltzmann model, distribution functions move on a D3Q27 lattice with the lattice velocity ${\bf c}_{i}$ defined as follows:

[TABLE]

where $c=\delta_{x}/\delta_{t}$, $\delta_{x}$ is the lattice spacing and $\delta_{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 (He and Luo, 1997).

In this model, two immiscible fluids are represented as pseudo red and blue fluids, respectively. The distribution function, $f_{i}^{k}$, 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 $f_{i}=f_{i}^{r}+f_{i}^{b}$. The evolution of the distribution function is expressed by the following lattice Boltzmann equation:

[TABLE]

where ${\bf x}$ and $t$ are the position and time, respectively. The collision operator $\Omega_{i}^{k}$ is made up of three sub-operators Tölke et al. (2002):

[TABLE]

where $(\Omega_{i}^{k})^{(1)}$ is the single-phase collision operator, $(\Omega_{i}^{k})^{(2)}$ is the perturbation operator, and $(\Omega_{i}^{k})^{(3)}$ is the recoloring operator. Using the MRT operator, the single-phase collision operator can be written as

[TABLE]

The density of the fluid $k$ is given by

[TABLE]

The total fluid density is given by $\rho=\sum_{k}\rho_{k}$, and the total momentum is defined as

[TABLE]

where ${\bf F}$ is the body force. Note that Eq. (6) indicates that the local velocity is modified to incorporate the spatially varying body force (Guo et al., 2002). In Eq. (2), ${\bf M}$ and ${\bf K}$ are, respectively, the $27\times 27$ transformation and relaxation matrices. In this paper, Dirac’s bra-ket notation is employed. Here, the “bra” operator $\bra{f}$ denotes a row vector along each lattice-velocity direction, i.e., $(f_{1},f_{2},\dots,f_{27})$, and the “ket” operator $\ket{f}$ denotes a column vector, i.e., $(f_{1},f_{2},\dots,f_{27})^{\rm T}$, where the superscript “${\rm 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 (2011) and Premnath and Banerjee (2011) presented the formulation independently in 2011. Dubois and Lallemand (2011) 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 (2011) 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. (2015) provided the following orthogonal moment set for the D3Q27 lattice:

[TABLE]

where $|{\bf c}_{i}|=\sqrt{c_{ix}^{2}+c_{iy}^{2}+c_{iz}^{2}}$. 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 ${\bf K}$ is the diagonal matrix given by (Suga et al., 2015)

[TABLE]

where the elements $0<s_{i}<2$ represent both the hydrodynamic and non-hydrodynamic relaxation parameters. The hydrodynamic parameters satisfy the following relations:

[TABLE]

where $\nu$ and $\zeta$ are the kinematic viscosity and the bulk viscosity, respectively. In this paper, we use the optimized parameters proposed by Suga et al. (2015): $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. (2013a) is used in this paper:

[TABLE]

In the case of $\Phi_{i}^{k}=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 (Leclaire et al., 2013a, 2014, 2015). 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 (2014).

The weights, $w_{i}$, are those of a standard D3Q27 lattice He and Luo (1997):

[TABLE]

Moreover, in the D3Q27 lattice, one can derive

[TABLE]

and

[TABLE]

where $\otimes$ is the tensor product and the symbol “$:$” indicates tensor contraction; $\bar{\nu}$ is the kinematic viscosity interpolated by (Liu et al., 2012; Tölke et al., 2002)

[TABLE]

Here, $\phi$ is the order parameter to distinguish the two components in a multicomponent flow, defined as (Liu et al., 2012)

[TABLE]

The values of the order parameter $\phi=1,-1$, and [math] correspond to a purely red fluid, a purely blue fluid, and the interface, respectively (Tölke et al., 2002). In the framework of D3Q27 lattice, the tensor ${\bf G}_{k}$ in Eq. (14) is defined as

[TABLE]

As established in (Grunau et al., 1993), the density ratio between the fluids, $\gamma$, must be taken into account as follows to obtain a stable interface:

[TABLE]

where the superscript “0” indicates the initial value of the density at the beginning of the simulation (Leclaire et al., 2013a). In each homogeneous phase region, the pressure of the fluid $k$ is given by

[TABLE]

for the D3Q27 lattice. This corresponds to an isothermal equation of state. In this paper, we choose $\alpha_{b}=8/27$, in which $c_{s}^{b}=1/\sqrt{3}$ (Leclaire et al., 2014; Saito et al., 2016).

The term $\Ket{F}$ in Eq. (4) represents the discrete forcing term accounting for the body force ${\bf F}$. In the MRT framework, the forcing term reads as (Yu and Fan, 2010)

[TABLE]

where ${\bf I}$ is a unit matrix, $\ket{F}=(F_{1},F_{2},\dots,F_{27})^{{\rm T}}$, and $\ket{F^{\prime}}=(F^{\prime}_{1},F^{\prime}_{2},\dots,F^{\prime}_{27})^{{\rm T}}$ is given by

[TABLE]

Eqs. (20) and (21) reduce to Guo et al.’s original forcing scheme (Guo et al., 2002) when using a single-relaxation time Yu and Fan (2010)

To model the interfacial tension, Liu et al. (2012) derived a generalized perturbation operator based on the CSF (Brackbill et al., 1992) concept, and the work of Reis and Phillips (2007) is employed to obtain the interfacial tension:

[TABLE]

Eq. (22) satisfies the correct form of the interfacial-tension 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:

[TABLE]

In this model, the interfacial tension can be directly given by

[TABLE]

where we assumed that $A=A_{r}=A_{b}$, $\tau$ is the relaxation time. Parameter $A$ controls the strength of interfacial tension, $\sigma$.

Although the perturbation operator, $(\Omega_{i}^{k})^{(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 Latva-Kokko and Rothman (2005); Halliday et al. (2007); Leclaire et al. (2012)

[TABLE]

where

[TABLE]

and $f_{i}^{(e)}$ is evaluated using Eq. (11), a zero velocity, and

[TABLE]

In the present model, the following continuity and Navier–Stokes equations can be derived via Chapman–Enskog analysis (Liu et al., 2012; Guo et al., 2002; Krüger et al., 2017)

[TABLE]

where

[TABLE]

is the viscous stress tensor with the shear viscosity $\nu$ given by Eq. (9) and the bulk viscosity $\zeta$ given by Eq. (10); $p=p_{r}+p_{b}$ is the pressure. In Eq. (30), the term $\nabla\cdot{\bf S}$ arises from the perturbation operator given by Eq. (22) and is equivalent to the interfacial force based on the CSF concept (Liu et al., 2012); the capillary stress tensor, ${\bf S}$, is given by

[TABLE]

For the computation of the gradient operator for an arbitrary function $\chi$, the following second-order isotropic central scheme (Liu et al., 2012; Guo et al., 2011; Liang et al., 2014; Lou et al., 2012)

[TABLE]

is adopted.

All the simulations in this paper are carried out in lattice units. In this paper, $\delta_{x}$ and $\delta_{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. (2013b).

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

III.1 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\times 85\times 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:

[TABLE]

where $W=4$ and $r=\sqrt{(x-x_{c})^{2}+(y-y_{c})^{2}+(z-z_{c})^{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\times 10^{-3}$, and gravity is neglected throughout the simulations. The parameters are $A$ in Eq. (24) and initial droplet radius $R$ (see Table 1). Periodic boundary conditions are imposed on all sides of the computational domain.

The Laplace equation in three dimensions is given by

[TABLE]

where $\Delta 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. (2012). 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 1 summarizes the simulation parameters and the evaluated error. The error $E$ is calculated as (Liu et al., 2012)

[TABLE]

where $\sigma_{th}$ and $\sigma_{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 (Saito et al., 2016). The maximum spurious velocities $|{\bf u}|_{max}$ are $1.2\times 10^{-2}$ for D3Q19 and $5.8\times 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.

III.2 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 (2007), except for the interfacial-tension coefficient. A droplet with a prolate-spheroid shape is placed at the center of a computational domain discretized by a $41\times 41\times 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 ($\nu=\nu_{r}=\nu_{b}$), and gravity is neglected throughout the simulations. The parameters are $\nu$ and $A$ (see Table 2). Parameter $A$ is the same as the ones provided in Sec. III.1.

The analytical solution of Miller and Scriven (1968) is used for comparison with the computed time periods. The frequency of the $n$-th mode is given by

[TABLE]

where $\omega_{n}$ is Lamb’s natural resonance frequency (Lamb, 1945)

[TABLE]

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 $\chi$ is given by

[TABLE]

We are only interested in the second mode $(n=2)$ here. The analytical expression for the time period is obtained by $T_{th}=2\pi/\omega_{2}$.

Figure 6 shows the transient change of the oscillating droplet shapes with $\nu=3.333\times 10^{-3}$ and $A=1.0\times 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. (Premnath and Abraham, 2007)

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 $\nu=1.667\times 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 2. The error $E$ is calculated as

[TABLE]

where $T_{sim}$ is the measured oscillation period. Note that, in the case of $\nu=1.667\times 10^{-2}$ and $A=4.0\times 10^{-4}$, we cannot measure the oscillation period since no oscillation occurred. The maximum errors are 2.4% for $A=1.0\times 10^{-2}$ and 10.8% for $A=4.0\times 10^{-4}$, respectively. Throughout, the accuracy is better when the kinematic viscosity is lower. It is found that the accuracy of the low-interfacial-tension case is difficult to assess via droplet-oscillation tests; however, the tests show that the present lattice Boltzmann model can be applied to unsteady two-phase-flow simulations.

III.3 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 (Chandrasekhar, 1961). The Rayleigh–Taylor instability has received considerable attention owing to its extensive applications, e.g., in the fundamental process of melt jet breakup (Abe et al., 2006; Saito et al., 2017). To our knowledge, this is the first time a color-fluid model has been applied to the three-dimensional Rayleigh–Taylor instability.

We refer to the computational setup adopted in He et al. (1999b). 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 (He et al., 1999b), the single-mode initial perturbation is initially imposed as

[TABLE]

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

[TABLE]

with ${\bf g}=(0,0,-g)$. The gravitational acceleration $g$ is chosen to satisfy $(Wg)^{1/2}=0.04$ (He et al., 1999a). The computational domain is set to be $W\times W\times 4W$ with $W=128$. The domain’s size resulted in a $128\times 128\times 512$ lattice.

The Atwood number

[TABLE]

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

[TABLE]

We use three patterns of Reynolds numbers: ${\rm Re}=512$, $1,024$, and $5,120$.

Figure 9 shows the interfacial evolution of the Rayleigh–Taylor instability. The dimensionless time is given by

[TABLE]

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 ${\rm Re}=5,120$ case [Fig. 9(c)], the interface becomes especially complicated. The Kelvin–Helmholtz instability appears conspicuously as the Reynolds number increases.

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. (1999b) pointed out that, when ${\rm 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. (1999b), three-dimensional Rayleigh–Taylor simulations were carried out using the lattice Boltzmann (Wang et al., 2016) and phase-field methods (Lee and Kim, 2013). Wang et al. (2016) used the phase-field-based MRT lattice Boltzmann model; Lee and Kim (2013) directly solved the Navier–Stokes–Cahn–Hilliard equations. Here, we compare our results with previous works (He et al., 1999b; Lee and Kim, 2013; Wang et al., 2016). Figure 10(b) shows the time histories of the positions at ${\rm Re}=1,024$ and ${\rm 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 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. The difference between the lattice Boltzmann and phase-field methods is considered to arise from the different wall-boundary conditions implemented.

IV Liquid jet breakup

IV.1 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. (2016), except for the outflow boundary. In the initial state, the computational domain is filled with blue-particle-distribution functions, $f_{i}^{b}$, 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 (Succi, 2001) is implemented as a wall boundary condition. At the outflow, a convective boundary condition (Lou et al., 2013) is used, in which the following convective equation for the distribution functions

[TABLE]

is solved, where $N$ is the node on the outflow boundary. Following Lou et al. (2013), we set two ghost-nodes $N+1$ and $N+2$. The discretized form can be given by the first-order implicit scheme,

[TABLE]

where $\lambda=U_{c}(t+\delta_{t})\delta_{t}/\delta_{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 (Orlanski, 1976). Through some numerical tests, we determine that the local velocity is suitable to the present system, that is,

[TABLE]

where $u_{z}({\bf x},t)=u_{z}(x,y,z,t)$ is the $z$-direction component of fluid velocity ${\bf u}$.

The body force in Eq. (21) for the liquid-jet simulations is set as

[TABLE]

with ${\bf g}=(0,0,g)$. Eq. (50) means that the gravitational force acts only in the dispersed phase (Chen et al., 2014).

IV.2 Comparison with experiments

Using the lattice Boltzmann model developed in Sec. II and the boundary conditions provided in Sec. IV.1, 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 (Saito et al., 2016), and the results are compared. In the target experiments, a glycerin-water-mixture jet was injected into a silicon-oil pool. These test fluids were immiscible with each other. In (Saito et al., 2016), 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 Lin and Reitz (1998)

[TABLE]

The conditions of the experiments and simulations are summarized in Table 3 and in the dimensionless groups in Eqs. (51)–(55). The parameters for the target experiments and the present simulation match exactly. In Table 3, only the kinematic-viscosity ratio $\eta$ in the previous simulation differs from the others. This was because the previous D3Q19 MRT color-fluid model could not maintain the numerical stability when $\eta=4.2$. Note that the present lattice Boltzmann model remains stable under the condition in Table 3.

According to the grid refinement study in (Saito et al., 2016), we determine the inlet diameter $D_{j0}=25$ and the computational domain $8D_{j0}\times 8D_{j0}\times 20D_{j0}$. In this paper, the density of the dispersed phase $\rho_{j}(=\rho_{r}^{0})$ 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: $\rho_{c}=\rho_{b}^{0}=0.70$, $\sigma=1.1\times 10^{-3}$, $\nu_{j}=\nu_{r}=7.4\times 10^{-4}$, $\nu_{c}=\nu_{b}=1.8\times 10^{-4}$, and $g=4.7\times 10^{-5}$.

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

[TABLE]

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 $t^{*}=72$. Such a series of processes is similar to an experimental observation (Saito et al., 2016).

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 (Saito et al., 2016). The enhanced numerical stability of the model enables us to simulate more realistic conditions.

IV.3 Transition of breakup regimes

Recently, Saito et al. (2017) have classified jet breakup regimes in liquid-liquid systems based on observations. This was an extension of Ohnesorge’s classification for liquid-gas systems Ohnesorge (1936); Kolev (2005); McKinley and Renardy (2011). The classified breakup regimes are as 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 Saito et al. (2017):

[TABLE]

for Regimes I and II, and

[TABLE]

for Regimes II and III, where the Reynolds number ${\rm Re}$ is defined as in Eq. (53) and the Ohnesorge number ${\rm Oh}$ can be given by Eqs. (53) and (54) as ${\rm Oh}={\rm We}^{1/2}/{\rm 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 4. 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.2. 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 $\gamma$, viscosity ratio $\eta$, and Froude number ${\rm Fr}$, are set to be the same as in Table 3. In Case 1, ${\rm Re}$ is decreased by increasing the jet viscosity $\nu_{j}$ and decreasing jet velocity $u_{j0}$ while ${\rm 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, ${\rm Re}$ is fixed at the same value as in the Ref. case while ${\rm Oh}$ is increased by decreasing the interfacial tension $\sigma$. 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 4. We set the inlet diameter and velocity to $D_{j0}=15$ and $u_{j0}=0.05$, respectively. The other parameters are as follows: $\sigma=9.1\times 10^{-3}$, $\nu_{j}=\nu_{r}=1.6\times 10^{-3}$, $\nu_{c}=\nu_{b}=3.9\times 10^{-4}$, and $g=2.0\times 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 liquid-liquid systems (Saito et al., 2017), 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 so-called pinch-off behavior, which is similar to the scenario considered in Rayleigh’s breakup for a liquid column (Rayleigh, 1878). 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 4. We set the inlet diameter and velocity to $D_{j0}=30$ and $u_{j0}=0.1$, respectively. The other parameters are as follows: $\sigma=3.0\times 10^{-5}$, $\nu_{j}=\nu_{r}=8.8\times 10^{-4}$, $\nu_{c}=\nu_{b}=2.1\times 10^{-4}$, and $g=3.9\times 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 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 computational setup described in Sec. IV.1 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 dimensionless diagram proposed in Ref. (Saito et al., 2017).

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.1 well reproduce the breakup characteristics expected by the regime map proposed in Ref. (Saito et al., 2017).

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 (Dubois and Lallemand, 2011; Geier et al., 2015) was introduced. The choice of relaxation parameters optimized by Suga et al. (2015) greatly improved the numerical stability of the model. A generalized perturbation operator (Liu et al., 2012) and an enhanced equilibrium distribution function (Leclaire et al., 2013a) 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. (1999b). The positions of the bubble front, spike tip, and saddle point were measured for comparison with previous works using the lattice Boltzmann (He et al., 1999b; Wang et al., 2016) and phase-field methods (Lee and Kim, 2013). 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 spike-tip 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 (Saito et al., 2016). 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 (Saito et al., 2016) failed to simulate stably under the same conditions. The present results showed that the characteristic 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 (Saito et al., 2017), 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.1. For such a static case, at least, it is confirmed that the available maximum density ratio $\gamma~{}(=\rho_{r}/\rho_{b})$ was up to 1,000 at the unit kinematic-viscosity ratio ($\nu_{r}=\nu_{b}=1.0\times 10^{-3}$); the maximum kinematic-viscosity ratio $\eta~{}(=\nu_{r}/\nu_{b})$ was up to 1,000 at $\gamma=1.5$. From such a point of view, the authors believe that the three-dimensional color-fluid lattice Boltzmann model developed in this paper is generally suitable for various other applications within the aforementioned stability ranges.

Acknowledgements.

This work was supported by Mitsubishi Heavy Industries, Ltd. The authors are grateful to H. Sakaba and H. Sato for their invaluable support. S.S. is grateful to A. Kaneko and Y. Iwasawa for the helpful discussions. The support of JSPS KAKENHI Grant Number 16J02077 is also acknowledged.

Appendix A Transformation matrix

A transformation matrix in practical form is obtained by substituting the velocity set [Eq. (1)] into the moment set [Eq. (7)]:

[TABLE]

and its inverse is analytically given by

[TABLE]

where $c=1$ is assumed.

Appendix B Derivation of lattice-specific coefficients

The lattice-specific coefficients for the D3Q27 lattice are $\varphi_{i}^{k}$ in Eq. (13), $\Phi_{i}^{k}$ in Eq. (14), and $B_{i}$ in Eq. (23). Derivations of these coefficients are provided in this appendix.

First, we focus on $\varphi_{i}^{k}$ in the equilibrium function Eq. (11). The equilibrium function $f_{i}^{k(e)}$ can be chosen arbitrarily to satisfy the conservations of mass and momentum:

[TABLE]

The form of target equation here is

[TABLE]

To satisfy Eqs. (61) and (62), the required condition for $\varphi_{i}^{k}$ is

[TABLE]

As in (Reis and Phillips, 2007), we assume

[TABLE]

and

[TABLE]

where $r$ is some constant. By choosing $r=4$ from the relations of weight factors $w_{i}$ (i.e., $w_{2-7}/w_{8-19}=w_{8-19}/w_{20-27}=4$), one can obtain the following form

[TABLE]

Next, we derive coefficients of the additional term $\Phi_{i}^{k}$ in Eq. (11). Let us consider the conservation of color-blinded variables $\Phi_{i}=\sum_{k}\Phi_{i}^{k}$ for simplicity. From the conservation of the 0-th to third order moments in the equilibrium functions, $f_{i}^{(e)}=\sum_{k}f_{i}^{k(e)}$, the relations to be satisfied are as follows (Nourgaliev et al., 2003; Leclaire et al., 2013a)

[TABLE]

To satisfy Eqs. (68)–(71), we use a form similar to (Leclaire et al., 2013a):

[TABLE]

with

[TABLE]

where $A_{i}$ are the lattice-specific coefficients to be determined and $C$ is the arbitrary constants. We assume

[TABLE]

with $r=4$. The resultant form is

[TABLE]

with

[TABLE]

Finally, we move on $B_{i}$ in Eq. (23). Following Ref. (Liu et al., 2012), the conditions for $B_{i}$ to be satisfied are

[TABLE]

In addition to the relations (77)–(79), we assume

[TABLE]

with $r=4$. One can choose the following coefficients

[TABLE]

to satisfy Eqs. (77)–(79). Note that the choice of coefficients in Eq. (81) is somewhat arbitrary as in (Liu et al., 2012).