# Phase-Field-Based LBM Analysis of KHI and RTI in Wide Ranges of Density Ratio, Viscosity Ratio, and Reynolds Number

Authors: Xun Zhou; Bo Dong; Weizhong Li
Journal: International Journal of Aerospace Engineering
 (2020)

---

## Abstract

Numerous studies have elaborated the dominated roles of Kelvin-Helmholtz instability (KHI) and Rayleigh-Taylor instability (RTI) in the liquid sheet breakup and primary atomization. As for applications in aeronautics, the liquid-gas mixing generally occurs at the challenging conditions of a large density ratio and high Reynolds number. Hence, the evaluation of KHI and RTI under such challenging conditions will have great significance in better understanding the destabilizing mechanism of the liquid layer. To this end, a lattice Boltzmann multiple-relaxation-time (MRT) two-phase model, based on the conservative Allen-Cahn equation, is reconstructed for the present study. Preliminarily, the numerical stability and accuracy of this MRT model are tested by Laplace’s law under a large density ratio and high Reynolds number, including the sensitivity study to the values of mobility. Afterward, KHI and RTI are investigated in wide ranges of the Reynolds number, density ratio, and viscosity ratio. Numerical results indicate that the enhanced viscous force of light fluid with an increasing viscosity ratio notably suppresses the roll-ups of heavy fluid in KHI and RTI. As for the density ratio, it generally shows negative impacts on fluid-mixing in KHI and spike-spiraling in RTI. However, when the density ratio and the Reynolds number both arrive at high levels, the Kelvin-Helmholtz wavelets aroused by a dominated inertia force of heavy fluid trigger severe interface disintegration. The above results once more demonstrate the excellent ability of the present model in dealing with challenging conditions. Besides, the morphological characteristics of KHI and RTI at a high Reynolds number and large density ratio also greatly support the typical interface breakup mechanism observed in primary atomization.


---

## Body

## 1. Introduction

Among many fundamental and ubiquitous fluid phenomena in nature and engineering, the Kelvin-Helmholtz instability (KHI) and the Rayleigh-Taylor instability (RTI) have attracted much attention for their essential roles playing in the interface distortion and breakup [1]. To be specific, KHI occurs at a perturbed interface between two fluid flows with different tangential velocities, while RTI is aroused when a heavy fluid is accelerated by gravity against a light one [2, 3]. Recent studies have elaborated that these two instabilities dominate the liquid sheet breakup and primary atomization: KHI leading to the amplification of interface perturbations and formation of longitudinal waves, followed by RTI resulting in the formation of bulges on top of the wave crest and subsequent destabilization of the liquid sheet [4–6]. Generally, for applications in aeronautics, the liquid-gas mixing layer encountered in atomization owns a large density contrast, and the Reynolds number of gas happens to be much higher than that of liquid. In that perspective, it is necessary to conduct investigations of KHI and RTI by considering the effects of large density contrast and high Reynolds number simultaneously, which has great significance in better understanding the destabilizing mechanism of the liquid layer.In the past decades, extensive efforts have been devoted to probe KHI and RTI through theoretical analyses [7, 8], experiments [9, 10], and numerical simulations [11–15]. However, few studies deal with them under the conditions of large density contrast and high Reynolds number. In terms of KHI, Ceniceros and Roma [12] studied the long-time dynamics of KHI by a fully adaptive nonstiff method. They inspected the interface evolution at a high Reynolds number only for the cases of density- and viscosity-matched fluids. Rangel and Sirignano [14] investigated the effects of the density ratio and surface tension on the nonlinear growth of interface disturbance using a vortex-sheet discretization approach. They found out that the disturbance growth was suppressed by increasing the surface tension or the density ratio, and a bifurcation phenomenon was observed when density ratios were larger than 0.2. In a recent experimental study, Wan et al. [9] observed the single-mode KHI in a supersonic flow with the Atwood number being 0.81, and hydrodynamic simulations reported in their work can reproduce the resulting interface structure fairly well. Besides, a hydrodynamic model based on the evolution of KHI in shear viscous flows has been proposed in the work of Konovalov et al. [16] to specify the formation of nanostructures in materials subjected to the action of concentrated energy flows. As for RTI, the pioneer study was performed by Taylor [3] via theoretical analyses to explore the factors affecting the instability growth rate, followed by the experiment of Lewis [10] to test the theoretical conclusions. Then, in the work of Goncharov [17], the continuous bubble evolution of single-mode RTI at arbitrary Atwood numbers was provided by an analytical model from the earlier exponential growth to the nonlinear regime. Later, Wei and Livescu [13] investigated the growth of 2D single-mode RTI at low Atwood numbers by using direct numerical simulation. The mean quadratic growth was found at late times and sufficiently high Reynolds numbers.In the current work, we invoke the lattice Boltzmann method (LBM) to investigate KHI and RTI, given its distinct advantages in tracking interface evolution with large morphological deformations [18, 19]. To begin with KHI, Zhang et al. [20] used the two-phase LB model proposed by He et al. [21] to study this type of instability for density-matched incompressible fluids, focusing on the effects of surface tension at a low Reynolds number (Re=250). In order to assess KHI at a high Reynolds number, Fakhari and Lee [22] incorporated the multiple-relaxation-time (MRT) collision operator into a LB multiphase model. Their simulation, however, was limited to low density ratios with a Reynolds number up to 10000. Through an efficient discrete Boltzmann model, Gan et al. [23] investigated nonequilibrium and morphological characterizations of KHI in compressible flows, taking the effects of viscosity and heat conduction into account. Regarding the applications of LBM in RTI, most of the related studies account for incompressible fluids. For example, He et al. [21] proposed a LB multiphase model in the nearly incompressible limit, and it was applied to simulate RTI with a low density contrast at a moderate Reynolds number. Later, Liang et al. [24] studied RTI via a new phase-field-based MRT LB model in incompressible multiphase flow systems. The Reynolds number was increased to 30000 in their study.Based on the research listed above, the LBM has achieved great success in the fields of KHI and RTI. However, at present, it is still an open subject for the LBM with many challenges, especially those relevant to a large density contrast at a high Reynolds number. The existing multiphase LB models dealing with a large density contrast can be roughly classified into pseudo-potential-based type and phase-field-based type. For the pseudo-potential-based type, the models of Li et al. [25] and Xu et al. [26] are conducted at a density ratio in excess of 700. For the phase-field-based type, most of the models use the Cahn-Hilliard (CH) equation for interface tracking, such as those of Zu and He [19] and Yan and Zu [18]. Recently, Liang et al. [27] and Fakhari and Bolster [28] also apply the Allen-Cahn (AC) equation for interface tracking. Compared with the CH equation, the AC equation contains a lower-order diffusion term. Thus, according to the work of Chai et al. [29] and Wang et al. [30], the AC-based model theoretically has a higher numerical accuracy and stability in solving the index function and density field than the CH-based one. To this end, a lattice Boltzmann two-phase MRT model based on the conservative AC equation is reconstructed in this paper to enhance its numerical stability at a large density contrast and high Reynolds number. Then, it is adopted to investigate KHI and RTI by considering the effects of the viscosity ratio, density ratio, and Reynolds number. Special attention is paid on morphological characteristics under the challenging conditions of a high Reynolds number and large density ratio, which are rarely covered in the literature. The rest of this paper is therefore organized as follows: In Section 2, the present MRT model is introduced in detail. In Section 3, after the preliminarily evaluation, the present model is utilized to inspect KHI and RTI in wide ranges of the viscosity ratio, density ratio, and Reynolds number. In Section 4, some conclusions are summarized. It is worth noting that the postprocessing tool of Tecplot 360 is utilized to visualize all the numerical results in this paper.


## 2. Mathematical Method

### 2.1. The MRT LB Model for the Conservative Allen-Cahn Equation

With the MRT collision operator [29], the generalized evolution equation for the conservative AC equation [29, 31] can be written as
(1)fix+ciδt,t+δt=fix,t−Λijffj−fjeqx,t+δtFi−0.5ΛijfFjx,t,where fix,t and fieqx,t are the particle distribution function and its corresponding equilibrium distribution function for the order parameter ϕ at position x and time t, respectively. Here, the order parameter ϕ taking 1 and 0 is adopted to distinguish different fluids with the interface marked by the contour level of ϕ=0.5. ci is the discrete velocity of the ith direction, δt is the time step. Λijf stands for an element of the collision matrix Λf=M−1SfM, in which M is an orthogonal transformation matrix and Sf is a diagonal matrix given by (for the D2Q9 model)
(2)M=111111111−4−1−1−1−122224−2−2−2−21111010−101−1−110−20201−1−110010−111−1−100−20211−1−101−11−10000000001−11−1,Sf=diags0f,s1f,s2f,s3f,s4f,s5f,s6f,s7f,s8f.In Equation (1), Fi represents the forcing term in the velocity space [27, 32], and it is defined as
(3)Fi=ωici·∂tϕu+cs2λncs2,where ωi is the weight coefficient with the value given by ω0=4/9, ω1,⋯4=1/9, and ω5,⋯8=1/36, cs=c/3 is the sound speed, and u is the macroscopic velocity. n=∇ϕ/∣∇ϕ∣ denotes the unit vector normal to the interface. λ=4ϕ1−ϕ/W is a function of ϕ, and W is the interface thickness. The time derivative term ∂tϕu is introduced to eliminate the artificial term in the recovered equation.Through the orthogonal transformation matrixM, the right-hand side of Equation (1) can be projected onto the moment space via mf=Mf and mfeq=Mfeq [33, 34] as
(4)mf∗=mf−Sfmf−mfeq+δtI−Sf2F^,where f=f0,⋯f8T, feq=f0eq,⋯f8eqT, I is the unit tensor, and F^=MF=MF0,⋯,F8T is the forcing term in the moment space. Based on feq and F in the velocity space [27, 32], mfeq and F^ can be easily derived for the D2Q9 model as
(5)mfeq=ϕ1,−2,1,ux,−ux,uy,−uy,0,0T,(6)F^=000∂tϕux+cs2λnx−∂tϕux−cs2λnx∂tϕuy+cs2λny−∂tϕuy−cs2λny00,where ux and uy are the components of the macroscopic velocity u, nx and ny are the components of the unit vector n. Afterward, the streaming process is given as
(7)fix+ciδt,t+δt=fi∗x,t,where the postcollision distribution function fi∗ can be obtained via f∗=M−1mf∗.Applying the Chapman-Enskog expansion to Equation (1), the conservative AC equation [30, 35] can be recovered correctly as
(8)∂ϕ∂t+∇·ϕu=∇·M∇ϕ−λn,where the mobility M is determined by
(9)M=cs2τf−0.5δt,in which τf is the dimensionless relaxation time of fj, and it is related to the diagonal matrix Sf with s3f=s5f=1/τf. In the present model, the order parameter is calculated as
(10)ϕ=∑ifi,and the fluid density ρ should be consistent with the order parameter by taking the linear interpolation as
(11)ρ=ϕρA−ρB+ρB,where ρA and ρB stand for the densities of two different fluids corresponding to the order parameter ϕA=1 and ϕB=0, respectively.


### 2.2. The MRT LB Model for the Navier-Stokes Equations

In order to improve the numerical stability, the lattice Boltzmann equation for the incompressible NS equations is combined with the MRT collision operator as(12)gix+ciδt,t+δt=gix,t−Λijggj−gjeqx,t+δtRi−0.5ΛijgRjx,t,where gi and gieq are the density distribution function and its corresponding equilibrium distribution function, respectively. Ri is the forcing term in the velocity space. Λijg is an element of the collision matrix Λg=M−1SgM, in which Sg is a diagonal matrix given by (for the D2Q9 model)
(13)Sg=diags0g,s1g,s2g,s3g,s4g,s5g,s6g,s7g,s8g.Similar to Equation (3), the right-hand side of Equation (12) in the moment space reads
(14)mg∗=mg−Sgmg−mgeq+δtI−Sg2R^where mg=Mg, mgeq=Mgeq, and R^=MR are the corresponding matrices in the moment space with g=g0,⋯g8T, geq=g0eq,⋯g8eqT, and R=R0,⋯,R8T. Inspired by the BGK model of Liang et al. [27], the matrices mgeq and R^ for the D2Q9 model are evaluated as
(15)mgeq=06p+3ρu2−9p−3ρu2ρux−ρuxρuy−ρuyρux2−uy2ρuxuy,R^=ux∂xρ+uy∂yρ0−ux∂xρ−uy∂yρGx−GxGy−Gy2ux∂xρ−uy∂yρ/3ux∂yρ+uy∂xρ/3,where p is the hydrodynamic pressure. Gx and Gy are the components of the total force G=Fs+Fe+Fa. Here, Fs=−σαW∣∇ϕ∣∇ϕ∇·n is the surface tension force with the formulation recommended by Kim [36] and Ren et al. [37], where σ is the surface tension coefficient. Fe is the possible body force. At the interface, phase-field-based LB models fail to satisfy the continuity equation [38] and, hence, an additional interfacial force Fa=qau is introduced with the term qa determined by
(16)qa=MρA−ρB∇2ϕ−∇·λn.Ultimately, withg∗=M−1mg∗, the streaming process is achieved as
(17)gix+ciδt,t+δt=gi∗x,t.Based on the density distribution functiongi, the macroscopic quantities u and p can be evaluated as
(18)u=1ρ−0.5qa∑icigi+δt2Fs+Fe,p=cs21−ω0∑i≠0gi+δt2u·∇ρ−ω0ρu22cs2.Besides, through the Chapman-Enskog analysis, the kinematic viscosityν is determined by
(19)ν=cs2τg−0.5δt,where τg=1/s7g=1/s8g. Note that the value of viscosity is usually not uniform in a two-phase flow system. To make it smoothly across the interface, the popular treatments are supposed that the viscosity is a linear or inverse linear function of the order parameter [18, 21, 39]. In this work, the linear form is adopted,
(20)ν=ϕνA−νB+νA,where νA and νB are the kinematic viscosities of two different fluids corresponding to the order parameter ϕA=1 and ϕB=0, respectively. Hence, the elements s7g and s8g of the diagonal matrix Sg are not uniform in the whole system due to its relation with the kinematic viscosity. Regarding the derivative terms in the present model, the explicit Euler scheme is adopted to calculate the temporal derivative in Equation (6) [40], while the gradient and the Laplacian operator are determined by second-order isotropic central schemes used in the studies of Yan and Zu and Zu and He [18, 19].


## 2.1. The MRT LB Model for the Conservative Allen-Cahn Equation

With the MRT collision operator [29], the generalized evolution equation for the conservative AC equation [29, 31] can be written as
(1)fix+ciδt,t+δt=fix,t−Λijffj−fjeqx,t+δtFi−0.5ΛijfFjx,t,where fix,t and fieqx,t are the particle distribution function and its corresponding equilibrium distribution function for the order parameter ϕ at position x and time t, respectively. Here, the order parameter ϕ taking 1 and 0 is adopted to distinguish different fluids with the interface marked by the contour level of ϕ=0.5. ci is the discrete velocity of the ith direction, δt is the time step. Λijf stands for an element of the collision matrix Λf=M−1SfM, in which M is an orthogonal transformation matrix and Sf is a diagonal matrix given by (for the D2Q9 model)
(2)M=111111111−4−1−1−1−122224−2−2−2−21111010−101−1−110−20201−1−110010−111−1−100−20211−1−101−11−10000000001−11−1,Sf=diags0f,s1f,s2f,s3f,s4f,s5f,s6f,s7f,s8f.In Equation (1), Fi represents the forcing term in the velocity space [27, 32], and it is defined as
(3)Fi=ωici·∂tϕu+cs2λncs2,where ωi is the weight coefficient with the value given by ω0=4/9, ω1,⋯4=1/9, and ω5,⋯8=1/36, cs=c/3 is the sound speed, and u is the macroscopic velocity. n=∇ϕ/∣∇ϕ∣ denotes the unit vector normal to the interface. λ=4ϕ1−ϕ/W is a function of ϕ, and W is the interface thickness. The time derivative term ∂tϕu is introduced to eliminate the artificial term in the recovered equation.Through the orthogonal transformation matrixM, the right-hand side of Equation (1) can be projected onto the moment space via mf=Mf and mfeq=Mfeq [33, 34] as
(4)mf∗=mf−Sfmf−mfeq+δtI−Sf2F^,where f=f0,⋯f8T, feq=f0eq,⋯f8eqT, I is the unit tensor, and F^=MF=MF0,⋯,F8T is the forcing term in the moment space. Based on feq and F in the velocity space [27, 32], mfeq and F^ can be easily derived for the D2Q9 model as
(5)mfeq=ϕ1,−2,1,ux,−ux,uy,−uy,0,0T,(6)F^=000∂tϕux+cs2λnx−∂tϕux−cs2λnx∂tϕuy+cs2λny−∂tϕuy−cs2λny00,where ux and uy are the components of the macroscopic velocity u, nx and ny are the components of the unit vector n. Afterward, the streaming process is given as
(7)fix+ciδt,t+δt=fi∗x,t,where the postcollision distribution function fi∗ can be obtained via f∗=M−1mf∗.Applying the Chapman-Enskog expansion to Equation (1), the conservative AC equation [30, 35] can be recovered correctly as
(8)∂ϕ∂t+∇·ϕu=∇·M∇ϕ−λn,where the mobility M is determined by
(9)M=cs2τf−0.5δt,in which τf is the dimensionless relaxation time of fj, and it is related to the diagonal matrix Sf with s3f=s5f=1/τf. In the present model, the order parameter is calculated as
(10)ϕ=∑ifi,and the fluid density ρ should be consistent with the order parameter by taking the linear interpolation as
(11)ρ=ϕρA−ρB+ρB,where ρA and ρB stand for the densities of two different fluids corresponding to the order parameter ϕA=1 and ϕB=0, respectively.


## 2.2. The MRT LB Model for the Navier-Stokes Equations

In order to improve the numerical stability, the lattice Boltzmann equation for the incompressible NS equations is combined with the MRT collision operator as(12)gix+ciδt,t+δt=gix,t−Λijggj−gjeqx,t+δtRi−0.5ΛijgRjx,t,where gi and gieq are the density distribution function and its corresponding equilibrium distribution function, respectively. Ri is the forcing term in the velocity space. Λijg is an element of the collision matrix Λg=M−1SgM, in which Sg is a diagonal matrix given by (for the D2Q9 model)
(13)Sg=diags0g,s1g,s2g,s3g,s4g,s5g,s6g,s7g,s8g.Similar to Equation (3), the right-hand side of Equation (12) in the moment space reads
(14)mg∗=mg−Sgmg−mgeq+δtI−Sg2R^where mg=Mg, mgeq=Mgeq, and R^=MR are the corresponding matrices in the moment space with g=g0,⋯g8T, geq=g0eq,⋯g8eqT, and R=R0,⋯,R8T. Inspired by the BGK model of Liang et al. [27], the matrices mgeq and R^ for the D2Q9 model are evaluated as
(15)mgeq=06p+3ρu2−9p−3ρu2ρux−ρuxρuy−ρuyρux2−uy2ρuxuy,R^=ux∂xρ+uy∂yρ0−ux∂xρ−uy∂yρGx−GxGy−Gy2ux∂xρ−uy∂yρ/3ux∂yρ+uy∂xρ/3,where p is the hydrodynamic pressure. Gx and Gy are the components of the total force G=Fs+Fe+Fa. Here, Fs=−σαW∣∇ϕ∣∇ϕ∇·n is the surface tension force with the formulation recommended by Kim [36] and Ren et al. [37], where σ is the surface tension coefficient. Fe is the possible body force. At the interface, phase-field-based LB models fail to satisfy the continuity equation [38] and, hence, an additional interfacial force Fa=qau is introduced with the term qa determined by
(16)qa=MρA−ρB∇2ϕ−∇·λn.Ultimately, withg∗=M−1mg∗, the streaming process is achieved as
(17)gix+ciδt,t+δt=gi∗x,t.Based on the density distribution functiongi, the macroscopic quantities u and p can be evaluated as
(18)u=1ρ−0.5qa∑icigi+δt2Fs+Fe,p=cs21−ω0∑i≠0gi+δt2u·∇ρ−ω0ρu22cs2.Besides, through the Chapman-Enskog analysis, the kinematic viscosityν is determined by
(19)ν=cs2τg−0.5δt,where τg=1/s7g=1/s8g. Note that the value of viscosity is usually not uniform in a two-phase flow system. To make it smoothly across the interface, the popular treatments are supposed that the viscosity is a linear or inverse linear function of the order parameter [18, 21, 39]. In this work, the linear form is adopted,
(20)ν=ϕνA−νB+νA,where νA and νB are the kinematic viscosities of two different fluids corresponding to the order parameter ϕA=1 and ϕB=0, respectively. Hence, the elements s7g and s8g of the diagonal matrix Sg are not uniform in the whole system due to its relation with the kinematic viscosity. Regarding the derivative terms in the present model, the explicit Euler scheme is adopted to calculate the temporal derivative in Equation (6) [40], while the gradient and the Laplacian operator are determined by second-order isotropic central schemes used in the studies of Yan and Zu and Zu and He [18, 19].


## 3. Results and Discussion

### 3.1. Static Droplet

Before the simulation of KHI and RTI, the present model is preliminarily tested by Laplace’s law as well as the sensitivity study to the values of mobility. Initially, a static droplet with the radiusR=25 surrounded by the gas phase is placed at the center of a periodic domain with 150×150 grid points. To mimic the challenging condition of the large density ratio and high Reynolds number, the density ratio and kinematic viscosity of liquid and gas are set as ρA/ρB=1000 and νA=νB=0.002, respectively. According to Laplace’s law, the theoretical pressure difference between liquid and gas is Δp=σ/R at steady-state. Figure 1(a) compares the theoretical and numerical values of pressure difference for various surface tension coefficients, which are in good agreement. In the simulation, we fix s0f=s7f=s8f=1.0, s1f=s2f, and s4f=s6f as usual. It is found out that 0.75≤s1f=s2f≤1.2 and 0.6≤s4f=s6f≤1.3 all can give satisfying results. For simplicity, s1f=s2f=s4f=s6f=1.1 is set. Regarding the relaxation elements in Sg, similar values are adopted as those in Ref. [22], except s0g=s3g=s5g=1.0 and s7g=s8g related to kinematic viscosity.Figure 1
Validation and sensitivity study. (a) Comparison of the pressure difference between theoretical and numerical results. (b) Maximum kinetic energy for various mobility.
(a)(b)For the importance of mobility in the present model, a sensitivity study to its values is necessary. Figure1(b) depicts the maximum kinetic energy Emax in the system for various mobilities. Here, the kinetic energy and the dimensionless time are defined as E=ρu2/2 and t∗=tσ/ρBνBR, respectively. It is noteworthy that the present MRT model can successfully damp high-frequency oscillations of the kinetic energy at the early stages. Besides, the undervalue of mobility allows less damping time and lower kinetic energy, which also means smaller spurious currents (10−9) in the system. Taking the numerical stability into account as well, the mobility is selected as a reasonable value of M=0.01 in the following study.


### 3.2. Kelvin-Helmholtz Instability

In this section, KHI in a two-phase incompressible flow system is considered with the schematic diagram depicted in Figure2. The shear flow of two immiscible fluids with densities ρA and ρB takes place in a 2D channel of size (λ×2λ), whose top and bottom walls move in opposite directions. Four dimensionless groups are chosen to regulate the present issue: density ratio, viscosity ratio, Weber number We=4U02λ/νA, and Reynolds number Re=2U0λ/νA, where U0 is the streaming velocity of the top and bottom walls. To minimize the compressibility effects, the Mach number is set to satisfy Ma=U0/cs≤0.05. Besides, the dimensionless time is measured as t∗=2U0t/λ.Figure 2
Schematic diagram for Kelvin-Helmholtz instability.As specified in Ref. [12], we prescribed a vortex sheet in the computational domain that the interface is initially perturbed by a uniformly concentrated vorticity distribution ω0x,y in the form of the Dirac δ function as
(21)ω0x,y=δhhx,y,with
(22)hx,y=yλ+0.01sin2πx+yλ,where hx,y is the interface location in the dimensionless form at t∗=0, and the Dirac δ function is approximated with the following smoothed 4-point cosine function [41],
(23)δhh=π+2sinπ2h+1/4−2sinπ2h−1/4/4π,∣h∣<1.5,5π−2π∣h∣−4sinπ2∣h∣−1/4/8π,1.5≤∣h∣≤2.5,0,∣h∣>2.5.Based on the vorticity distribution in the whole domain, we can find the streaming functionΔψ by solving the Poisson equation Δψ=−ω0, and subsequently, obtain the initial velocity field via u0=∇×ψ. In the simulation, the periodic boundary conditions are applied in the streamwise direction, while the top and bottom walls with a streaming velocity are implemented by the general bounce-back scheme [42]. The distribution profile of the order parameter is initialized by
(24)ϕx,y=0.5−0.5tanh2y−λ−0.01λsin2πx/λW,which enables the value of the order parameter to be smooth across the interface. First, the density- and viscosity-matched fluids are considered for grid resolution study by comparing the results with those of the AMR scheme [12]. As shown in Figure 3, a symmetric interface roll-up is gradually formed with Re=10000 and We=400, and both vorticity contours and interface locations agree well with the AMR results based on the grids of 512×1024.Figure 3
Vorticity contours (color) and interface locations (red line) of AMR scheme [12] and present model with Re=10000, We=400, ρA=ρB, and νA=νB.The further study focuses on the shear-layer instability of fluids with a fixed density contrast but different viscosity ratios (rν=νB/νA=1, 2, 4, and 8). The Reynolds number is fixed at Re=10000 with Weber number being We=400. Figure 4 shows the corresponding results of temporal interface evolutions. As it can be seen, the interface roll-ups are no longer symmetric compared to those of the density-matched fluids, and such interface asymmetry caused by density contrast was also reported in the work of Tauber et al. [43]. Besides, by increasing the viscosity ratio from 1 to 8, the roll-ups of the growing mixing layer are suppressed. This phenomenon can be reflected by the time-varying kinetic energy in the system. Here, the kinetic energy along x- and y-directions are defined as Ex=ρux2/2 and Ey=ρuy2/2, respectively, and their peak values (Ex∣max and Ey∣max) versus the dimensionless time t∗ are portrayed in Figure 5. It is observed that both Ex∣max and Ey∣max, which represent the interacting strength of two different fluids [44], reduce with the increasing of the viscosity ratio, thus leading to the weakened interface roll-ups.Figure 4
Snapshots of interface evolutions for different viscosity ratios (Re=10000, We=400, and ρA/ρB=2).Figure 5
Temporal evolutions of the peak kinetic energy for different viscosity ratios.
(a)(b)We continue with our investigations by taking the effects of the density ratio and Reynolds number into account. In this part, three density ratios (rρ=ρA/ρB=10, 100, and 1000) are examined at different Reynolds numbers varying from 2000 to 10000, which intends to mimic the operating condition in prefilming primary atomization. As shown in Figure 6, a thick tip of the heavy fluid is stretched into a narrow ligament at the interface under low or moderate density ratios (rρ≤100). When the Reynolds number is high enough (Re=10000), the ligament tends to detach from the heavy fluid. However, such a breakup pattern is not observed in the density-matched fluid system. This is due to that the inertial force of the heavy fluid comes to dominate the viscous force with the increasing density ratio and Reynolds number, causing a distorted interface and its subsequent disintegration. On the other hand, once the density ratio and Reynolds number both reach a high level, small Kelvin-Helmholtz wavelets are formed on large-scale waves of the perturbed interface, which trigger the chaotic interface evolution as shown in Figure 6(b) [45]. These severe topological changes at the interface are quite in line with those observed in the primary breakup of liquid sheets [1, 46].Figure 6
Snapshots of KHI for various density ratios and Reynolds numbers (We=400, νA=νB, and t∗=4.0).From a different perspective, the effects of the density ratio and Reynolds number on the interface disturbance are quantitatively assessed by the averaged order parameter (ϕavg) with its definition given by ϕavg=∫0λϕx,ydx/λ. Taking the case of Re=2000 and rρ=10 as an example, Figure 7 depicts the profiles of ϕavg against the y-direction at different times, in which the y-coordinate is measured by y∗=y+λ/2λ. As seen, the profile keeps to be smooth at the initial time, and then, it oscillates due to the mixing of two different fluids. Therefore, the symbol dmix marked in Figure 7 can be roughly regarded as the width of the mixing layer. We then record the values of dmix under different density ratios and Reynolds numbers in Table 1. Evidently, the width of the mixing layer increases with the Reynolds number, yet a reversed trend is observed by increasing the density ratio. These results are consistent with those displayed in Figure 6. Furthermore, a linear analysis has been carried out for KHI in wide ranges of the density ratio and Reynolds number. It is found out that the interface perturbation has an exponential growth at the initial linear increasing stage of KHI, which greatly depends on the density ratio and the Reynolds number. To be specific, the perturbation growth is accelerated by increasing the Reynolds number but it is suppressed by increasing the density ratio, which is in line with the numerical results. However, especially in the case of the large density ratio and high Reynolds number, the perturbation grows nonlinearly to form small Kelvin-Helmholtz wavelets at the later stage of KHI. The linear theory fails in this nonlinear regime, and the numerical results in this paper become useful alternatives.Figure 7
Averaged order parameter profiles against the y-direction at different dimensionless times (Re=2000, We=400, and rρ=10).Table 1
Widths of the mixing layer for different density ratios and Reynolds numbers.
Redmixrρ=10rρ=50rρ=100rρ=100020000.2800.2500.2030.16250000.4270.3920.3070.266100000.4780.4670.3410.269


### 3.3. Rayleigh-Taylor Instability

Figure8 sketches the other incompressible two-phase flow system involved in this paper, where a layer of heavy fluid with density ρA is located on top of the light one with density ρB. As stated in Ref. [21], any disturbance at the interface will be accelerated by gravity to produce downward-moving spikes of heavy fluid and upward-moving bubbles of light fluid. This is the so-called RTI, a crucial type of instability that is responsible for the interface destabilization. In this section, we pay special attention to RTI regarding fluids with different viscosity ratios and density ratios in a wide range of the Reynolds number. The simulation is implemented in a 2D domain of size d×4d with periodic and no-slip boundary conditions in the streamwise and normal-wall directions, respectively. The numerical results are presented in terms of dimensionless parameters, including the viscosity ratio, Atwood number At=ρA−ρB/ρA+ρB, and Reynolds number Re=dgd/νA. Here, g is the gravity acceleration. The time is normalized as t∗=t/d/g.Figure 8
Schematic diagram for Rayleigh-Taylor instability.Numerical study is started with two benchmark cases for grid resolution study and further verification. In the first case of viscosity-matched fluids, the initial order parameter profile, that smoothly across the interface with an amplitude of0.1d, is given as
(25)ϕx,y=0.5+0.5tanh2y−2d−0.1dcos2πx/dW.As specified in Ref. [21], the key parameters are fixed as Re=2048, At=0.5, and gd=0.04. After that, the temporal evolution of the perturbed interface is shown in Figure 9 based on the grids of 512×2048. It is observed that the heavy fluid is downward accelerated by gravity, and it gradually penetrates into the light fluid. At t∗=2.0, the front-end of the heavy fluid evolves to be an umbrella-like shape. As the light fluid is upward pressured, the velocity difference between two fluids triggers KHI, leading to the counter-rotating vortices of the heavy fluid (t∗=3.0). At a later time (t∗=5.0), these two unstable vortices arouse a pair of secondary vortices at the tips of the roll-ups. The obtained interface evolution compares well with those reported in Refs. [21, 24]. Front-positions of the spike (hs) and bubble (hb) measured in units of d also agree with the results of He et al. [21], as shown in Figure 10.Figure 9
Snapshots of interface evolution forRe=2048 and At=0.5.
(a)(b)(c)(d)(e)Figure 10
Comparison of the bubble and spike front-positions between present results and those of He et al. [21].The first test confirms the grid resolution and ability of the present model in capturing the feature of RTI under a low density ratio, but the reliability for those with large density ratios at high Reynolds numbers is not sufficiently verified. Thus, in the second benchmark case, we intend to collect more evidence by comparing our simulation with the analytical data [47]. It is argued in Refs. [21, 47] that the slightly perturbed interface has an exponential growth in the early stage,
(26)a=a0eαt,where α is the growth rate, and a is the temporal interface amplitude with its initial value of a0. For the viscosity-matched fluids with negligible surface tension, the growth rate will be a function of Atwood number and wavenumber (k=2π/d). For presenting the results, the dimensionless growth rate and wavenumber are defined as α∗=α/g2/ν3 and k∗=k/g/ν23, respectively. Figure 11 depicts the results of α∗ obtained from the present model at At=0.998 with 0.01≤k∗≤2.0, as well as those predicted by linear theory [47]. In this case, the grids are 512×2048 as well, but the initial interface amplitude is set with a much small value of a0=0.01d. Noting that in this case, At=0.998 represents a large density ratio of ρA/ρB=1000, and 0.01≤k∗≤2.0 indicates the Reynolds number varying in the range of 5.6≤Re≤15749.6. Compared with the numerical results in Ref. [48], our results are in better agreement with the analytical data [47]. The above results once more demonstrate the ability of the present model in dealing with challenging conditions and also affirm the grid resolution for the following study.Figure 11
Comparison of the dimensionless growth rate obtained from the numerical methods (symbol) with those predicted by linear theory [47] (solid line).In the above two benchmark cases, the simulation focuses on viscosity-matched fluids. However, the actual fluids generally own a viscosity contrast. Hence, we assess the effects of viscosity ratio (rν=νB/νA) on the interface evolution of RTI in this part. To reduce the influence of density ratio, the Atwood number is fixed at At=0.5. The numerical study then considers four viscosity ratios, including rν=1, 2, 4, and 8 at Re=256 and Re=2048, respectively.Figure12 shows the snapshots of the interface shape of different viscosity ratios for Re=256 at t∗=5.0. Compared with those in Figure 13 for Re=2048, the heavy fluid at a low Reynolds number still rolls up as two side spikes, but the spiral vortices are not observed. By increasing the viscosity ratio, it is found out that roll-ups of the heavy fluid and the subsequent spiral vortices are all suppressed. This is due to the enhanced viscous force of the light fluid that reduces the strength of KHI, which is consistent with the results reported in Section 3.2. As plotted in Figures 14 and 15, the viscosity ratio also demonstrates distinct impacts on bubble position and velocity. Here, the bubble velocity ub is measured in units of Atgd/1+At. As seen, the motion of the bubble is restrained with the increasing viscosity ratio and such effects are more evident at a low Reynolds number. The possible reason is that the decrease of vortical effects stemmed from KHI results in a poor acceleration of bubbles. Such an adverse impact on the bubble motion happens to be more remarkable when the viscous force dominates the inertia force at a low Reynolds number (Re=256).Figure 12
Snapshots of interface shape for different viscosity ratios (Re=256, At=0.5, and t∗=5.0).
(a)(b)(c)(d)Figure 13
Snapshots of interface shape for different viscosity ratios (Re=2048, At=0.5, and t∗=5.0).
(a)(b)(c)(d)Figure 14
Temporal variations of the normalized bubble position and velocity for different viscosity ratios (Re=256 and At=0.5).
(a)(b)Figure 15
Temporal variations of the normalized bubble position and velocity for different viscosity ratios (Re=2048 and At=0.5).
(a)(b)Ultimately, the effects of density ratio (rρ=ρA/ρB) on the interface evolution of RTI are investigated with the viscosity ratio fixed at rν=2. The following study examines four different density ratios (rρ=3, 10, 100, and 1000) at Re=2048 and Re=10000. As shown in Figure 16, the snapshots of interface at Re=2048 for different density ratios are captured when the spike-front is close to the bottom wall. For the low density ratios (rρ=3 and 10), we can observe the umbrella-like spikes with spiral vortices appearing at the tails of the side roll-ups. By increasing the density ratio to rρ=, the roll-ups of heavy fluid shrink to be small fingers, which is in line with the evolving trend predicted in Ref. [21] and the behavior of KHI at moderate density ratios [22]. Once the density ratio reaches a sufficiently large value (rρ=1000), the umbrella-like spike-front evolve to be a rocket-like configuration with small serrated waves growing on the spike sides. These small serrated waves are born out of Kelvin-Helmholtz wavelets, which are exactly the characteristics of KHI at large density ratios [45]. Such topological interface deformations are more pronounced at a higher Reynolds number (Re=10000). As depicted in Figures 17(c) and 17(d), the small-scale ligaments are stretched and then tend to detach from the spikes due to the dominated effects of inertia force under present challenging conditions. This kind of severe breakup pattern is generally observed in the liquid jet at the large density ratio and high Reynolds number [1, 48]. Likewise, we compare the temporal variations of the normalized bubble position and velocity for different density ratios in Figure 18. Apparently, the bubble of light fluid under a larger density ratio grows much faster than those under lower density ratios.Figure 16
Snapshots of interface shape for different density ratios (Re=2048 and rν=2).
(a)(b)(c)(d)Figure 17
Snapshots of interface shape for different density ratios (Re=10000 and rν=2).
(a)(b)(c)(d)Figure 18
Temporal variations of the normalized bubble position and velocity for different density ratios (Re=10000 and rν=2).
(a)(b)


## 3.1. Static Droplet

Before the simulation of KHI and RTI, the present model is preliminarily tested by Laplace’s law as well as the sensitivity study to the values of mobility. Initially, a static droplet with the radiusR=25 surrounded by the gas phase is placed at the center of a periodic domain with 150×150 grid points. To mimic the challenging condition of the large density ratio and high Reynolds number, the density ratio and kinematic viscosity of liquid and gas are set as ρA/ρB=1000 and νA=νB=0.002, respectively. According to Laplace’s law, the theoretical pressure difference between liquid and gas is Δp=σ/R at steady-state. Figure 1(a) compares the theoretical and numerical values of pressure difference for various surface tension coefficients, which are in good agreement. In the simulation, we fix s0f=s7f=s8f=1.0, s1f=s2f, and s4f=s6f as usual. It is found out that 0.75≤s1f=s2f≤1.2 and 0.6≤s4f=s6f≤1.3 all can give satisfying results. For simplicity, s1f=s2f=s4f=s6f=1.1 is set. Regarding the relaxation elements in Sg, similar values are adopted as those in Ref. [22], except s0g=s3g=s5g=1.0 and s7g=s8g related to kinematic viscosity.Figure 1
Validation and sensitivity study. (a) Comparison of the pressure difference between theoretical and numerical results. (b) Maximum kinetic energy for various mobility.
(a)(b)For the importance of mobility in the present model, a sensitivity study to its values is necessary. Figure1(b) depicts the maximum kinetic energy Emax in the system for various mobilities. Here, the kinetic energy and the dimensionless time are defined as E=ρu2/2 and t∗=tσ/ρBνBR, respectively. It is noteworthy that the present MRT model can successfully damp high-frequency oscillations of the kinetic energy at the early stages. Besides, the undervalue of mobility allows less damping time and lower kinetic energy, which also means smaller spurious currents (10−9) in the system. Taking the numerical stability into account as well, the mobility is selected as a reasonable value of M=0.01 in the following study.


## 3.2. Kelvin-Helmholtz Instability

In this section, KHI in a two-phase incompressible flow system is considered with the schematic diagram depicted in Figure2. The shear flow of two immiscible fluids with densities ρA and ρB takes place in a 2D channel of size (λ×2λ), whose top and bottom walls move in opposite directions. Four dimensionless groups are chosen to regulate the present issue: density ratio, viscosity ratio, Weber number We=4U02λ/νA, and Reynolds number Re=2U0λ/νA, where U0 is the streaming velocity of the top and bottom walls. To minimize the compressibility effects, the Mach number is set to satisfy Ma=U0/cs≤0.05. Besides, the dimensionless time is measured as t∗=2U0t/λ.Figure 2
Schematic diagram for Kelvin-Helmholtz instability.As specified in Ref. [12], we prescribed a vortex sheet in the computational domain that the interface is initially perturbed by a uniformly concentrated vorticity distribution ω0x,y in the form of the Dirac δ function as
(21)ω0x,y=δhhx,y,with
(22)hx,y=yλ+0.01sin2πx+yλ,where hx,y is the interface location in the dimensionless form at t∗=0, and the Dirac δ function is approximated with the following smoothed 4-point cosine function [41],
(23)δhh=π+2sinπ2h+1/4−2sinπ2h−1/4/4π,∣h∣<1.5,5π−2π∣h∣−4sinπ2∣h∣−1/4/8π,1.5≤∣h∣≤2.5,0,∣h∣>2.5.Based on the vorticity distribution in the whole domain, we can find the streaming functionΔψ by solving the Poisson equation Δψ=−ω0, and subsequently, obtain the initial velocity field via u0=∇×ψ. In the simulation, the periodic boundary conditions are applied in the streamwise direction, while the top and bottom walls with a streaming velocity are implemented by the general bounce-back scheme [42]. The distribution profile of the order parameter is initialized by
(24)ϕx,y=0.5−0.5tanh2y−λ−0.01λsin2πx/λW,which enables the value of the order parameter to be smooth across the interface. First, the density- and viscosity-matched fluids are considered for grid resolution study by comparing the results with those of the AMR scheme [12]. As shown in Figure 3, a symmetric interface roll-up is gradually formed with Re=10000 and We=400, and both vorticity contours and interface locations agree well with the AMR results based on the grids of 512×1024.Figure 3
Vorticity contours (color) and interface locations (red line) of AMR scheme [12] and present model with Re=10000, We=400, ρA=ρB, and νA=νB.The further study focuses on the shear-layer instability of fluids with a fixed density contrast but different viscosity ratios (rν=νB/νA=1, 2, 4, and 8). The Reynolds number is fixed at Re=10000 with Weber number being We=400. Figure 4 shows the corresponding results of temporal interface evolutions. As it can be seen, the interface roll-ups are no longer symmetric compared to those of the density-matched fluids, and such interface asymmetry caused by density contrast was also reported in the work of Tauber et al. [43]. Besides, by increasing the viscosity ratio from 1 to 8, the roll-ups of the growing mixing layer are suppressed. This phenomenon can be reflected by the time-varying kinetic energy in the system. Here, the kinetic energy along x- and y-directions are defined as Ex=ρux2/2 and Ey=ρuy2/2, respectively, and their peak values (Ex∣max and Ey∣max) versus the dimensionless time t∗ are portrayed in Figure 5. It is observed that both Ex∣max and Ey∣max, which represent the interacting strength of two different fluids [44], reduce with the increasing of the viscosity ratio, thus leading to the weakened interface roll-ups.Figure 4
Snapshots of interface evolutions for different viscosity ratios (Re=10000, We=400, and ρA/ρB=2).Figure 5
Temporal evolutions of the peak kinetic energy for different viscosity ratios.
(a)(b)We continue with our investigations by taking the effects of the density ratio and Reynolds number into account. In this part, three density ratios (rρ=ρA/ρB=10, 100, and 1000) are examined at different Reynolds numbers varying from 2000 to 10000, which intends to mimic the operating condition in prefilming primary atomization. As shown in Figure 6, a thick tip of the heavy fluid is stretched into a narrow ligament at the interface under low or moderate density ratios (rρ≤100). When the Reynolds number is high enough (Re=10000), the ligament tends to detach from the heavy fluid. However, such a breakup pattern is not observed in the density-matched fluid system. This is due to that the inertial force of the heavy fluid comes to dominate the viscous force with the increasing density ratio and Reynolds number, causing a distorted interface and its subsequent disintegration. On the other hand, once the density ratio and Reynolds number both reach a high level, small Kelvin-Helmholtz wavelets are formed on large-scale waves of the perturbed interface, which trigger the chaotic interface evolution as shown in Figure 6(b) [45]. These severe topological changes at the interface are quite in line with those observed in the primary breakup of liquid sheets [1, 46].Figure 6
Snapshots of KHI for various density ratios and Reynolds numbers (We=400, νA=νB, and t∗=4.0).From a different perspective, the effects of the density ratio and Reynolds number on the interface disturbance are quantitatively assessed by the averaged order parameter (ϕavg) with its definition given by ϕavg=∫0λϕx,ydx/λ. Taking the case of Re=2000 and rρ=10 as an example, Figure 7 depicts the profiles of ϕavg against the y-direction at different times, in which the y-coordinate is measured by y∗=y+λ/2λ. As seen, the profile keeps to be smooth at the initial time, and then, it oscillates due to the mixing of two different fluids. Therefore, the symbol dmix marked in Figure 7 can be roughly regarded as the width of the mixing layer. We then record the values of dmix under different density ratios and Reynolds numbers in Table 1. Evidently, the width of the mixing layer increases with the Reynolds number, yet a reversed trend is observed by increasing the density ratio. These results are consistent with those displayed in Figure 6. Furthermore, a linear analysis has been carried out for KHI in wide ranges of the density ratio and Reynolds number. It is found out that the interface perturbation has an exponential growth at the initial linear increasing stage of KHI, which greatly depends on the density ratio and the Reynolds number. To be specific, the perturbation growth is accelerated by increasing the Reynolds number but it is suppressed by increasing the density ratio, which is in line with the numerical results. However, especially in the case of the large density ratio and high Reynolds number, the perturbation grows nonlinearly to form small Kelvin-Helmholtz wavelets at the later stage of KHI. The linear theory fails in this nonlinear regime, and the numerical results in this paper become useful alternatives.Figure 7
Averaged order parameter profiles against the y-direction at different dimensionless times (Re=2000, We=400, and rρ=10).Table 1
Widths of the mixing layer for different density ratios and Reynolds numbers.
Redmixrρ=10rρ=50rρ=100rρ=100020000.2800.2500.2030.16250000.4270.3920.3070.266100000.4780.4670.3410.269


## 3.3. Rayleigh-Taylor Instability

Figure8 sketches the other incompressible two-phase flow system involved in this paper, where a layer of heavy fluid with density ρA is located on top of the light one with density ρB. As stated in Ref. [21], any disturbance at the interface will be accelerated by gravity to produce downward-moving spikes of heavy fluid and upward-moving bubbles of light fluid. This is the so-called RTI, a crucial type of instability that is responsible for the interface destabilization. In this section, we pay special attention to RTI regarding fluids with different viscosity ratios and density ratios in a wide range of the Reynolds number. The simulation is implemented in a 2D domain of size d×4d with periodic and no-slip boundary conditions in the streamwise and normal-wall directions, respectively. The numerical results are presented in terms of dimensionless parameters, including the viscosity ratio, Atwood number At=ρA−ρB/ρA+ρB, and Reynolds number Re=dgd/νA. Here, g is the gravity acceleration. The time is normalized as t∗=t/d/g.Figure 8
Schematic diagram for Rayleigh-Taylor instability.Numerical study is started with two benchmark cases for grid resolution study and further verification. In the first case of viscosity-matched fluids, the initial order parameter profile, that smoothly across the interface with an amplitude of0.1d, is given as
(25)ϕx,y=0.5+0.5tanh2y−2d−0.1dcos2πx/dW.As specified in Ref. [21], the key parameters are fixed as Re=2048, At=0.5, and gd=0.04. After that, the temporal evolution of the perturbed interface is shown in Figure 9 based on the grids of 512×2048. It is observed that the heavy fluid is downward accelerated by gravity, and it gradually penetrates into the light fluid. At t∗=2.0, the front-end of the heavy fluid evolves to be an umbrella-like shape. As the light fluid is upward pressured, the velocity difference between two fluids triggers KHI, leading to the counter-rotating vortices of the heavy fluid (t∗=3.0). At a later time (t∗=5.0), these two unstable vortices arouse a pair of secondary vortices at the tips of the roll-ups. The obtained interface evolution compares well with those reported in Refs. [21, 24]. Front-positions of the spike (hs) and bubble (hb) measured in units of d also agree with the results of He et al. [21], as shown in Figure 10.Figure 9
Snapshots of interface evolution forRe=2048 and At=0.5.
(a)(b)(c)(d)(e)Figure 10
Comparison of the bubble and spike front-positions between present results and those of He et al. [21].The first test confirms the grid resolution and ability of the present model in capturing the feature of RTI under a low density ratio, but the reliability for those with large density ratios at high Reynolds numbers is not sufficiently verified. Thus, in the second benchmark case, we intend to collect more evidence by comparing our simulation with the analytical data [47]. It is argued in Refs. [21, 47] that the slightly perturbed interface has an exponential growth in the early stage,
(26)a=a0eαt,where α is the growth rate, and a is the temporal interface amplitude with its initial value of a0. For the viscosity-matched fluids with negligible surface tension, the growth rate will be a function of Atwood number and wavenumber (k=2π/d). For presenting the results, the dimensionless growth rate and wavenumber are defined as α∗=α/g2/ν3 and k∗=k/g/ν23, respectively. Figure 11 depicts the results of α∗ obtained from the present model at At=0.998 with 0.01≤k∗≤2.0, as well as those predicted by linear theory [47]. In this case, the grids are 512×2048 as well, but the initial interface amplitude is set with a much small value of a0=0.01d. Noting that in this case, At=0.998 represents a large density ratio of ρA/ρB=1000, and 0.01≤k∗≤2.0 indicates the Reynolds number varying in the range of 5.6≤Re≤15749.6. Compared with the numerical results in Ref. [48], our results are in better agreement with the analytical data [47]. The above results once more demonstrate the ability of the present model in dealing with challenging conditions and also affirm the grid resolution for the following study.Figure 11
Comparison of the dimensionless growth rate obtained from the numerical methods (symbol) with those predicted by linear theory [47] (solid line).In the above two benchmark cases, the simulation focuses on viscosity-matched fluids. However, the actual fluids generally own a viscosity contrast. Hence, we assess the effects of viscosity ratio (rν=νB/νA) on the interface evolution of RTI in this part. To reduce the influence of density ratio, the Atwood number is fixed at At=0.5. The numerical study then considers four viscosity ratios, including rν=1, 2, 4, and 8 at Re=256 and Re=2048, respectively.Figure12 shows the snapshots of the interface shape of different viscosity ratios for Re=256 at t∗=5.0. Compared with those in Figure 13 for Re=2048, the heavy fluid at a low Reynolds number still rolls up as two side spikes, but the spiral vortices are not observed. By increasing the viscosity ratio, it is found out that roll-ups of the heavy fluid and the subsequent spiral vortices are all suppressed. This is due to the enhanced viscous force of the light fluid that reduces the strength of KHI, which is consistent with the results reported in Section 3.2. As plotted in Figures 14 and 15, the viscosity ratio also demonstrates distinct impacts on bubble position and velocity. Here, the bubble velocity ub is measured in units of Atgd/1+At. As seen, the motion of the bubble is restrained with the increasing viscosity ratio and such effects are more evident at a low Reynolds number. The possible reason is that the decrease of vortical effects stemmed from KHI results in a poor acceleration of bubbles. Such an adverse impact on the bubble motion happens to be more remarkable when the viscous force dominates the inertia force at a low Reynolds number (Re=256).Figure 12
Snapshots of interface shape for different viscosity ratios (Re=256, At=0.5, and t∗=5.0).
(a)(b)(c)(d)Figure 13
Snapshots of interface shape for different viscosity ratios (Re=2048, At=0.5, and t∗=5.0).
(a)(b)(c)(d)Figure 14
Temporal variations of the normalized bubble position and velocity for different viscosity ratios (Re=256 and At=0.5).
(a)(b)Figure 15
Temporal variations of the normalized bubble position and velocity for different viscosity ratios (Re=2048 and At=0.5).
(a)(b)Ultimately, the effects of density ratio (rρ=ρA/ρB) on the interface evolution of RTI are investigated with the viscosity ratio fixed at rν=2. The following study examines four different density ratios (rρ=3, 10, 100, and 1000) at Re=2048 and Re=10000. As shown in Figure 16, the snapshots of interface at Re=2048 for different density ratios are captured when the spike-front is close to the bottom wall. For the low density ratios (rρ=3 and 10), we can observe the umbrella-like spikes with spiral vortices appearing at the tails of the side roll-ups. By increasing the density ratio to rρ=, the roll-ups of heavy fluid shrink to be small fingers, which is in line with the evolving trend predicted in Ref. [21] and the behavior of KHI at moderate density ratios [22]. Once the density ratio reaches a sufficiently large value (rρ=1000), the umbrella-like spike-front evolve to be a rocket-like configuration with small serrated waves growing on the spike sides. These small serrated waves are born out of Kelvin-Helmholtz wavelets, which are exactly the characteristics of KHI at large density ratios [45]. Such topological interface deformations are more pronounced at a higher Reynolds number (Re=10000). As depicted in Figures 17(c) and 17(d), the small-scale ligaments are stretched and then tend to detach from the spikes due to the dominated effects of inertia force under present challenging conditions. This kind of severe breakup pattern is generally observed in the liquid jet at the large density ratio and high Reynolds number [1, 48]. Likewise, we compare the temporal variations of the normalized bubble position and velocity for different density ratios in Figure 18. Apparently, the bubble of light fluid under a larger density ratio grows much faster than those under lower density ratios.Figure 16
Snapshots of interface shape for different density ratios (Re=2048 and rν=2).
(a)(b)(c)(d)Figure 17
Snapshots of interface shape for different density ratios (Re=10000 and rν=2).
(a)(b)(c)(d)Figure 18
Temporal variations of the normalized bubble position and velocity for different density ratios (Re=10000 and rν=2).
(a)(b)


## 4. Conclusions

Recent studies have revealed that KHI and RTI dominate the processes of liquid sheet breakup and primary atomization. Generally, the liquid and gas own a large density contrast and high Reynolds number in these applications. Hence, the investigations of KHI and RTI by considering the effects of the large density contrast and high Reynolds number have great significance in better understanding the destabilizing mechanism of the liquid layer. Unfortunately, the related investigations of KHI and RTI are rarely conducted in such challenging conditions. To this end, the morphological characteristics of KHI and RTI are inspected in wide ranges of the density ratio, Reynolds number, and viscosity ratio in this paper. First, a lattice Boltzmann MRT model based on the conservative AC equation is reconstructed with necessary modifications, and its numerical ability at challenging conditions is well verified together with sensitive study to the values of mobility. Afterward, the numerical simulations are implemented, and the corresponding results indicate that the viscosity ratio and density ratio tend to suppress the interface disturbance, while the Reynolds number favors the fluids interpenetrating. Nonetheless, when the density ratio and the Reynolds number are both increased to large values, the aroused Kelvin-Helmholtz wavelets trigger severe interface evolutions due to the dominated inertia force of heavy fluid. Noting that the aforementioned results of KHI and RTI, particularly, under a large density ratio and high Reynolds number greatly support the typical mechanisms observed in the primary atomization process.



---

*Source: 8885226-2020-06-19.xml*