Influence of Cross Perturbations on Turbulent Kelvin–Helmholtz Instability

: Kelvin–Helmholtz instability has been studied extensively in 2D. This study attempts to address the influence of turbulent flow and cross perturbation on the growth rate of the instability and the development of mixing layers in 3D by means of direct numerical simulation. Two perfect gases are considered to be working fluids moving as opposite streams, inducing shear instability at the interface between the fluids and resulting in Kelvin–Helmholtz instability. The results show that cross perturbation affects the instability by increasing the amplitude growth while adding turbulence has almost no effect on the amplitude growth. Furthermore, by increasing the turbulence intensity, a more distinct presence of the inner flow can be seen in the mixing layer of the two phases, and the presence of turbulence expands the range of high-frequency motion significantly due to turbulence structures. The results give a basis for which 3D Kelvin–Helmholtz phenomena should be further investigated using numerical simulation for predictive modeling, beyond the use of simplified 2D theoretical models.


Introduction
Kelvin-Helmholtz instability (KHI) is a classic fluid dynamical problem relevant to many industrial applications and natural phenomena, such as interfacial instability of two flow streams in fuel injection systems, turbulence transition, and breaking waves in the ocean.The phenomena was first described by Helmholtz in 1868 [1].Kelvin-Helmholtz instability was first studied through experiments in the late 1960s [2].This instability was then explored assuming predominantly 2D behavior, linearizing the problem using linear stability theory [3,4].Conventionally, the development of the interfacial wave is addressed through a linearized theory of small perturbations of the two-dimensional base flow.The linear stability theory that leads to the theoretical prediction of the most unstable frequency was first carried out assuming inviscid flows [5].This theory was developed for the inviscid case, showing that the introduction of body forces such as gravity and surface tension tend to dampen instabilities, giving a minimum wave number for which instability will result.This is described in depth in the next section.
The theory was then extended to viscous flow.Some early work in KHI was carried out by Ozgen and others which described the effects of viscosity on the stability characteristics, as well as the effects of gravity and surface tension (Praturi et al. [6], Ozgen et al. [7], Sivadas et al. [8], Otto et al. [9], and Obergaulinger [10]).Comparison of theory to experiments was continued with Yoshikawa for large contrast in viscosity [11].Experimental results were compared to linear stability theory by Matas in 2011 [5].Experimental analysis of KHI was continued by Karabyoglu and Petrarolo in engineering applications of shear flow [12][13][14].While it is important to characterize the viscous effects, this study focuses on the contribution of the initial disturbances in the inviscid case in order to compare more directly to linear stability theory.
With the rise of computational power and numerical simulations, direct numerical simulations (DNSs) of interfacial instabilities became viable.Direct numerical simulation has been a great tool in furthering the understanding of this instability phenomenon.Lee and Kim studied KHI with a 2D DNS using a phase field [15].The use of DNS for solving KHI is explored in Atmakidis, which validates different DNS methods against linear stability theory [16].Zhang in 2018 used a 2D DNS to present the effect of three fluid layers on KHI [17].A parametric study of KHI behavior was also developed using 2D simulation results [18].Recent simulations have provided high-resolution details of the interfacial instability development, the interaction between the interfacial wave and the gas stream, and the formation of the two-phase mixing layer [19].Here, the 'mixing layer' is defined as the region near the fluid interface where one fluid is thoroughly interspersed with the other fluid [20].Therefore, it is mechanically mixed, and no diffusion is assumed between fluids.
After validation of KHI simulations in 2D, this phenomenon can be investigated in more detail in 3D.Linear stability theory deals with simplifications that restrict it to a 2D problem, while turbulence is a predominantly 3D phenomenon.Rogers and Moser 1991 [21], Martinez 2006 [22], Soetrisno 1990 [23], and Wang 2016 [24] analyzed 3D KHI simulations.The study of turbulence in relation to KHI has only recently been explored [20,25].
Kelvin-Helmholtz instability occurs when two parallel fluid streams have a difference in velocities, which induces shear instability at the interface between the fluids.This instability creates roll-up structures that eventually dissipate due to turbulent effects.Although the simple Kelvin-Helmholtz instability behavior has long been an active field of study, comparatively less is known about how these Kelvin-Helmholtz instability structures can be changed by modifying initial flow conditions in three dimensions.This would be helpful in developing predictive modeling of such instabilities in turbulent conditions.Other parametric studies have been conducted on the effects of Kelvin-Helmholtz instability, especially in 2D [26], while experimental and simulation studies have been conducted over the years to understand the roll-up structures that develop as a result [12][13][14].However, the interest of this study lies in the initiation of Kelvin-Helmholtz instability and which parameters, such as interface shape and velocity profiles, can change the development of Kelvin-Helmholtz instability features in three dimensions; this can only be studied with DNSs.This would be useful in predictive models since the 2D theoretical description might not correctly describe the intricacies of instability development with the addition of cross-wise disturbances.Parametric studies tend to focus on laminar cases, while the application of Kelvin-Helmholtz shear instability in a turbulent environment can be complex compared to the 2D theory.This study looks at the addition of three-dimensional interface shape changes as well as the addition of turbulence in the stream-wise driving velocity.In this study, a DNS is conducted to evaluate the effect of cross perturbation in combination with turbulent structures on the instability mechanism and the mixing layer between two streams.Here, two perfect gases are considered as working fluids moving as opposite streams.

Linear Stability Theory
Linear stability theory for Kelvin-Helmholtz was first introduced by Chandrasekhar [3] and is also explained by Drazin [4].The theory is developed assuming incompressible, inviscid flow restricted to velocity only in the x-direction and assuming sinusoidal perturbation of the interface.Linearization of the equations assumes a small initial perturbation, with small defined as k λ < 0.02.Using the method of normal nodes, the solution is therefore assumed to have the following exponential form: where y is the height of the interface from the center line, η 0 is the initial amplitude of the sinusoidal perturbation of the interface, ω is the frequency and comes from the roots of the characteristic equation, and k is the wave number.
where the first term in the square root is the contribution due to shear velocity, the second is the contribution of gravity g, and the third term is the contribution of surface tension σ.
Gravity and surface tension are stabilizing, so a minimum wave number, k min , must be exceeded for the instability to occur.Conversely, this means that for the inviscid case neglecting body forces, the result should always be unstable.
The exponential growth of the initial instability is given by the real part of the ω equation: with the amplitude for the Kelvin-Helmholtz simulations in this study given by the following expression: This solution only holds at the very early stages of the simulation when the sinusoidal shape of the interface is maintained.As soon as cresting of the waves occurs, the amplitude growth slows down, and linear theory no longer applies.This analysis, however, gives a quantifiable way to assess whether the Kelvin-Helmholtz simulations are valid for the VOF solver.

Isotropic Turbulence
To explore the influence of a turbulent initial velocity field on the growth of Kelvin-Helmholtz instability, the 3D case is given a shear velocity with the addition of isotropic synthetic turbulence of the stream-wise velocity.This turbulent field was generated using the Kraichnan method presented by Saad [27].This method initiates the velocity in each cell using a parsed series of 5000 terms.The fluctuation in the x-direction velocity u(x) is given by the following equation: where α is a scaling constant, κ e is the wave number where energy is maximum, and u ′ is the root mean square value of the velocity fluctuations and is adjusted depending on the desired turbulence intensity.
where L is the integral length scale, defined by L ≈ 0.746834/κ e , with κ e , and ν is the kinematic viscosity of the fluid.

Methodology and Simulation Domain
A compressible volume of fluid (VOF) method is used to solve the Euler Equations, neglecting viscous diffusion.
where ⃗ S is the source term vector, p is pressure, and H is total enthalpy defined by This advects the conserved values (αρ g , ρ, ρE, ρ⃗ u), in which α represents the volume fraction of the gas phase (or the fluid of interest in a 2-gas case), ρ g is the density of the gas phase in the domain, ρ is the density of the bulk fluid in any cell within the domain, E is total energy, and ⃗ u is the velocity within the domain.The VOF code is part of an open-source software ABLATE, developed to conduct DNS simulations of the combustion process inside hybrid rocket motors [28].It uses explicit solvers and a stratified flow model to improve scalability, as this allows for a multiphase solution without the reconstruction of the interface.The code is validated in multiple research papers [29,30].
The case considered here is an inviscid shear-driven instability due to the velocity difference between two streams.The coordinate system used is a Cartesian system, with the flow directed in the x-direction, the interface height measured in the y-direction, and the depth defined as the z-direction, as shown in Figure 1.The conventional perturbation is in the y-direction, while the cross perturbation is in the z-direction.The initial interface for the classic Kelvin-Helmholtz case is given by |y| = 0.25 + 0.01sin(10πx) + 0.01sin(10πz) (10) where the second sinusoidal term acts as the additional cross perturbation.The resulting interface shape is shown in Figure 2. The wave amplitude is measured as half the vertical distance from wave crest to trough, thus the initial amplitude A 0 = 0.01 m in the above case, with wave number ω r = 10π.Both fluids are perfect gases (γ = 1.4,R = 287 J kg•K ).The density in the inner region is ρ 1 = 1 kg/m 3 and that of the outer region is ρ 2 = 2 kg/m 3 , with uniform pressure p = 2.5 Pa everywhere.The velocity is 0.5 m/s in the inner region and −0.5 m/s in the outer region.The velocity of ±0.5 m/s was chosen so that the case would be in the regime for instability.According to linear stability theory, for inviscid cases with no additional body forces, any wave number will produce instability, and the regime is always unstable.However, due to numerical diffusion caused by the mesh resolution, a sufficient velocity is necessary for the Kelvin-Helmholtz instability to develop.The 3D domain size is 1 × 1 × 0.2 m with a mesh of 700 × 700 × 100, in streamwise, vertical, and spanwise directions, respectively.Periodic boundary conditions are applied to all boundaries of the domain.In the cases with turbulent initial conditions, an isotropic synthetic turbulent field is generated using synthetic isotropic turbulence of the Kraichnan method [27], with intensities of I = 2% and 4%, where I = u U .All cases are simulated until t = 1 s in order to see the breakdown of the interface into turbulence.

Grid Convergence Study
In order to evaluate the sufficiency of grid sizes to capture the features of the instability at the interface, a convergence study is conducted as the mesh size is increased on a domain of 0.2 m × 1 m × 0.2 m; three computations are conducted on a sequence of three meshes with various grid resolutions: coarse, medium, and fine.These are defined in Table 1.For the generated isotropic turbulence, the integral length scale is L ≈ 0.746834/κ e , with κ e being the wave number where energy is maximum in the energy spectrum [27].Here the isotropic turbulence is generated with an integral length scale of L ≈ 0.014 m, which is significantly larger than the coarse grid size.Therefore, all three grids resolve the integral length scale.Figure 3 shows the interface contours for various meshes, taken at t = 0.1 s.As the contours indicate, the numerical solution tends toward a unique shape as the grid resolution is increased.It is apparent that the coarse grid does not resolve the details of the interface where the curvature at the interface begins to develop, x ≈ 0.04 m.In contrast, the medium and fine meshes resolve curvature adequately; therefore, the medium grid is selected to conduct simulations.

Results
Figure 4 shows the 2D cross section of the 3D simulation results in the x-y plane (z = 0).This shows the expected roll-up behavior of the sinusoidal initial interface, forming wave crests at t = 0.2 s and continuing to swirl at t = 0.3 s.These swirling structures roll onto themselves at t = 0.5 s, as shown in the bottom left image in Figure 5.These later diffuse and break up at t = 0.75 s and t = 1.0 s, as shown in the bottom center and bottom right images in Figure 5.The amplitude as well as kinetic energy was measured at each time step to compare the results between cases.Q-criterion contours were also shown for t = 0.5 s to compare the swirling behavior between the cases with and without cross perturbation.

Effect of Cross Perturbation
Figure 6 shows the growth of wave amplitudes over time for one case with and one without cross perturbation, as well as cases with cross perturbation and the addition of I = 2% and I = 4% turbulence intensity in the initial streamwise velocity.The initial amplitude in the domain is 0.01 m.As predicted by linear stability theory [3,4], the initial growth is exponential, from t = 0 to t = 0.2 s.This is the time at which the wave begins to crest and linear theory breaks down, after which growth decelerates and the transition to turbulence begins, as shown in Figure 4.The definition of 'mixing' as stated by Ling et al. [20] is adopted which refers to the immiscible fluids creating a layer near the interface where one fluid is dispersed into the other fluid phase.At t = 0.2 s, the wave enters the nonlinear phase.The momentum drives the swirling of the vortex instead of contributing to amplitude growth, and eventually, the two-phase mixing layer starts to develop.As evident from the difference in slopes in Figure 6, cross perturbation affects the instability by increasing the amplitude growth, while adding turbulence has almost no effect on the amplitude growth.To better understand the mechanism of instability in the presence of cross perturbation, iso-surfaces of the Q-criterion are plotted in Figure 7 colored by vorticity in the spanwise direction, ω z with units [rad/s].Note that vorticity is negligible in other directions.The Q-criterion allows one to visualize a vortex, and it is defined as where Ω and S are rotation and strain rate tensors, respectively.The counter-rotating vortices of upper and lower waves are apparent from the red and blue contours of spanwise vorticity, respectively.The roll wave structures are thickened especially near the peak of waves, leading to an early initiation of instability due to cross perturbation.On the other hand, the contours reveal that the vorticity magnitude is decreased compared to the no cross-perturbation case.This could be due to the increased spanwise convection caused by the cross-stream, which also stretches the vortices.The spanwise perturbation also seems to break up the tip of roll-up vortices into multiple edges.Overall, the structures are strongly affected by cross perturbation, leading to the initiation of instability.

Effect of Turbulence
To evaluate the effect of initial turbulence intensity on the interfacial waves, the contours of volume fraction between two gasses are plotted in Figure 8.The turbulence does not affect the initial roll-up (t = 0.5 s); however, the initial development of mixing at t = 0.75 s shows the presence of the inner phase in the mixing is intensified as the turbulence intensity increases.As stated before, the 'mixing layer' is defined as the region near the fluid interface where one fluid is thoroughly interspersed with the other fluid.Consequently, in the mixing layer at t = 1 s, a distinct fraction of the inner gas appears as shown in the boxed regions.Note that the shown boxed regions consist of about 110 × 100 cells at t = 0.75 s and about 50 × 50 cells at t = 1 s, which are significantly resolved compared to the integral length scale, L ≈ 0.014 m.Therefore, the intensified colors of the volume fraction are clearly due to the presence of turbulent structures.By increasing the turbulence intensity from I = 2% to I = 4%, a more distinct presence of the inner flow can be seen in the mixing layer of the two phases.These observations are the motivation for evaluating the power spectral density (PSD) of density fluctuations in the mixing layer, which describes the power of fluctuations as a function of frequency [31].The probe is located along the initial sinusoid, which develops into the center of one of the vortices and eventually positions inside the layer of mixing fluids.As Figure 5a shows, without turbulence, the spectra are limited to low-frequency motions due to extraction and contractions in the mixing layer.However, turbulence expands the range of high-frequency motions significantly due to turbulence structures, as shown in Figure 5a.In addition, it is apparent that the magnitude of the power spectra increases with turbulence intensity.Figure 5b shows the time evolution of the kinetic energy, KE = KE + KE ′ , where KE = 1 2 U 2 and KE ′ = 1/2 u 2 + v 2 + w 2 where u, v, and w are components of velocity fluctuation in the x, y, and z directions, respectively At t ≈ 0.3 s, the kinetic energy values drop to a local minimum because the Kelvin-Helmholtz's billow structure is partly destroyed, as was also observed by Liu et al. [25].Then, KE increases, suggesting the emergence of turbulence.There is a slight increase in the magnitude of KE with an intensity of 2% compared to the no turbulence case.The KE evolution, however, is higher for the turbulence intensity of I = 4%, while it is at a steady rate contrary to the no turbulence and I = 2% cases.This could be due to higher values of the turbulent kinetic energy, KE ′ , associated with the turbulent stage of the flow evolution for I = 4%.

Discussion
The effects of cross perturbation and turbulence intensity in the initial conditions are evaluated for shear-driven two-phase Kelvin-Helmholtz instability to characterize their effect on the development of the mixing layer.The appearance of turbulent structures in the gas stream near the interface influences the shape of the interfacial wave over time.This suggests the combined impact of cross perturbation and variation in the velocity field may influence the interface shape in three dimensions over time, deviating from the two-dimensional theoretical models.This is significant in exposing the need for a more detailed parametric study of 3D effects to properly model shear instability for future highperformance simulations.The results show that cross perturbation affects the instability by increasing the amplitude growth while adding turbulence has almost no effect on the amplitude growth.Linear stability theory explains that the addition of vorticity in the initial flow field should not affect the stability in the inviscid case [4], therefore the results found here agree with the theory.The significance is the change in amplitude growth rate caused by the cross perturbation in the initial interface.If these instability growth rates are affected by the initial conditions in 3D, the accuracy of solutions to real-world problems can be significantly improved if more research is conducted on the trends found in varying the interface contour in the z-direction.
Another observation made is by increasing the turbulence intensity a more distinct presence of the inner flow can be seen in the mixing layer of the two phases.Analysis of PSD reveals that the presence of turbulence expands the range of high-frequency motion significantly due to turbulence structures.The evolution of kinetic energy shows a slight increase as the turbulence intensity is increased.This study can be extended to explore how variations in the initial perturbation modify the interface instability process.As simulations are becoming more prevalent and detailed, it would be helpful to break down the individual phenomena in these flows to have better parametric controls on the initial conditions.Continued study of Kelvin-Helmholtz instability with the theory extended to 3D would be a rich resource for future model development for highly accurate DNSs.

Figure 3 .
Figure 3. Interface contour for various mesh resolutions at t = 0.1 s.

Figure 6 .
Figure 6.Amplitude of interface over time.

Figure 7 .
Figure 7. Q-criterion threshold for cases without turbulence at t = 0.5 s.

Figure 8 .
Figure 8. Contours of volume fraction for various turbulence intensity at different times.

Table 1 .
Grid sizes for KHI convergence.