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

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 spikespiraling 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.


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][5][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.
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 singlemode 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 multiplerelaxation-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-fieldbased 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. With the MRT collision operator [29], the generalized evolution equation for the conservative AC equation [29,31] can be written as

Mathematical Method
where f i ðx, tÞ and f eq i ðx, 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. c i is the discrete velocity of the ith direction, δt is the time step. Λ International Journal of Aerospace Engineering which M is an orthogonal transformation matrix and S f is a diagonal matrix given by (for the D2Q9 model) In Equation (1), F i represents the forcing term in the velocity space [27,32], and it is defined as where ω i is the weight coefficient with the value given by ω 0 = 4/9, ω 1,⋯4 = 1/9, and ω 5,⋯8 = 1/36, c s = c/ ffiffi ffi 3 p 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 matrix M, the right-hand side of Equation (1) can be projected onto the moment space via m f = Mf and m eq f = Mf eq [33,34] as where f = ð f 0 ,⋯f 8 Þ T , f eq = ðf eq 0 ,⋯f eq 8 Þ T , I is the unit tensor, and F̂= MF = MðF 0 ,⋯,F 8 Þ T is the forcing term in the moment space. Based on f eq and F in the velocity space [27,32], m eq f and F̂can be easily derived for the D2Q9 model as m eq f = ϕ 1,−2, 1, u x ,−u x , u y ,−u y , 0, 0 where u x and u y are the components of the macroscopic velocity u, n x and n y are the components of the unit vector n. Afterward, the streaming process is given as where the postcollision distribution function f * i can be obtained via f * = M −1 m * f . Applying the Chapman-Enskog expansion to Equation (1), the conservative AC equation [30,35] can be recovered correctly as where the mobility M is determined by in which τ f is the dimensionless relaxation time of f j , and it is related to the diagonal matrix S f with s f 3 = s f 5 = 1/τ f . In the present model, the order parameter is calculated as and the fluid density ρ should be consistent with the order parameter by taking the linear interpolation as where ρ A and ρ B stand for the densities of two different fluids corresponding to the order parameter ϕ A = 1 and ϕ B = 0, respectively.

The MRT LB Model for the Navier-Stokes Equations.
In order to improve the numerical stability, the lattice 3 International Journal of Aerospace Engineering Boltzmann equation for the incompressible NS equations is combined with the MRT collision operator as where g i and g eq i are the density distribution function and its corresponding equilibrium distribution function, respectively. R i is the forcing term in the velocity space. Λ g ij is an element of the collision matrix Λ g = M −1 S g M, in which S g is a diagonal matrix given by (for the D2Q9 model) Similar to Equation (3), the right-hand side of Equation (12) in the moment space reads where m g = Mg, m eq g = Mg eq , andR = MR are the corresponding matrices in the moment space with g = ðg 0 ,⋯g 8 Þ T , g eq = ðg eq 0 ,⋯g eq 8 Þ T , and R = ðR 0 ,⋯,R 8 Þ T . Inspired by the BGK model of Liang et al. [27], the matrices m eq g andR for the D2Q9 model are evaluated as where p is the hydrodynamic pressure. G x and G y are the components of the total force G = F s + F e + F a . Here, F s = −σα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. F e 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 F a = q a u is introduced with the term q a determined by Ultimately, with g * = M −1 m * g , the streaming process is achieved as Based on the density distribution function g i , the macroscopic quantities u and p can be evaluated as Besides, through the Chapman-Enskog analysis, the kinematic viscosity ν is determined by where τ g = 1/s g 7 = 1/s g 8 . 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, 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 s g 7 and s g 8 of the diagonal matrix S g 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].

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 4 International Journal of Aerospace Engineering static droplet with the radius R = 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. For the importance of mobility in the present model, a sensitivity study to its values is necessary. Figure 1(b) depicts the maximum kinetic energy E max in the system for various mobilities. Here, the kinetic energy and the dimensionless time are defined as E = ρjuj 2 /2 and t * = tσ/ðρ B ν B RÞ, 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.

Kelvin-Helmholtz Instability.
In this section, KHI in a two-phase incompressible flow system is considered with the schematic diagram depicted in Figure 2. 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 = 4U 2 0 λ/ν A , and Reynolds number Re = 2U 0 λ/ν A , where U 0 is the streaming velocity of the top and bottom walls. To minimize the compressibility effects, the Mach number is set to satisfy Ma = U 0 /c s ≤ 0:05. Besides, the dimensionless time is measured as t * = 2U 0 t/λ.
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 ω 0 ðx, yÞ in the form of the Dirac δ function as with where hðx, 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],  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 u 0 = ∇ × ψ. 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 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.
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 densitymatched 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 E x = ρu 2 x /2 and E y = ρu 2 y /2, respectively, and their peak values (E x|max and E y|max ) versus the dimensionless time t * are portrayed in Figure 5. It is observed that both E x|max and E y|max , which represent the interacting strength of two different fluids [44], reduce with the increas-ing of the viscosity ratio, thus leading to the weakened interface roll-ups.
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].
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, yÞdx/λ. 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 d mix marked in Figure 7 can be roughly regarded as the width of the mixing layer. We then record the values of d mix 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 8 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  International Journal of Aerospace Engineering 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

Rayleigh-Taylor Instability.
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 = ffiffiffiffiffi ffi dg p d/ν A . Here, g is the gravity acceleration. The time is normalized as t * = t/ ffiffiffiffiffiffiffi d/g p . 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 of 0:1d, is given as As specified in Ref. [21], the key parameters are fixed as Re = 2048, At = 0:5, and ffiffiffiffiffi ffi gd p = 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 counterrotating 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 (h s ) and bubble (h b ) measured in units of d also agree with the results of He et al. [21], as shown in Figure 10.
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 AMR results Figure 3: Vorticity contours (color) and interface locations (red line) of AMR scheme [12] and present model with Re = 10000, We = 400, 7 International Journal of Aerospace Engineering the slightly perturbed interface has an exponential growth in the early stage,

Present results
where α is the growth rate, and a is the temporal interface amplitude with its initial value of a 0 . For the viscositymatched 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 α * = α/ ffiffiffiffiffiffiffiffiffi g 2 /ν 3 p and k * = k/ ffiffiffiffiffiffiffiffiffi g/ν 2 3 p , 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 a 0 = 0:01d. Noting that in this case, At = 0:998 represents a   [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.
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. Figure 12 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 u b is measured in units of ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Atgd/ð1 + AtÞ p . 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.    Figure 11: Comparison of the dimensionless growth rate obtained from the numerical methods (symbol) with those predicted by linear theory [47] (solid line). 9 International Journal of Aerospace Engineering 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).
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 Figure 12: Snapshots of interface shape for different viscosity ratios (Re = 256, At = 0:5, and t * = 5:0). Figure 13: Snapshots of interface shape for different viscosity ratios (Re = 2048, At = 0:5, and t * = 5:0).

10
International Journal of Aerospace Engineering 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.

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.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.