Color-gradient lattice Boltzmann model with nonorthogonal central moments: Hydrodynamic melt-jet breakup simulations

We develop a lattice Boltzmann (LB) model for immiscible two-phase flow simulations with central moments (CMs). This successfully combines a three-dimensional nonorthogonal CM-based LB scheme [A. De Rosis, Phys. Rev. E 95, 013310 (2017)] with our previous color-gradient LB model [S. Saito, Y. Abe, and K. Koyama, Phys. Rev. E 96, 013317 (2017)]. Hydrodynamic melt-jet breakup simulations show that the proposed model is significantly more stable, even for flow with extremely high Reynolds numbers, up to $O(10^6)$. This enables us to investigate the phenomena expected under actual reactor conditions.


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, and understanding the breakup of liquid jets has been a topic of significant interest for more than a century. Since the pioneering works of Plateau [1] and Rayleigh [2], this subject has been extensively studied both theoretically, experimentally, and numerically [3][4][5][6]. In the linear theory framework, the liquid jet breakup problem is described in terms of the density ratio γ, viscosity ratio η, Reynolds number Re, Weber number We, and Froude number Fr as follows [4]: where ρ is the density, ν is the kinematic viscosity, u j0 is the jet velocity, D j0 is the jet inlet diameter, σ is the interfacial tension, and g is the acceleration due to gravity. The subscript j and c denote the dispersed and the continuous phases, respectively. At low injection velocities, drops form directly at the nozzle, while at higher velocities a liquid jet issues from the nozzle and then breaks into various droplet patterns. Discovering when these regimes occur is of significant interest in the study of liquid-jet breakup. Ohnesorge [7] classified his results into four breakup regimes: dripping (0), varicose (I), sinuous (II), and atomization (III) [8,9]. He also provided a map of these regimes for liquid jets in a gas in terms of the Reynolds number Re and the Ohnesorge number Oh, where Oh = We 1/2 /Re and can thus be calculated using Eqs. (3) and (4). Following Ohnesorge's work, there has been much research on this subject (see, e.g., [10][11][12]), most of which has focused on liquid-gas systems (liquid jets into gaseous atmospheres). Jet breakup in liquid-liquid systems (liquid jets into other liquids) has not been investigated as extensively.
Recently, Saito et al. [13] used a series of observations to classify the jet breakup regimes in liquid-liquid systems, as shown in Fig. 1(a), extending Ohnesorge's classification scheme for liquid-gas systems [7][8][9]. This classification largely follows Ohnesorge's one, but it further divides Regime II into two sub-regimes: sinuous without entrainment (IIa) and sinuous with entrainment (IIb). On the basis of observations and phenomenological considerations, they derived the following flow-transition criteria [13]: for Regimes I and II, and for Regimes II and III. These criteria can be used to predict the breakup regimes of immiscible liquid-liquid jets based on their initial parameters. We have successfully reproduced these flow-transition criteria in numerical simulations [14] based on their multiple-relaxation-time (MRT) lattice Boltzmann (LB) two-phase flow model, as illustrated in Fig. 1(b).  [13]. (b) Map of jet breakup regimes in immiscible liquid-liquid systems [13,14] Liquid-liquid-jet systems can be found in several fields, e.g., chemical processing [15][16][17] and CO 2 storage in oceans [18,19]. In the nuclear engineering field, it is important to fully understand the interactions between melt and coolant when designing nuclear reactor safety. As a result, the dispersion of liquid metal in water has been extensively investigated in the literature [20][21][22][23], going all the way back to G.I. Taylor's classic experiment with mercury and water [24]. In these experiments, hightemperature melt and water are often used to simulate the core melt materials and coolant. Using the experimental facility at University of Tsukuba (UT), Matsuo et al. [23,25] injected molten alloy with a melting point of 78 o C (U-Alloy78) into a water pool, using high-speed visualization to help understand the mechanism behind melt-jet breakup in water. The experiments of Magallon et al. [26], called FARO-TERMOS (FT), are unique in that they used a liquid sodium coolant: they poured pure molten UO 2 into a pool of liquid sodium. The fragment size analysis showed that fine fragments were generated by interactions between the molten UO 2 and the liquid sodium, but they obtained little physical insight into the breakup behavior because liquid sodium is not transparent. In this paper, our simulation targets are these two melt-coolant experiments [23,26]; we call them the UT and FT experiments, respectively.
The complexity of the phenomena involved in meltcoolant interactions means it is difficult to understand all the mechanisms simultaneously. Investigating the hydrodynamic interactions separately will thus help us to better understand the fundamental melt-jet breakup processes. Several researchers have already attempted to simulate jet breakup behavior using the volume-offluid method (see, e.g., [27][28][29]). In this paper, we use the LB method for multiphase flows, which has come to be recognized as a powerful tool for analyzing complex fluid dynamics, including multicomponent and multiphase flows [30]. The LB method is a mesoscopic simulation approach that lies between microscopic particlebased (e.g., molecular dynamics) and macroscopic Navier-Stokes-based methods.
CFD methods, which are based on the Navier-Stokes equations, the LB method, which uses mesoscopic kinetic equations, has several advantages, such as making it easy to incorporate mesoscale physics like interfacial breakup and coalescence. In addition, the computational cost of simulating realistic fluid flows is far more reasonable than with particle-based methods (e.g., molecular dynamics). Two-phase or multiphase LB models can be divided into four categories, namely, color-gradient [31,32], pseudopotential [33,34], free-energy [35,36], and meanfield [37] models. This is not an exhaustive classification; for instance, the latter two model types are sometimes called phase-field models [38] since the Cahn-Hilliard (or similar) interface tracking equations can be derived from them. For further details about multiphase LB models, interested readers can refer to several comprehensive review papers [30,[38][39][40][41][42] and references therein. This paper focus on color-gradient (CG) models, as they have many strengths for simulating multiphase or multicomponent flows, including strict mass conservation for each fluid and flexibility in adjusting the interfacial tension [43]. They also do not require us to use the static drop test to determine the interfacial tension, as this can be obtained directly without further analysis or assumptions. In addition, CG models exhibit very low dissolution for tiny droplets or bubbles [42].
CG models, often called R-K models, were first developed by Gunstensen et al. [31], who extended Rothman and Keller's two-component lattice gas automata model [44]. Later, Grunau et al. [32] enabled density and viscosity ratios to be introduced by modifying the forms of the distribution functions. Latva-Kokko and Rothman [45] then replaced Gunstensen's maximization-recoloring step with a formulaic segregation algorithm. Instead of widening the interface, Latva-Kokko-Rothman's recoloring algorithm solves two issues with the previous CG models, namely, the lattice-pinning problem and spurious velocities. Reis and Phillips [46] extended the model to a common two-dimensional nine-velocity (D2Q9) lattice and modified the perturbation operator to correctly recover the Navier-Stokes equations. Leclaire et al. [47] demonstrated that combining Latva-Kokko-Rothman's recoloring operator [45] with Reis-Phillips' perturbation operator [46] greatly improves the numerical stability and accuracy of the solutions over a wide parameter range. Using an isotropic gradient operator also enhances the numerical stability and accuracy [48]. Liu et al. [49] derived a generalized perturbation operator using the phase-field (or order parameter), and formulated the CG model in three dimensions. Leclaire et al. [50] generalized the CG model to two and three dimensions. For interested readers, Leclaire's MATLAB scripts will help to understand how to code the CG model [51] The so-called BGK [52] approximation refers to this simplest form of the collision operator, which forces all populations to relax towards an equilibrium state with the same rate. Despite its simplicity and phenomenal popularity, the BGK LB method is known to suffer from numerical instability under high-Re (low-viscosity) conditions. One way to overcome this issue is to modify the collision operator [53]. For example, MRT collision operators [54][55][56] have been widely used, even for multiphase flows, to enhance numerical stability and accuracy and reduce spurious currents near the interface. Later, Geier et al. [57] proposed a new collision operator based on the relaxation of central moments (CMs), that can be obtained by shifting the lattice directions according to the local fluid velocity. Many studies have developed this approach to fully exploit the properties of CM-based schemes (see, e.g., De Rosis [58] and references therein). For multiphase flows, Lycett-Brown and Luo [59] first introduced CMs into the pseudopotential multiphase LB model. Leclaire et al. [60] also introduced CMs into the CG model with unit density ratio. Very recently, De Rosis et al. [61] formulated a CM-based LB scheme for coupled Cahn-Hilliard-Navier-Stokes equations. De Rosis has consistently adopted nonorthogonal CMs [62][63][64][65][66], which are characterized by straightforward derivation and easy practical implementation. Moreover, his analytical formulation is very general, as it can be extended to any lattice velocity space.
In this paper, we present a three-dimensional CG LB model and apply it to hydrodynamic simulations of meltjet breakup. Sec. II describes the formulation of this LB model. Sec. III uses numerical tests on static droplets to evaluate the proposed model. Sec. IV applies the model to simulating melt-jet breakup under the conditions of the U T and FT experiments. Finally, Sec. V concludes this paper.

II. METHODOLOGY
The LB model presented here is based on our previous MRT CG model [14]. The most significant difference between the current and previous models is the introduction of De Rosis' nonorthogonal CMs [62,63]. In the current three-dimensional LB model, the distribution functions move on a three-dimensional 27-velocity (D3Q27) lattice [67]. We adopt a lattice speed c = δ x /δ t = 1, where x and t are the lattice spacing and time step, respectively. The lattice velocities c i = [|c ix , |c iy , |c iz ] are defined as follows: where i (= 0, 1, . . . , 26) represents the lattice-velocity directions and the superscript " " is the transpose operator. Here, we employ Dirac's bracket notation, where the "bra" operator ·| denotes a row vector along one of the lattice-velocity directions and the "ket" operator |· denotes a column vector. The model represents two immiscible fluids as red and blue fluids. Distribution functions f k i represent the fluids k, where k = r and b denote "red" and "blue," respectively, and i is the lattice-velocity direction. The total distribution function is defined as f i = f r i + f b i , and the evolution is expressed by the following LB equation: where x = [x, y, z] and t are the position and time, respectively. The collision operator Ω k i is made up of three sub-operators [68]: where (Ω k i ) (1) , (Ω k i ) (2) , and (Ω k i ) (3) are the single-phase collision, perturbation, and recoloring operators, respectively. As in Ref. [50], the single-phase and perturbation operators are applied using the color blind distribution function f i .
In this paper, we adopt the general MRT (GMRT) framework [69,70] to describe the single-phase collision operator with nonorthogonal CMs, due to the simplicity of its relationship to the MRT and SRT collision operators. It should be noted that Fei et al. [70] propose a simplified version of De Rosis' nonorthogonal CMs [62,63], showing a significantly reduced computational cost. In the GMRT framework, the single-phase collision operator can be written as where M, N, and K are the transformation, shift [69][70][71], and relaxation matrices, respectively. 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, in Eq. (13), the local velocity has been modified to incorporate the spatially varying body force [72]. To model the single-phase collision operator [Eq. (11)], we use the nonorthogonal CMs proposed by De Rosis [63], namely, where |c ix = |c ix − u x , |c iy = |c iy − u y , and |c iz = |c iz − u z . All the components of the vector |c i | 0 are equal to 1. The transformation matrix M, whose components are constant, transforms the distribution functions into raw moments. The shift matrix N, a lowertriangular matrix with components given by the macroscopic velocity u, transforms the raw moments into CMs and can be written as N = TM −1 . The practical forms of M, M −1 , N, and N −1 are given in Appendix A. The relaxation matrix K is a diagonal matrix given by where the elements are the moments' relaxation times. If s 0 = . . . = s i = . . . = s 6 , the model reduces to the BGK (single-relaxation-time) model. The subscripts represent the orders (e.g., 2 means second-order). The parameters s 2ν and s 2b satisfy the following relations: where ν and ζ are the kinematic and bulk viscosities, respectively. Here, we use s 0 = s 1 = 0 and s 2b = s 3 = s 4 = s 5 = s 6 = 1.
For the single-phase collision operator, we use the following enhanced equilibrium distribution function [73] in three-dimensions [14]: If Φ i = 0, Eq. (18) reduces to the standard form of an equilibrium distribution function up to third order. Using Eq. (18) improves the Galilean invariance of the variable density and viscosity ratios under the assumption of a small pressure gradient [73][74][75]. The weights, w i , are those of a standard D3Q27 lattice [76], as follows: In addition, for a D3Q27 lattice, we can derive and where ⊗ is the tensor product, ":" represents tensor contraction, andν is the kinematic viscosity, which interpolates between the red and blue viscosities ν r and ν b via the following harmonic mean [50,77,78]: Here, φ is the order parameter that distinguishes the two components in the multicomponent flow, defined as [43] where the superscript "0" indicates the initial density value at the beginning of the simulation [73]. The order parameter value φ = 1, −1, and 0 correspond to a purely red fluid, a purely blue fluid, and the interface between the two, respectively [68]. In the D3Q27 lattice framework, the tensor G in Eq. (21) is defined as As established in Ref. [32], in order to obtain a stable interface, we must take the fluid density ratio γ into account, which is defined as follows: The fluid pressures are given by an isothermal equation of state for the D3Q27 lattice: whereᾱ interpolates between α r and α b as follows [50] In this paper, we set The term |F in Eq. (11) is a discrete forcing term that accounts for the body force F. In the GMRT framework [69], it is where I is the unit matrix, |F = (F 0 , F 1 , . . . , F 26 ) , and |F = (F 0 , F 1 , . . . , F 26 ) is given by Equations (28) and (29)  To model the interfacial tension, we use the generalized perturbation operator derived in Ref. [49], based on the idea of continuum surface force (CSF) [81], and follow [46] to obtain the interfacial tension as follows: Equation (30) takes the correct form for an interfacial tension force in the Navier-Stokes equations when the lattice-specific variables B i are chosen correctly. We have derived the following B i values for the D3Q27 lattice framework: In this model, the interfacial tension can be given directly by where τ is the relaxation time and we have assumed that The parameter A controls the interfacial tension strength σ.
Although the perturbation operator (Ω k i ) (2) generates the interfacial tension, it does not guarantee the two fluids are immiscible. To promote phase segregation and maintain the interface, we apply the following recoloring operators [45,47,82] ( ( where θ i is the is the angle between ∇φ and c i , defined by Here, we set the parameter β to 0.7 to reproduce the correct interfacial behavior with as narrow an interface as possible [49,82,83]. For the current model, we can derive the following continuity and Navier-Stokes equations via Chapman-Enskog analysis [49,72,84] ∂ρ ∂t where is the viscous stress tensor, with D = 3 in the three dimensions; the shear viscosity ν is given by Eq. (16) and the bulk viscosity ζ is given by Eq. (17). In Eq. (37), the ∇ · S term arises from the perturbation operator given by Eq. (30) and, according to the CSF idea, is equivalent to the interfacial force [49]. The capillary stress tensor S is given by The solutions of the present model with CMs [63] are consistent with the Navier-Stokes equations to second order in diffusive scaling [59,62,85,86] with the body [72] and interfacial [49] forces.
To compute the gradient operator for an arbitrary function χ, we adopt the following second-order isotropic central scheme [49,[87][88][89]: In this paper, we set δ x and δ t to 1, as is usual in LB simulations. Although the above formulation focuses on two-component systems, it should also be straightforward to implement this model for systems with three or more components [90].

III. STATIC DROPLET TESTS
In this section, we carry out static droplet tests to evaluate whether the interfacial tension predicted by Eq. (32) is correct for various density ratios. We discretized the computational domain as a 100 × 100 × 100 lattice and immersed a static red droplet of radius R in a blue fluid. The initial density fields for each phase were as follows: where Here (x c , y c , z c ) is the center of the computational domain. We set the kinematic viscosity ratio to be 1, with each phase's kinematic viscosity being 1/6, and set the parameters A in Eq. (32) to 0.01 and the initial droplet radius R to 25. We neglect gravity throughout the simulations and imposed periodic boundary conditions on all sides of the computational domain.
The three-dimensional Laplace equation is given by where ∆p is the pressure difference across the droplet interface. We evaluated the pressure for each phase using Eq. (26) and measured it after 100 000 iterations via the procedure used by Leclaire et al. [47]. Table I summarizes the simulation parameters and resulting errors E, which were calculated as [49] where σ th and σ Lap are the interfacial tensions predicted by Eq. (32) and measured using the Laplace equation [Eq. (43)], respectively. Note that the droplets are always spherical at equilibrium, indicating that numerical stability was maintained for all density ratios. These tests confirm that the CM-based CG model described in Sec. II can predict the interfacial tension in static cases to within a maximum error of 0.40%. In addition, we measured the maximum spurious velocity |u max | in the domain at equilibrium, finding a maximum value of 1.22 × 10 −4 (Table I). This value is smaller than that of our MRT-based CG model [14] (|u max | = 5.8 × 10 −3 for γ = 1.5). These findings indicate that introducing the CMs into the single-phase collision operator can help to reduce the spurious velocity, contributing to enhance the numerical stability.

IV. JET BREAKUP SIMULATIONS
As mentioned in Sec. I, we will now simulate the following two experiments using the method presented in Sec. II: • UT (University of Tsukuba) experiments [23] (see Sec. IV B) • FT (FARO-TERMOS) experiments [26] (see Sec. IV C) In the following simulations, we neglect temperature changes and do not take phase-change effects (e.g., vaporization, condensation, or solidification) into account, meaning that these are strictly hydrodynamic simulations.
A. Setup Figure 3 illustrates the computational setup for our hydrodynamic melt-jet breakup simulations. The boundary conditions are the same as in Ref. [14]. Initially, the computational domain consists entirely of blue particle distribution functions f b i with zero velocity. The boundaries consist of an inflow boundary, wall boundaries, and an outflow boundary. There is a circular inflow boundary at the top of the domain, within (x − x c ) 2 + (y − y c ) 2 < (D j0 /2) 2 , where (x c , y c ) represents the center of the x-y plane. Here, the velocity u j0 is uniform, with corresponding equilibrium functions, and there are no artificial disturbances at this boundary. Wall boundaries cover the rest of the top and sides of the domain, with free-slip [67] boundary conditions. At the outflow boundary, we imposed a convective boundary condition [91], applying the following convective equation to the distribution functions: where N is the outflow boundary node. Following Lou et al. [91], we added two additional ghost nodes, N +1 and N + 2. The discretized form of the distribution functions can be given by the first-order implicit scheme where λ = U c (t + δ t )δ t /δ x . There are several possibilities for the convective velocity U c normal to the outflow boundary, such as the local, average, or maximum velocity [92]. After conducting some numerical tests, we determined that the local velocity was most suitable for the current system, namely, where We represented the body force in Eq. (29) as with g = (0, 0, g). This means that gravity only acts on the dispersed phase [41].

B. UT experiments
In the UT experiments [23], an alloy called U-Alloy78 (Osaka Asahi Co., Ltd.), with a melting point of 78 o C, was injected into a stagnant water pool under atmospheric pressure. We considered four cases with different nozzle diameters (D j0 = 7, 10, 15, and 20 mm) from Ref. [23] (see Table II). In these four cases, the melt and water temperatures were set to 270 o C and 70 o C, respectively, and the physical properties were as follows: ρ j = 8 183 kg/m 3 , ν j = 0.24 mm 2 /s, σ = 1.104 N/m, ρ c = 981 kg/m 3 , and ν c = 0.443 mm 2 /s. Note that they defined the jet velocity u j0 and inlet diameter D j0 as the contact velocity at the water surface and the nozzle diameter, respectively. In this case, we discretized the computational domain into an 8D j0 × 8D j0 × 40D j0 lattice with D j0 = 30, resulting in a total of 240 × 240 × 1 200 = 69 120 000 grid points being used in the simulation. We set the inlet velocity u j0 to 0.05, the jet density ρ j = ρ 0 r to 1, and the coolant density ρ c = ρ 0 b to 1/γ. We here mention the computational cost of the present simulations. For the simulations with 6 912 000 grids, it takes around 32 hours for 20 000 iterations with our computational environment. One way to reduce the cost may be using the reduced velocity model, such as D3Q15 or D3Q19 lattice models (see Refs. [63,70]).  simulation parameters used are shown (in lattice units) in the captions. All the simulations (Table II) were numerically stable, even for very low kinematic viscosities of O(10 −5 ). In all cases, the simulation results reproduce the qualitative interfacial behavior well, i.e., many fragments are generated as the jets penetrate each other, both around the leading edge and the side regions. However, there are still some differences in interfacial behavior between the experiments and the simulations, particularly for smaller D j0 values (e.g., at t = 0.15 s in Fig. 4 and t = 0.25 s in Fig. 5). One reason for this is the effect of gas entrapment: in the experiments, some of the gas was trapped in the water pool when the jet came in contact with the water surface, and the simulations did not take this into account. In contrast, the simulated shape of the interface agreed relatively well with the experiments for larger D j0 values (Figs. 6 and 7) because little gas was trapped in these cases. For Case 4, the details of the generated fragments and the flow field are shown in Fig. 8. We can find that the liquid jet column has large velocity, while the generated fragments has small velocity. In the snapshot at upstream region [ Fig. 8(b)], the fragments generate from the unstable liquid-jet interface. Most of the fragments in this region are stretched, which appear not to be spherical shapes. The velocity magnitude of stretched fragments are large, while that of spherical ones are small. In the snapshot at downstream region [ Fig. 8(c)], most of the fragments are spherical shapes with low velocity magnitude.
In order to make a quantitative comparison, we compared the evolution of the jet's leading edge over time for each condition, and the results are summarized in Fig. 9. This shows that the simulation results agree well with the experimental data [23]. In all cases, the experimental observations are ahead of the simulation early on, but the simulation passes the experiment in the later stages. This is due to differences in the inlet conditions between the experiments and simulations: in the experiments, a given mass of melt material is injected, so the injection velocity u j0 may decrease over time (i.e., is not constant), whereas the simulations assumed a constant injection velocity. This means that the simulations match the experiments particularly well in the early stages, during jet penetration.
Based on the evolution of the jet's leading edge over time, we can estimate the jet breakup length L, i.e., the length of the continuous liquid column emitted from the nozzle [3,4]. This is one of the metrics that characterize jet breakup behavior. Here, we estimated it via the procedure used by previous melt-jet experiments [23,93,94]. Table III compares the jet breakup lengths from Ref. [23] and with those in our simulations; the error E here is defined as where L exp and L sim are the breakup lengths obtained via the experiments and simulations, respectively. The  simulation's accuracy improves as D j0 increases, the error dropping from 31.4% for Case 1 to 10.8% for Case 3. Thus, these simulations were able to predict the experimental jet breakup length L to within a maximum error of 31.4%. As another quantitative comparison, we also evaluated the fragment diameters. The fragment diameter d was determined by the following procedures:

Convert the volume V into the equivalent spherical
diameter by d ≡ (6V /π) 1/3 Figure 10 shows histograms of the measured diameters; we have excluded the continuous liquid column from the nozzle from the calculation. From the figure, we can find that the smaller the nozzle diameter is, the higher the observation frequency of the droplet is. Except for D j0 = 7 mm [ Fig. 10(a)], all the distributions have a long tail to the right, similar to a log-normal distribution. Next, we compared the fragment diameters measured via the simulations with the experimental data and the predictions of hydrodynamic instability theories. We calculated the mass-median diameter as a fragment size metric, due to the shape of the distribution (Fig. 10). On the theoretical side, the first indicators we compared were the critical wavelengths of classical Rayleigh-Taylor (RT) and Kelvin-Helmholtz (KH) instabilities, λ cr,RT and λ cr,KH . Assuming a two-dimensional stratified geometry, these are [95] λ cr,RT = 2π σ (ρ j − ρ c )g 1 2 , Here, we used the jet velocity u j0 as the scaling velocity, assuming that the ambient fluid was stationary in Eq. (51). The second theoretical indicator we compared was the critical Weber number We cr , from which we obtained the critical droplet diameter d, as follows: The value of We cr depends on the assumptions made, so several values have been proposed, such as 12 (Pilch and Erdman [96]) and 18 (Matsuo et al. [25]). Previous studies [22,94,97] have pointed out that these theoretical quantities are related to the sizes of the fragments generated by jet breakup. Figure 11 compares the simulated and experimental fragment sizes. The aforementioned hydrodynamic instability theories [Eqs. (50)- (52)] are also shown in the graph. The error bars indicate the 95% confidence intervals for the median diameter, which we used to express the widths of the droplet size distributions even though they were not necessarily Gaussian [13]. As in the experiment [23], the simulation results are in agreement with both the KH [Eq. (51) (7), the jet's interfaces appear to be very unstable at the jet-side region both in the experiments and simulations. Since there is velocity difference between the jet and the ambient fluid, KH instability may be one of the reasons of the onset of interfacial instability. In addition, fragments generated at the tip of the jet (also at the jet side in some cases) are considered to be fragmented into smaller drops again due to the limitation of the critical Weber number. These simulations thus support the breakup mechanism proposed based on the experimental observations. Table IV compares the experimental and simulated median diameters; the error E is defined as where d m,exp and d m,sim are the experimental and simulated median diameters, respectively. The maximum error E was 41.9%, for Case 2.

C. FT experiments
The FT experiments were performed at the Joint Research Centre (JRC) in Ispra (Italy) by Magallon et al. [26]. Around 100-kg-scale of UO 2 melt was poured into a liquid-sodium pool. Two experiments, called T1 and T2, were carried out, with release diameters of 50 mm and 80 mm, respectively. In this paper, we focus on the T1 experiment. A notable feature of the FT experiments is that sodium in not transparent, so these simulations will hopefully help us to better understand the phenomena involved.
For these simulations, we discretized the computational domain into an 8D j0 × 8D j0 × 20D j0 lattice. Table V summarizes the simulation conditions, together with the corresponding dimensionless quantities, given by Eqs. (1)- (5). It should be noted that the Reynolds number in Table V (1.1 × 10 6 ) is extremely high from a multiphase LB perspective. The liquid-liquid-jet breakup flow-regime map [13] shown in Fig. 12 indicates we are in the atomization regime (Regime III) in this case. To examine the effect of grid resolution on jet breakup behavior, we considered two cases with different grids: a coarse case with D j0 = 30 and hence a 240 × 240 × 600  Figure 13 shows the simulation results for the coarse case. The left hand side [ Fig. 13(a)] shows the time evolution of the jet interface in real units, while the right hand side [ Fig. 13(b)] shows the calculated fragment size distribution at t = 0.125 s. The simulation parameters are given (in lattice units) in the caption. Even though the numerical conditions were somewhat extreme, numerical stability was maintained throughout the simulations These results show that using our CG LB model, based on nonorthogonal CMs, allowed Re to be increased significantly. In Fig. 13, most of the fragments are generated at the side of the jet due to entrainment, and they appear to be tiny compared with the inlet diameter D j0 . This is characteristic of the atomization regime (Regime III) [13]. This simulation therefore suggests that the jet state in the FT experiment should be similar to the atomization regime and, in fact, the debris collected in the T1 experiment was very fine (around 30-600 µm) [26]. That said, however, some unphysical points remain in terms of the fragment shapes and numerical dissolution. The resulting fragments were not spherical but irregular, and the jet's leading edge faced numerical dissolution at the later stage compared to that at the initial stages. In addition, Fig. 13(b) shows that the fragment-diameter distribution appears to be log-normal, but this pattern is interrupted near the minimum resolution, which we believe is due to the low spatial resolution.
To investigate the effect of spatial resolution, we then carried out a fine simulation, in which there were twice as many grid points in each direction. We used the same dimensionless quantities (Table V) as in the coarse case. Figure 14 shows the results of this simulation, with the simulation parameters again shown (in lattice units) in the caption. Compared with Fig. 13(a), the interfacial behavior seen in Fig. 14(a) is clearly different in terms of the fragment sizes and shapes: they are smaller and nearly spherical, and fragmentation now occurs around the jet's leading edge and side. The fragment-diameter distribution in Fig. 14(b) is also smoother than that in Fig. 13(b), resembling a log-normal distribution. Adopting a higher-resolution computational domain also im-  [26]. Here, we focus on their T1 experiment. The physical properties for the UO2 melt and liquid sodium were as follows [98]: ρj = 8 663 kg/m 3 , νj = 0.46 mm 2 /s, σ = 0.465 N/m, ρc = 856 kg/m 3 , and νc = 0.28 mm 2 /s. The dimensionless quantities, calculated using Eqs. (1)- (5), are also shown, including the density ratio γ, kinematic viscosity ratio η, Reynolds number Re, Weber number We, and Froude number Fr.  proved the numerical dissolution issue, as more physically reasonable results were obtained. In both cases (coarse and fine), the atomization regime appeared, implying that the grid resolutions used in this paper were sufficient to simulate the qualitative behavior of the entire jet. However, we can also observe that finer grids will be required to perform quantitative evaluations.

V. CONCLUSION
In this paper, we have extended the previous CG LB model [14] by introducing nonorthogonal CMs in three dimensions [63] within the GMRT framework [69,70]. Static droplet tests showed that our model can predict the interfacial tension for a range of density ratios (up to 1 000) to within a maximum error of 0.40% and also that it can greatly reduce the spurious velocity. We also applied our model to hydrodynamic melt-jet breakup simulations, targeting two different experiments, namely, the UT [23] and FT [26] experiments. Numerical simulations under corresponding conditions, including those equivalent to an actual reactor, demonstrated that our model was more stable. In particular, our model allowed the Reynolds number to be increased significantly, up to O(10 6 ). The UT simulations predicted the jet breakup length and median fragment-diameter to within maximum errors of 31.4% and 41.9%, respectively. The results of the FT simulation suggested that the jet in the experiment was in the atomization regime. To investigate the effect of grid resolution, the latter simulation was carried out at two grid resolutions, with a minimum spacing of ∆x = 0.83 mm. The results imply that finer grid resolutions will be required to evaluate the behavior accurately. Since the spatial resolution used was limited by our computing environment, we plan to use a higherperformance environment in our future work.  Re Oh Present UT simulations Previous simulations [14] FT experiment [26] FIG. 12. Flow regime region covered by the FT experiment together with the conditions our previous simulations ('o') [14]. Our current simulations ('×') involve substantially higher Re values compared with the previous ones.The present simulation condition is located extremely higher Re condition compared with our previous simulation. The hatched gray region shows the reactor conditions estimated in Ref. [13]; the FT experiment fell within this region. This diagram predicts that the breakup will be in the atomization regime (III). The conditions of UT simulations presented in Sec. IV B are also shown in the same map for reference. 099. The minimum spatial resolution in this case was ∆x = 1.67 mm. Although the computation is stable, some unphysical behavior can be seen: the fragments do not appear to be spherical, and the fragments that should appear around the jet's leading edge are not present. In addition, the fragment-diameter distribution appears to be log-normal but is interrupted near the minimum resolution ∆x.  Fig. 13(a), the interfacial behavior seen in Fig. 14(a) is clearly different in terms of the fragment sizes and shapes: the fragments are smaller and nearly spherical. In addition, the fragment-diameter distribution in Fig. 14(b) is smoother than that in Fig. 13(b), and the higher-resolution domain has also improved the numerical dissolution.