Numerical investigation of the three-dimensional secondary instabilities in the time-developing compressible mixing layer

TECHNICAL PAPERS

Numerical investigation of the three-dimensional secondary instabilities in the time-developing compressible mixing layer

Ricardo A. Coppola Germanos^I; Leandro Franco de Souza^II; Marcello A. Faraco de Medeiros^III

^[e-mail], University of São Paulo - USP. Escola de Engenharia de São Carlos- EESC 13566-590 São Carlos, SP, Brazil

^IIMember, ABCM, [e-mail], University of São Paulo USP. Instituto de Ciências Matemáticas e de Computação ICMC 13566-590 São Carlos, SP, Brazil

^IIISenior Member, ABCM, [e-mail], University of São Paulo USP. Escola de Engenharia de São Carlos EESC 13566-590 São Carlos, SP, Brazil

ABSTRACT

Mixing layers are present in very different types of physical situations such as atmospheric flows, aerodynamics and combustion. It is, therefore, a well researched subject, but there are aspects that require further studies. Here the instability of two-and three-dimensional perturbations in the compressible mixing layer was investigated by numerical simulations. In the numerical code, the derivatives were discretized using high-order compact finite-difference schemes. A stretching in the normal direction was implemented with both the objective of reducing the sound waves generated by the shear region and improving the resolution near the center. The compact schemes were modified to work with non-uniform grids. Numerical tests started with an analysis of the growth rate in the linear regime to verify the code implementation. Tests were also performed in the non-linear regime and it was possible to reproduce the vortex roll-up and pairing, both in two-and three-dimensional situations. Amplification rate analysis was also performed for the secondary instability of this flow. It was found that, for essentially incompressible flow, maximum growth rates occurred for a spanwise wavelength of approximately 2/3 of the streamwise spacing of the vortices. The result demonstrated the applicability of the theory developed by Pierrehumbet and Widnall. Compressibility effects were then considered and the maximum growth rates obtained for relatively high Mach numbers (typically under 0.8) were also presented.

Keywords: mixing layer, linear stability theory, compact finite-difference scheme, stretching grid, linear and non-linear growth rate

Introduction

The mixing layer, a classical phenomenon of hydrodynamic instability, comes in numerous practical situations such as aerodynamics, propulsion and environmental engineering. As a specific example, progress in space research is dependent on vehicles capable of carrying higher payloads into orbit and therefore requires the development of more efficient propulsion systems. Such systems involve combustion and turbulent mixing. In fact, a detailed understanding of the physics of mixing layers is essential for the development of new turbulence and mixing models. Compressibility is a key factor in many scenarios. Owing to these aspects, this field has constituted one of the main research themes in turbulence over the last thirty years. Nevertheless, there are aspects of the problem that need further investigation. The numerical investigations here presented are devoted to the regime governed by non-linearity.

The engineering research and modern design requirements pose great challenges in computer simulation to engineers and scientists who are called upon to analyze phenomena in continuum mechanics. However, the advent of the computational fluid dynamics technology has had a profound impact on fluid dynamics research. It has enabled careful and detailed investigations that were inconceivable a few decades ago. Nowadays, direct numerical simulations are generally employed to solve the Navier-Stokes equations and to study the physics of hydrodynamic instability.

The mixing-layer instability, often called Kelvin-Helmholtz instability, promotes the formation of spanwise vortex structures. Analytical study on this problem began at the end of the nineteenth century, when Rayleigh (1980) demonstrated that the profile has to have an inflexion point to be inviscidly unstable. Further analytical work also revealed the preference for amplification of two-dimensional disturbances (spanwise wavenumber equals zero) in incompressible flows (Squire, 1933). Numerical solution of the equations of Linear Stability Theory (LST) to the incompressible mixing layer was presented only much later, by Michalke (1964). The effect of compressibility on the mixing layer instability was investigated in laboratory experiments by Lessen, Fox and Zien (1965, 1966) and Gropengiesser (1970). The experimental work published by Birch and Eggers (1973) showed a reduction in the growth rate of the two-dimensional disturbances as the Mach number increased. Brown and Roshko (1974) found that the effect of density ratio alone could not be responsible for this reduction implying that the compressibility effect is the mechanism that leads to it.

Renewed interest in compressible mixing layers motivated experiments by Papamoschou and Roshko (1988) in which high-speed mixing layers with various velocity and density ratios were investigated. In this work, a parameter called convective Mach number is proposed:

where U₁ and U₂ are the free-stream velocities and c₁ and c₂ are the free-stream sound speeds. The subscripts 1 and 2 refer to the velocity at upper (y > 0) and lower (y < 0) free stream, respectively. Works involving the study of linear stability were presented by Sandham and Reynolds (1990) and Fortuné (2000). When the disturbances become large the non-linear effects must be taken into account. The theory associated with these effects is more complicated than its linear counterpart. Reviews about the theory can be found in Bayly (1988). Both theory and experiments show that these flow disturbances grow and saturate in a limit cycle pattern of co-rotating vortices. In turn, the vortices are themselves unstable to subharmonic disturbances, meaning that, if a subharmonic oscillation exists in the flow it will grow. It eventually results in a pairing of the vortices created by the fundamental disturbance. This phenomenon is an example of secondary instability.

Investigation of three-dimensional disturbances (spanwise wavenumber different from zero) triggering secondary instability was also carried out. In this scenario, large spanwise vortices, called rollers, represent the saturated state reached by the mixing layer instability. A major work in this field was produced by Pierrehumbert and Widnall (1982). They assumed a base flow consisting of a hyperbolic tangent profile with superposed Stuart vortices (Stuart, 1967). These are steady solutions of the incompressible Navier-Stokes equations. Two classes of instability were found: fundamental and subharmonic. In the fundamental case, two modes were found corresponding to vortex core deformations called respectively bulging and translative modes. The phase was the parameter that selected the type of mode. The translative mode studied by Pierrehumbert and Widnall (1982), and Sandham and Reynolds (1990) was the more unstable. In Sandham and Reynolds (1990) the wavelength used was roughly equal to the spacing between the streamwise vortices, as found in the experiments of Bernal and Roshko (1986). Pierrehumbert and Widnall (1982) also showed that for the fundamental mode maximum growth rate was obtained for a spanwise wavelength approximately 2/3 of the streamwise spacing of the vortices.

Pierrehumbert and Widnall (1982) considered that the three-dimensional fundamental mode was responsible for the initial three-dimensional deformation of the spanwise vortices leading to the generation of longitudinal vortex structures called braids. The braids appeared only at non-linear stages of the three-dimensional development of the rollers and as such were not directly accessible by the primary linear stability theory. For the three-dimensional subharmonic mode, Pierrehumbert and Widnall (1982) found another kind of secondary instability called helical pairing. However, the most amplified subharmonic waves were the two-dimensional ones. This corresponds to the two-dimensional pairing process most commonly observed in experiments. Ragab (1989) presents theoretical studies on secondary stability for compressible mixing layer. In his work, the study of the subharmonic instability found that for convective Mach numbers of about 0.4, the helical pairing mode was more amplified than the two-dimensional pairing mode. This is a further indication that the dominant character at high Mc is three-dimensional.

The main contribution of the current work consists in verifying the growth rate of the three-dimensional flow disturbances in the secondary instability. This was done for both compressible and incompressible flows, and the results were analyzed in the light of the theory developed by Pierrehumbert and Widnall (1982). Based on the works presented in this review, the current code was subjected to numerous tests to assure the correct implementation of the numerical methodology. These tests mainly involved the reproduction of Linear Stability Theory results. The initial sections give the main characteristics of both the code and the simulations. Then the results of tests for the time developing mixing layer are presented, first in the linear regime and then in the nonlinear regime. The conclusions are given in the last section. The introduction should contain information intended for all readers of the journal, not just specialists in its area. It should describe the problem statement, its relevance, significant results and conclusions from prior work and objectives of the present work.

Nomenclature

xi = Cartesian coordinate in (x,y,z) direction, m t = time, s c = sound speed, m/s M = Mach number, dimensionless M_c = convective Mach number, dimensionless Re = Reynolds number, dimensionless Pr = Prandtl number, dimensionless U = free stream velocity, m/s T = free stream temperature, K u_i = velocity component in (x,y,z) direction, m/s p = pressure, N/m ² E_t = total energy, J/(Kg K) e = internal energy, J/(Kg K) q_i = heat transfer heat, J/(Kg K) c_v = specific heat at const. volume, J/(Kg K) c_p = specific heat at const. pressure, J/(Kg K) k = thermal conductivity, J/(Kg K) A = amplitude of disturbance, m f = generic function f' = derivative of the function h = grid spacing N = number of grid points dx = grid spacing in the streamwise direction, m dy = grid spacing in the normalwise direction, m dz = grid spacing in the spanwise direction, m L_x = wavelength in the streamwise direction, m L_y = wavelength in the normalwise direction, m L_z = wavelength in the spanwise direction, m y_c = center of computational domain in the y-direction, m H = length of computational domain in the y-direction, m LHS = left hand side RHS = right hand side a,b,c = coefficients used in the numerical methods Greek Symbols µ = dynamic viscosity, kg/(m s) ρ = density, kg/m ³ γ = ratio of specific heats, dimensionless τ = viscous stress tensor δ_Ω0 = vorticity thickness of the initial profile α = wavenumber of the disturbance in the x-direction, 1/m β_n = wavenumber of the disturbance in the z-direction where n corresponds to the index of spanwise wavenumber, 1/m ω = frequency of the disturbance flow, 1/m σ = decay rate of the disturbance in the normal direction φ = phase between a two-dimensional and an oblique wave η = transformation of y coordinate κ, λ = coefficients used in the numerical methods Superscripts * = dimensionless variables - = base flow variables ´ = flow disturbance variables ^ = filtered function Subscripts mim = minimum value of any variable max = maximum value of any variable 1 = refers to the upper (y > 0) free stream or dominant mode 2 = refers to the lower (y < 0) free stream or subharmonic mode

Submission

In the current study, the governing equations are the compressible Navier-Stokes equations. The continuity equation is

The momentum equations for the velocity components in the streamwise,normal and spanwise directions can be written as

and the energy equation is

where x_i is the Cartesian coordinates (x,y,z), t is the time, ui is the velocity components (u,v,w), ρ is the density and p is the pressure. The total energy is given by E_T = ρ(e + u² + v² + w² /2) and the primitive variable e is the internal energy. The non-dimensional constitutive relations for a Newtonian fluid and Fourier heat conduction are

And

where the Prandtl number is defined as P_r =c_pµ/k , γ is the ratio of specific heat, c_p is the specific heat at constant pressure and k is the thermal conductivity. The Reynolds number of the flow is defined as

The dimensional variables (ρ₁,U₁, µ ) are the density, the velocity and the dynamic viscosity of the base flow at the upper (y > 0) free stream. The parameter δ_Ω0 d is the vorticity thickness of the initial velocity profile given by

where the subscripts 1 and 2 refer to the velocity at the upper (y > 0) and lower (y < 0) free stream, respectively. The perfect-gas law for non-dimensional pressure and temperature is

And

These conservation equations were presented with the following non-dimensionalization scheme:

where α is the streamwise wavenumber and ω is the radial frequency of the imposed disturbance. The superscript ^* indicates dimensional parameters.

In problems of hydrodynamic instability the variables are often decomposed into two parts: the base flow and a small disturbance flow. This decomposition, adopted in the present work, can be written in the following way

where ( ^_ ) indicates the base flow and ( ^' ) refers to the disturbance flow. In the current analysis, the base flow is invariant in the streamwise direction and the components (v,w) of the base velocity are null. Sometimes this decomposition is used in the code implementation. However, in the current implementation the total variables were used.

Methodology

Base Flow

Many different velocity profiles have been proposed in the literature to model a mixing layer. The profile used in the current numerical work corresponds to a hyperbolic tangent function (Fortuné, 2000), in which the upper part travels to the right and the lower to the left. The base velocity profile was defined as

For this profile, the vorticity thickness is given by Eq. (8). This profile has often been used, partially because it is analytical. Consequently, the flow speed can be calculated exactly for each value of y. The derivatives of this function are known as well. A disadvantage of this profile is that it satisfies the momentum and energy equations only approximately. Therefore, the simulations exhibit a transient while the profile relaxes to the exact solution. As this transient is very short, the use of Eq. (13) does not affect the results at the later non-linear stages. Indeed, initial base flows that satisfy the conservation equation only approximately are commonly used (Sandham and Reynolds, 1990).

A uniform pressure was assumed for the initial base flow. The initial mean temperature profile was calculated from a solution of the compressible boundary layer energy equation assuming unity Prandtl number (White, 1974). For the antisymmetric mean velocity profile considered here and with equal free stream temperatures, the general relation of Crocco-Busemann leads to

This expression also involves the assumption that the flow is parallel. With a constant density along the velocity profile, the convective Mach number (M_c) of this flow is equal to the free stream Mach number (M).

Initial Condition

Along with the base flow profiles, the flow disturbances have to be defined. Numerical errors engender perturbations which are sufficient to trigger the instabilities in the mixing layer, but for very accurate codes such as the current one, this process takes a very long time. The development of the instability and the formation of vortex structures can be largely accelerated and better controlled by the introduction of flow disturbances.

In this work, the flow disturbances were exponential functions in the y-direction. These functions are convenient because they decrease really quickly and satisfy the boundary conditions of no disturbances at y ±∞ . In the x-direction the disturbance was periodic. All these issues were used to define v′ . For two-dimensional disturbances, u′ was calculated from the continuity equation for incompressible flow. This is only true for very small Mach numbers, but it did not affect the long time solution. The flow disturbances could then be written as

where A is the amplitude, α is the streamwise wavenumber and σ is the decay rate of the disturbance in the y-direction. The subscripts 1 and 2 refer to the fundamental and subharmonic modes respectively.

Based on similar hypotheses, the three-dimensional initial flow disturbances used were

where β_n represents the spanwise wavenumber, the lowest being for n = 0 , and φ is the phase difference between the two-dimensional wave and the pair of three-dimensional, oblique waves.

Boundary Condition

Periodic boundary conditions were used in the stream-and span-wise directions. In the normal direction, an unbounded boundary condition would be necessary. Thus boundary conditions were adapted to simulate an infinite domain, even though the computational domain was finite. In the current simulations, the normal component of velocity in the free stream was set to zero. This corresponds to the impermeability condition. For the other primitive variables the first derivatives in the normal direction were set to zero. This boundary condition is often referred to as free-slip. For a sufficiently large distance from the mixing layer, this boundary condition produces accurate results.

Numerical Method

Simulations performed in the current work required high accuracy of the numerical solution (Lele, 1992; Souza, 2003). Compact finite-difference schemes were then adopted. The compact schemes for spatial derivatives are extremely attractive when explicit time advancement schemes are used. The most popular compact scheme for spatial derivatives, also called Padé scheme, is the symmetric sixth order version. In the current implementation this approximation was used only for the interior points. This scheme is symmetric and does not exhibit dissipative errors.

For compact schemes, the finite-difference approximations to the derivatives of the function are expressed as a linear combination of the given function values and derivatives on a set of nodes. First, a uniformly spaced mesh was considered where the nodes are indexed by i, which varies from 1 to N. In these schemes, the value of f′ is dependent on all the others nodal values. In general, compact (also called implicit) schemes are significantly more accurate for short length scales than non-compact schemes (Collatz, 1966; Kopal, 1961). Lele (1992) emphasizes the importance of using methods of high-order and shows schemes with 6^th order approximation for the interior of the mesh.

Here, the following 6^th order compact (implicit) derivatives were used

In the equation above, h is the grid spacing. This equation corresponds to the approximation used in the streamwise and spanwise directions, for all points.

The compact approximations for the second derivative are similar to that of the first. These derivatives that correspond to the viscous terms in the governing equations may involve the evaluation of successive first derivatives. When a spectral method is used, there is no loss of accuracy if these derivatives are computed in such a way. However, with finite-difference methods one finds that two applications of a first derivative give a significantly worse representation at high wavenumbers than a direct second derivative computation. The spectral analysis carried out by Lele (1992) shows that for the first derivative the modified wavenumber at high wavenumbers goes to zero. The solution for this problem is to calculate the second derivatives directly, avoiding these problems and providing more accuracy. This approach, however, is non-conservative, but for the case studied here, it did not bring any difficulty. The scheme for second derivatives adopted in this work was

The approximation above was used for the stream and span-wise directions, which have periodic boundary conditions.

In the normal direction, a 2^nd order approximation at the boundaries was used. Tests were performed with this and other schemes. The results showed that the numerical scheme is more stable with 2^nd order approximations than the high order schemes. It is emphasized that, for a sufficiently long distance from the shear zone, it is possible to use low order without reducing the overall accuracy of the simulations.

For the interior points in the y-direction, the same compact scheme described above was used for first and second derivatives, Eq. (17) and Eq. (18) respectively. At the boundary i = 1, the scheme for first derivative can be written as

For the points next to the boundary, i = 2

The second derivatives, at the boundary i = 1, were discretized as

For the points next to the boundary, i = 2

For the points at the opposite boundary, i = N and i = N-1, similar approximations were used.

Table 1 shows the stencil size and the truncation error for first and second derivatives.

The time-advancement of the computational variables ( ρ ,ρu_i ,E ) was obtained by a 4^th order Runge-Kutta method. The scheme used here works in four steps (Ferziger and Peric, 1997). The combination of these steps results in a order accurate Where algorithm in time.

The non-linear terms in the Navier-Stokes equation can produce aliasing errors. In order to remove this error, a high-order compact filter was implemented (Lele, 1992). The numerical filter was applied after the last step of the Runge-Kutta scheme. This filter consists in recalculating the distribution of the primitive variables and is of 4^th order accuracy as below

Implementation of the filtering schemes on domains with nonperiodic boundaries requires the near boundary nodes to be treated separately. Therefore, near boundary explicit formulas were needed:

The truncation error for these methods is shown in Tab. 2

Grid Stretching

The governing equations can be transformed from a Cartesian coordinate system to any general orthogonal coordinate system. Hence a non-uniformly spaced computational grid in the physical plane can be transformed into uniformly spaced grid in the computational plane.

Simple transformations can be used to cluster grid points in regions of large gradients, where more resolution is required. Here, in the points near the free stream the problem permits less resolution, while at the center of the domain greater resolution is needed. Therefore, it is possible to stretch the grid and decrease the overall number of points used while maintaining the overall accuracy.

In this work, a formula from Anderson, Tannehill and Pletcher (1984) giving a constant stretching was adopted. τ is a stretching parameter that varies from zero (no stretching) to large values. It produces the most refinement near y = yc, which is located at the center of the grid

Where

In Eq. (27), the parameter H is the size of the computational domain. In order to apply this transformation to the governing equations, partial derivatives have to be taken. For the first derivatives we have

where

This relation should be used in combination with the compact schemes. Applying relation (28) to the derivative approximations, the following tridiagonal compact scheme for the first derivatives was obtained

where The parameters κ , a e b can be defined according to the compact finite-difference scheme used. A similar procedure can be applied for the second derivatives

where

Equation (31) can be rewritten as

Using the relation and replacing the term between square brackets in Eq. (33), one has

Therefore, applying relation (34) to the second derivative approximation

Bogey et al. (2000) recommend that the grid stretching should not exceed 1.8% in order to avoid problems with space derivatives. One basic assumption of these methods is that the mesh must be sufficiently smooth so that ∂/∂η and ∂²/ ∂η ² can be calculated without appreciable loss in the overall accuracy.

Computational Tests

In this section, results from direct numerical simulation of the compressible Navier-Stokes equations are used to test and verify the current code. Simulations of two-and three-dimensional flow disturbances in linear and non-linear regime were performed and the results were compared with existing theories.

Linear Stability

This section presents results of the evolution of two-dimensional sinusoidal disturbances in a mixing layer in the regime governed by the linear theory. Simulations were performed and the growth rate for uniform and non-uniform grids were compared to theoretical ones in order to analyze the code implementation and efficiency of the grid stretching. Theoretical results were extracted from Sandham and Reynolds (1990).

An important aspect to be considered in these simulations is the treatment of the vertical diffusion. This diffusion increases the width of the mixing layer during the simulation, which implies a variation of the mean flow over time. The strategy adopted here to avoid this diffusion was to cancel the vertical diffusion terms for the base flow

= 0.

First the mixing layer problem was simulated with the Euler equations and two tests were performed. These tests were carried out using a disturbance with only one mode, that is, only one streamwise wavenumber. The wavenumbers of the disturbances were close to the values of maximum amplification as predicted by linear theory. The mesh for both uniform and non-uniform grids had 40 x 80 points in the x-and y-directions respectively. The grid spacing used in the streamwise direction with a uniform mesh was dx ≈ 0.38 for M_c = 0.4 and dx ≈ 0.60 for M_c = 0.8. For the uniform mesh, in the normalwise direction the grid spacing utilized was dy ≈ 0.30. For the non-uniform mesh in the y-direction, the minimum and maximum mesh spacing were dy_min ≈ 0.0089 and dy_max ≈ 1.67. The initial amplitude of the disturbance given by Eq. (15) was approximately 10^-6. This amplitude ensured that the phenomenon started in the regime governed by the linear theory. The non-dimensional time step (dt) of these simulations was 10^-3 .

Figure 1 shows the growth rate of the unstable wave as a function of non-dimensional time for M_c = 0.4. The vertical coordinate is in logarithm scale. The dashed line shows the result for a non-uniform grid with a stretching parameter τ= 12 . The growth rate from this case was about 0.30 which is very close to the theoretical results (solid line) of approximately 0.31. The dashed-dotted line is the numerical result for a uniform grid. The amplification rate obtained was about 0.27. The disagreement might be attributed to the low resolution in the interior of the computational domain generated by poor accuracy of the uniform grid spacing simulation in the rotational region.

Figure 2 presents results for M_c = 0.8. Similar to the previous results, the dashed-dotted line gives at low accurate amplification rate of about 0.11 for uniform grid. The theoretical results for this convective Mach number give an amplification rate of about 0.14. The solid line represents the simulation with the use of grid stretching. The growth rate obtained in this simulation was 0.125. This result is closer to the theoretical prediction.

Simulations of the Navier-Stokes equations rather than the Euler equations were also performed. The grid and the initial amplitude of the disturbance adopted were the same as that in the previous inviscid simulations. The Reynolds number was 500. Again the wavenumbers selected corresponded to maximum amplification according to the theory.

In Fig. 3 the dashed-dotted line shows the time evolution of a two-dimensional disturbance at M_c = 0.4. For a grid with constant spacing the amplification rate was about 0.26. This is an underestimation compared to the theoretical results that give an amplification rate of about 0.28 (solid line). The dashed line shows a growth rate of approximately 0.27 with the use of grid stretching. Although this result still underestimates the theoretical values, it is closer to the linear analysis.

Figure 4 shows the time evolution for M_c = 0.8. The same analysis was made. The growth rate obtained for the uniform mesh was approximately 0.07, while for a mesh with grid stretching it was about of 0.09. Both simulations were close to the analytical results which give a growth rate about 0.11, but once more the grid stretching improved the result.

According to the results presented above, the stretching used in the y coordinate significantly improved the accuracy by clustering the points in the region of interest. This improvement is extremely necessary for the three-dimensional simulations that require a much larger number of points in the mesh. In the next section the code verification was extended to the non-linear regime.

Two-Dimensional Non-Linear Stability

The full compressible Navier-Stokes equations were used to carry out analysis in the non-linear regime. Again the evolution of a small two-dimensional disturbance in the mixing layer was simulated. The previous simulations confirmed that the judicious use of non-uniform grids provides better results and, consequently, the simulations here were performed with this numerical technique. The idea was to reproduce some classical phenomena in the nonlinear regime. These comparisons were mainly qualitative, but provided confidence that the code was correct.

The theory associated with the non-linear effects is more complex than the linear theory. Reviews of the subject can be found in Bayly (1988). The mesh for these problems had a dimension of 60 x 100 in the x-and y-directions respectively. The computational domain in the streamwise direction was L_x = 4π/0.82 ≈ 15. In the normal direction the computational domain was L_y = 24. The grid spacing used in the streamwise direction with a uniform mesh was dx = 0.25. For the nonuniform mesh in the y-direction the minimum and maximum mesh spacing were dy_min = 0.07 and dy_max = 0.7. The initial amplitude of the disturbance was approximately 10^-3. This amplitude ensured that the phenomenon started in the regime governed by the linear theory. The time step dt of these simulations was 10^-2 . The disturbance was composed of only one mode. Both theory and experiments show that in a mixing layer the disturbance does not grow to infinite. Instead, these disturbances saturate in a limit cycle pattern of co-rotating vortices.

Figure 5 shows the evolution of the two-dimensional disturbances in time. Initially, in the linear regime, the disturbances are very small and display a sinusoidal pattern. Next, the fundamental mode grows and saturated vortices are formed. After the vortices reach the limit cycle oscillation, they are dissipated due to viscous effects.

Figure 6 shows the evolution of the mixing layer with the introduction of both a fundamental and a subharmonic mode. In the initial stage, the growth of a fundamental wave up to vortex saturation was observed. In turn, the vortices are themselves unstable to a subharmonic disturbance. The result is a pairing of vortices. It is important to emphasize that in the system there is no mechanism for the production of subharmonic disturbances, but only for amplification. After the fundamental mode saturates, the subharmonic mode grows and two of the primary structures begin to rotate around each other. Apparently the vortices in each picture were perfectly identical, which testify to the accuracy of the code. After that, the pairing occurs between these two vortices and one large vortex results.

Results for M_c = 0.6 are shown in Fig. 7. The first frame shows the flow structure at a non-dimensional time equals to 55. Here, the two-dimensional disturbances excited by the fundamental mode are not evident. It can be observed in the sequence for Mc = 0.4 that at the same time the flow already presented a sinusoidal pattern. In other words, the disturbance at M_c = 0.4 grew faster than the disturbance at M_c = 0.6. Proceeding with the analysis for M_c = 0.6, it can be observed in the second frame that the fundamental mode was visible for a non-dimensional time equals to 73. The next frame, at time about 85, presents the structures of co-rotating vortices before a merging. The results for M_c = 0.4 show a similar phenomenon at an earlier time of 65. The disturbance excited by the subharmonic perturbation then becomes evident, leading to the pairing which occurs at time around 117. For M_c = 0.4 the same phenomenon occurred at a time of about 100. In other words, the pairing for M_c = 0.6 took a longer time to occur.

The fact that the growth rates are not the same makes the comparison of the results of Figures 6 and 7 more difficult. Possibly, another interesting way of analyzing these results is to use a normalized time that would take into account the difference in growth rates. This can be defined as the product ω_it . Since in both simulations the initial disturbance amplitude was identical, they would also have identical amplitudes throughout the linear regime for identical values of ω_it . The values of ω_it for the frames shown were calculated and are also given in the figures. In this normalized time frame the figures suggest that the vortex pairing in the process is comparatively faster for the higher Mach number case. This may be related to the fact that the pairing has a dynamics that is more related to the vortices configuration than to the linear regime growth rates.

The tests carried out for two-dimensional flow disturbances in the non-linear regime showed that the current code can recover the classical secondary instability phenomenon of the time-developing mixing layer.

Three-Dimensional Non-Linear Stability

The analysis was extended to the three-dimensional simulations. Once again the code was tested and the effect of convective Mach number in the non-linear evolution of various combinations of unstable waves was considered. In particular, the temporal development of a combination of a two-dimensional and a pair of oblique waves was performed in all simulations reported in this section.

The secondary instability process studied by Pierrehumbert and Widnall (1982) presents two different modes known as bulging and translative modes. The effect of the phase φ in Eq. (16) was here considered, first with the objective of recovering the secondary instabilities phenomenon. In this analysis the initial step of the flow disturbance presented a sinusoidal behavior. The initial amplitude was A₁ = 10^-3 for the two-dimensional wave and A₂ = 10^-4 for the pair of oblique waves. Two simulations were run: one for phase φ= 0 to recover the development of the bulging mode and another with φ =π /2 , for the translative mode. The Reynolds number was 500. The mesh adopted in these simulations was 60 x 100 x 60 that corresponded to the number of points in the (x,y,z) directions. The computational domain in the streamwise and spanwise direction was L_x = L_z = 4π / 0.82 ≈ 15. In the normal direction the computational domain was L_y = 24. The grid spacing used in the streamwise and spanwise direction with a uniform mesh was dx = dz = 0.25. The non-uniform grid in the y-direction ranged from dy_min = 0.07 to dy_max = 0.7. Pressure is a good way to identify large-scale vortical structure in the flow and, therefore, all figures show pressure isosurfaces. The streamwise wavenumber for all simulations was close to that of the largest growth rate from linear stability analysis.

Figure 8 presents an iso-surface of pressure for Mc = 0.4 showing the development of the bulging mode. In this mode the core of the spanwise vortices has a diameter that varies sinusoidal in the spanwise direction. The structure developed from the translative mode of the secondary instability is shown in Fig. 9. This mode was more amplified and gave a structure composed of a vortex core that oscillates in amplitude and locations in the spanwise direction.

Tests were further extended to simulate the three-dimensional subharmonic secondary instability (Fig. 10). Because of the three-dimensional subharmonic mode, the vortex cores are shifted alternately above and below the plane y = 0. One can observe the generation of spanwise phase dislocations via localized pairing and a coherent three-dimensional structure. The instability is helical in the sense that it causes neighboring vortex tubes to twist around one another. The helical pairing has been experimentally confirmed by Chandrsuda et al. (1978).

Figure 11 shows results of a simulation for M_c = 0.8 with the evolution of a translative mode. One can found a weakening of the spanwise structure, which develops larger amplitude of spanwise oscillation. One can also observe the development of oblique vortices in the region between two spanwise rollers, where at low convective Mach number the streamwise braid vortices were formed. Therefore, a change in the developed large-scale structure was observed as the convective Mach number was increased, with vortical regions oriented in a more oblique manner at the higher convective Mach numbers. Based on these results, one can conclude that the results of the three-dimensional simulation corresponded to the results obtained from the theory.

Numerical Results

The main goal of the current work was to verify which three-dimensional mode was the most amplified in the flow-field after the two-dimensional structure of vortices saturated. In addition, an analysis of the compressibility effects on the amplification rate of the oblique waves was also carried out. As discussed in the literature, the development of three-dimensionality in a nominally two dimensional mixing layer is often attributed to secondary instability effects as proposed by Pierrehumbert and Widnall (1982). The model proposed by them was concerned with the saturated state of the mixing layer, when the fundamental disturbances have developed into spanwise rollers. Therefore, in order the check the model, the excitation procedure had to be considered with care. The approach used here was as follows. First, a two dimensional wave was excited and allowed to grow up to saturation. Only at this stage the flow was excited with three-dimensional seeds. The procedure is illustrated in Fig. 12.

The lowest spanwise wavenumber was defined to be β₁/α= 0.25 . Within the nomenclature assumed in this paper, α represents the streamwise wavenumber of the two-dimensional wave introduced to trigger the fundamental instability. In the test presented here, 12 oblique modes were used in total. It means that the highest oblique mode β₁₂/α was equal to 3.0. Moreover, the oblique mode β₆/α= 3/ 2 = 1.5 corresponds to that of maximum amplification according to the theory given by Pierrehumbert and Widnall (1982). The main goal here was to verify this theoretical prediction.

The initial amplitude for all waves was A₁ = 10^-3. The Reynolds number selected was 1000. The mesh adopted in the simulations was 15 x 65 x 180 that corresponds to number of points in the (x,y,z) directions. The grid spacing used in the stream-and span-wise direction which had a uniform mesh was dx = 0.54 and dz = 0.17 respectively. The non-uniform grid in the y-direction ranged from dy_min = 0.19 to dy_max = 1.85. The use of grid stretching in the current simulations was essential due the large matrix utilized to solve the three-dimensional equations. The stretching provided a good resolution at the center of the computational domain.

Figure 13 presents the growth rate at M_c = 0.1 for various spanwise wavenumbers ( β_n ). The streamwise wavenumber of the disturbances were close to the values of maximum amplification predicted by linear theory for the respective Mach. In this figure, it can be seen that the most amplified wavenumber corresponds to β /α= 1.5. This result confirmed the theory given by Pierrehumbert and Widnall (1982). The results in Fig. 14 showed that at M_c = 0.4 the most amplified mode corresponds to β /α= 1.

Summary and Conclusions

In this work the numerical simulation of a time-developing mixing layer was performed. The governing equations were the compressible Navier-Stokes equations. A 6^th order compact finite-difference scheme was used for discretizing the spatial derivatives. The scheme adopted was also time accurate, using a 4^th order Runge-Kutta scheme. In order to remove short length scales, a 4^th order compact filter was applied. Periodic boundary condition was implemented in the x-and z-directions. A free-slip boundary condition was used in the normal direction. The simulations were run with a grid stretching in the y coordinate. This technique was utilized both to improve code performance and to remove sound waves produced by both the formation of vortical structures and the pairing.

First computational tests were performed to verify the code. In the linear regime, it was possible to obtain growth rates closer to the theoretical results with the use of stretching in the y coordinate for the same number of points. The testes covered a number of M_c. After that, tests involving the two-dimensional secondary instability were also performed. In this case, a sub-harmonic disturbance was excited to reproduce the vortex pairing. Simulations of the flow with two-dimensional waves and a pair of oblique waves were also run. The outcome reproduced the behavior expected according to the study carried out by Pierrehumbert and Widnall (1982). The main idea was to recover the classical modes of the secondary stability governed by the phase difference between a two-dimensional wave and a pair of oblique waves, namely, the bulging and the translative mode. In the nonlinear regime, the comparison with theoretical results was done on a qualitative basis. Nevertheless, simulations of the time evolution obtained for two-and three-dimensional instabilities compared well with others works presented in the literature (Pierrehumbert and Widnall, 1982; Sandham and Reynolds, 1989, 1990; Fortuné, 2000).

The contribution of the current work consisted in studying the growth rate of the flow disturbances in the three dimensional secondary instability of the mixing layer. The purpose was to verify which oblique wave corresponded to the most amplified and how the Mach number affected the problem. In the simulations, it was important that the disturbances that trigged the secondary instability were introduced after the fundamental two-dimensional disturbances reached a saturated level. Initially, the simulations were run for low Mach number (typically under 0.1). In this regime the numerical simulations were compared to the theoretical results extracted from the study carried out by Pierrehumbert and Widnall (1982). Based on their theory the most amplified spanwise wavenumber corresponds to 3/2 of the streamwise one. The numerical results obtained here confirm in a more systematic way the validity of this theory. Simulations were extended to investigate the compressibility effects on the amplification rates. For M_c = 0.4 it was found that the most amplified three-dimensional wave corresponded to β /α=1. Simulations for M_c = 0.8 showed the same trend. It appears that in the secondary instability regime, this flow may present a more two-dimensional character as the Mach number increases, inasmuch as the preferred spanwise wavenumber reduced as the Mach increased. This is interesting and contrary to the primary instability, for which the higher spanwise wavenumbers dominate as the Mach number increases.

Acknowledgments

This work was funded by the Air Force Office of Scientific Research (AFOSR/USA) under grant number FA9550-07-1-0055 and by the State of São Paulo Research Support Foundation (FAPESP/Brazil).

Paper accepted March, 2009.

Technical Editor: Francisco R. Cunha

• Anderson, D.A., Tannehill, J.C., and Pletcher, R.H., 1984, "Computational Fluids Mechanics and Heat Transfer", Hemisphere Publishing Corporation, pp. 1-591.
• Bayly, B.J., Orszag, S.A. and Herbert, T., 1988, "Instability mechanisms in shear-flow transition", Anne. Rev. Fluid Mechanics, Vol. 20, pp. 487-526.
• Bernal, L. and Roshko, A., 1986, "Streamwise vortex structure in plane mixing", Journal of Fluid Mecanics, Vol. 170, pp. 429-525.
• Birch, S.F. and Eggers J.M., 1973, "A Critical Review of the Experimental Data for Developed Free Turbulent Shear Layers", NASA SP. 321, pp. 11-40.
• Blumen, W., Drazin, P.G. and Billings, D.F., 1975, "Shear layer instability of an inviscid compressible fluid. Part 2", Journal of Fluid Mecanics, Vol. 71, pp. 305-316.
• Bogey, C., Bailly, C. and Juve, D., 2000, "Numerical simulation of sound generated by vortex pairing in a mixing layer", AIAA Journal, Vol. 38, pp. 2210-2218.
• Brown, G.L. and Roshko, A., 1974, "On density effects and large structure in turbulent mixing layer", Journal of Fluid Mecanics, Vol. 64, pp. 775-816.
• Canuto, C., Hussaini, M.Y., Quarteroni, A. and Zang, T.A., 1987, "Spectral Methods in Fluid Dynamics", Springer-Verlag, New York, pp. 1-557.
• Chandrsuda, C., Mehta, R.D., Weir, A.D. and Bradshaw, P., 1978, "Effect of free-stream turbulence on large structure in turbulent mixing layers", Journal of Fluid Mecanics, Vol. 85, pp. 693-.
• Collatz, L., 1966, "The Numerical Treatment of Differential Equations", Springer-Verlag, New York.
• Ferziger, J. H. and Peric, M., 1997, "Computational Methods for Fluid Dynamics", Springer-Verlag, New York.
• Fortuné, V., 2000, "Étude par Simulation Numérique Directe du Rayonnement Acoustique de Couches de Mélangue Isothermes et Anisothermes", PhD. thesis, Université de Poitier.
• Gropengiesser, H., 1970, "Study on the Stability of Boundary Layers and Compressible Fluids", NASA TT, F-12786.
• Harris, P.J. and Fasel, H.F., 1997, "Numerical Investigation of Transitional Compressible Plane Wakes", PhD. Thesis Submitted to the Faculty of the Department of Aerospace and Mechanical Engineering, University of Arizona.
• Kopal, Z., 1961, "Numerical Analysis", Wiley, New York.
• Lele, S.K., 1992, "Compact finite difference schemes with spectral-like resolution", Journal of Computational Physics, Vol. 103, pp. 16-42.
• Lesieur, M., 1997, "Turbulence in Fluids, (3^rd Edition)", Ed. Kluwer Academic Publisher.
• Lessen, M., Fox, J.A. and Zien, H.N., 1965, "On the inviscid of the laminar mixing of two parallel streams of a compressible fluid", Journal of Fluid Mecanics, Vol. 23, Part 2, pp. 355-367.
• Lessen, M., Fox, J.A. and Zien, H.N., 1965, "Stability of the laminar mixing of two parallel streams with respect to supersonic disturbances", Journal of Fluid Mecanics, Vol. 25, pp. 737-742.
• Lin, S.J. and Corcos, G.M., 1984, "The mixing layer: deterministic models of a turbulence flow. Part 3. The effect of a plane strain on the dynamics of streamwise vortices", Journal of Fluid Mecanics, Vol. 141, pp. 139-178.
• Metcalfe, R., Orszag, S., Brachet, M., Menon, S. and Riley, J., 1987, "Secondary instability of a temporally growing mixing layer", Journal of Fluid Mecanics, Vol.184, pp. 207-243.
• Michalke, A., 1964, "On the Instability of the Hyperbolic-Tangent Velocity Profile", DVL-Institut für Turbulenzforschung, Berlin, pp. 543-556.
• Papamoschou, D. and Roshko, A., 1988, "The Compressible Turbulent Shear Layer: An Experimental Study", J. Fluid Mech, Vol. 197, pp. 453-477.
• Pierrehumbert, R.T. and Widnall, S.E., 1982, "The two-and three-dimensional Instabilities of a spatially periodic shear layer", Journal of Fluid Mecanics, Vol. 114, pp. 59-82.
• Ragab, S.A. and Wu, J.L., 1989, "Linear subharmonic instabilities of periodic compressible mixing layers", AIAA Paper, No. 89-0039.
• Rayleigh, L., 1880, "On the stability, or instability, of certain fluid motions", Proceedings of London Mathematics Society, Vol. 11, pp. 57-70.
• Sandham, N.D. and Reynolds, W.C., 1989, "Compressible mixing layer: linear theory and direct simulation", AIAA Journal, Vol. 28, No. 4, pp. 618-624.
• Sandham, N.D. and Reynolds, W.C., 1990, "Three-dimensional simulations of large eddies in the compressible mixing layer", Journal of Fluid Mecanics, Vol. 224, pp. 133-158.
• Souza, L.F., Mendonça, M.T. and Medeiros, M.A.F., 2002, "Assessment of different numerical schemes and grid refinement for hydrodynamic stability simulations", Proceedings of 9^th Brazilian Congress of Thermal Sciences and Engineering, ABCM.
• Souza, L.F., Mendonça, M.T., Medeiros, M.A.F. and Kloker, M., 2002, "Analysis of Tollmien-Schlichting waves propagation on a flat plate with a Navier-Stokes solver", Proceedings of 9^th Brazilian Congress of Thermal Sciences and Engineering, ABCM.
• Souza, L.F., 2003, "Instabilidade Centrífuga e Transição para Turbulência em Escoamentos Laminares sobre Superfície Côncava", PhD. thesis, Instituto Tecnológico de Aeronáutica.
• Squire, H.B., "On the stability of three-dimensional distribution of viscous fluid between parallel walls", Proceedings of Royal Society of London -A, Vol. 142, pp. 621-628.
• Strang, G., 1988, "Linear Algebra and its Applications", Brooks/Cole.
• Stuart, J.T., 1967, "On finite amplitude oscillations in laminar mixing layers", Journal of Fluid Mecanics, Vol. 29, pp. 417-440.
• Thompson, K.W., 1987, "Time Dependent Boundary Conditions for Hyperbolic Systems", J. Comput. Physics, Vol. 68, pp. 1-24.
• Vichnevetsky, R. and Bowles, J.B., 1982, "Fourier Analysis of Numerical Appoximations of Hyperbolic Equations", SIAM, Philadelphia.
• White, F.M., 1974, "Viscous Fluid Flow, (1^st Edition)", McGraw-Hill, New York.

Publication Dates

• Publication in this collection
  26 Aug 2009
• Date of issue
  June 2009

History

• Received
  Mar 2009