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.

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 [46]. 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 [1115]. 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 (). 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 where and are the particle distribution function and its corresponding equilibrium distribution function for the order parameter at position and time , respectively. Here, the order parameter taking 1 and 0 is adopted to distinguish different fluids with the interface marked by the contour level of . is the discrete velocity of the direction, is the time step. stands for an element of the collision matrix , in which is an orthogonal transformation matrix and is a diagonal matrix given by (for the D2Q9 model)

In Equation (1), represents the forcing term in the velocity space [27, 32], and it is defined as where is the weight coefficient with the value given by , , and , is the sound speed, and is the macroscopic velocity. denotes the unit vector normal to the interface. is a function of , and is the interface thickness. The time derivative term is introduced to eliminate the artificial term in the recovered equation.

Through the orthogonal transformation matrix , the right-hand side of Equation (1) can be projected onto the moment space via and [33, 34] as where , , is the unit tensor, and is the forcing term in the moment space. Based on and in the velocity space [27, 32], and can be easily derived for the D2Q9 model as where and are the components of the macroscopic velocity , and are the components of the unit vector . Afterward, the streaming process is given as where the postcollision distribution function can be obtained via .

Applying the Chapman-Enskog expansion to Equation (1), the conservative AC equation [30, 35] can be recovered correctly as where the mobility is determined by in which is the dimensionless relaxation time of , and it is related to the diagonal matrix with . 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 and stand for the densities of two different fluids corresponding to the order parameter and , 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 where and are the density distribution function and its corresponding equilibrium distribution function, respectively. is the forcing term in the velocity space. is an element of the collision matrix , in which 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 , , and are the corresponding matrices in the moment space with , , and . Inspired by the BGK model of Liang et al. [27], the matrices and for the D2Q9 model are evaluated as where is the hydrodynamic pressure. and are the components of the total force . Here, is the surface tension force with the formulation recommended by Kim [36] and Ren et al. [37], where is the surface tension coefficient. 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 is introduced with the term determined by

Ultimately, with , the streaming process is achieved as

Based on the density distribution function , the macroscopic quantities and can be evaluated as

Besides, through the Chapman-Enskog analysis, the kinematic viscosity is determined by where . 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 and are the kinematic viscosities of two different fluids corresponding to the order parameter and , respectively. Hence, the elements and of the diagonal matrix 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 radius surrounded by the gas phase is placed at the center of a periodic domain with 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 and , respectively. According to Laplace’s law, the theoretical pressure difference between liquid and gas is 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 , , and as usual. It is found out that and all can give satisfying results. For simplicity, is set. Regarding the relaxation elements in , similar values are adopted as those in Ref. [22], except and related to kinematic viscosity.

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 in the system for various mobilities. Here, the kinetic energy and the dimensionless time are defined as and , 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 () in the system. Taking the numerical stability into account as well, the mobility is selected as a reasonable value of 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 Figure 2. The shear flow of two immiscible fluids with densities and takes place in a 2D channel of size (), 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 , and Reynolds number , where is the streaming velocity of the top and bottom walls. To minimize the compressibility effects, the Mach number is set to satisfy . Besides, the dimensionless time is measured as .

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 in the form of the Dirac function as with where is the interface location in the dimensionless form at , 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 , and subsequently, obtain the initial velocity field via . 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 and , and both vorticity contours and interface locations agree well with the AMR results based on the grids of .

The further study focuses on the shear-layer instability of fluids with a fixed density contrast but different viscosity ratios (1, 2, 4, and 8). The Reynolds number is fixed at with Weber number being . 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 - and -directions are defined as and , respectively, and their peak values ( and ) versus the dimensionless time are portrayed in Figure 5. It is observed that both and , 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.

We continue with our investigations by taking the effects of the density ratio and Reynolds number into account. In this part, three density ratios (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 (). When the Reynolds number is high enough (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 () with its definition given by . Taking the case of and as an example, Figure 7 depicts the profiles of against the -direction at different times, in which the -coordinate is measured by . 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 marked in Figure 7 can be roughly regarded as the width of the mixing layer. We then record the values of 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.

3.3. Rayleigh-Taylor Instability

Figure 8 sketches the other incompressible two-phase flow system involved in this paper, where a layer of heavy fluid with density is located on top of the light one with density . 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 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 , and Reynolds number . Here, is the gravity acceleration. The time is normalized as .

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 , is given as

As specified in Ref. [21], the key parameters are fixed as , , and . After that, the temporal evolution of the perturbed interface is shown in Figure 9 based on the grids of . It is observed that the heavy fluid is downward accelerated by gravity, and it gradually penetrates into the light fluid. At , 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 (). At a later time (), 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 () and bubble () measured in units of 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 the slightly perturbed interface has an exponential growth in the early stage, where is the growth rate, and is the temporal interface amplitude with its initial value of . For the viscosity-matched fluids with negligible surface tension, the growth rate will be a function of Atwood number and wavenumber (). For presenting the results, the dimensionless growth rate and wavenumber are defined as and , respectively. Figure 11 depicts the results of obtained from the present model at with , as well as those predicted by linear theory [47]. In this case, the grids are as well, but the initial interface amplitude is set with a much small value of . Noting that in this case, represents a large density ratio of , and indicates the Reynolds number varying in the range of . 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.

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 () on the interface evolution of RTI in this part. To reduce the influence of density ratio, the Atwood number is fixed at . The numerical study then considers four viscosity ratios, including 1, 2, 4, and 8 at and , respectively.

Figure 12 shows the snapshots of the interface shape of different viscosity ratios for at . Compared with those in Figure 13 for , 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 is measured in units of . 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 ().

Ultimately, the effects of density ratio () on the interface evolution of RTI are investigated with the viscosity ratio fixed at . The following study examines four different density ratios (3, 10, 100, and 1000) at and . As shown in Figure 16, the snapshots of interface at for different density ratios are captured when the spike-front is close to the bottom wall. For the low density ratios (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 , 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 (), 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 (). 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.

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.

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.

Acknowledgments

One of the authors gratefully acknowledges insightful discussions with Dr. Hong Liang. This work was supported by the National Natural Science Foundation of China under Grant No. (51776031) and the Key Project of National Science Foundation of Liaoning Province of China under Grant No. (20170540182).