Abstract

Dynamics of Coherent Vortices in Mixing Layers Using Direct Numerical and Large-eddy Simulations

Jorge H. Silvestrini

Departamento de Matemática Pura e Aplicada -Universidade Federal do Rio Grande do Sul, 91501-970 Porto Alegre RS, Brasil,

e-mail: silvestr@mat.ufrgs.br

Introduction

Since coherent vortices play a crucial role in mass, heat and momentum transport in geophysical and industrial turbulent flows, their identification has been one of the main objectives of research in turbulence theory in the last years. To be characterized as coherent three conditions are required (Lesieur, 1997):

These coherent vortices are normally called by the name of the hydrodynamic instability which originated them (Kelvin-Helmholtz vortices, Görtler vortices), or by their orientation (streamwise vortices), or by their form (hairpins, lambda vortices). Here is presented some numerical evidences of streamwise vortices, their origin and evolution, in transitional and turbulent mixing layers.

Two numerical techniques for the simulation of turbulent flows were used: Direct Numerical Simulation (DNS) and Large-Eddy Simulation (LES). In DNS, all turbulent scales should be simulated explicitly, in three dimensions of space, from the integral scale l down to the Kolmogorov scale (n3/e)^1/4 for free shear flows (in wall shear flows the dissipative scale is the viscous thickness n/u_*). This implies high-order schemes, small time steps, very fine 3D grids and, in practice, low Reynolds numbers, since it can be proved that the total number of degrees of freedom to represent a turbulent flow is of the order of R_l³. This is the main restriction of DNS to simulate turbulent flows of practical interest.

Since, in general, one is interested in the large scales of the flow, which contains most of the information about momentum and heat transfer, LES strategy consists in simulating, explicitly and in three dimensions, all motion larger than a certain cut-off scale. The smaller scales are modeled through a sub-grid model.

Full pseudo-spectral methods (for temporal mixing layers) and pseudo-spectral methods combined with high order compact finite differences methods (for spatial developing mixing layers) are used to solve the incompressible Navier-Stokes equations. Within this context, it is firstly presented, the general formalism of LES carried out in spectral space for the case where periodicity may be assumed in the three directions i.e. the temporal case. Extension to flows where only two directions may be assumed as periodic i.e. the spatial mixing layer, is then briefly described. The subgrid-scale modelisation strategy is explained and some subgrid-scale models defined in spectral and physical space are described. In particular, the new spectral-dynamic viscosity model, is presented. Finally some DNS and LES results of temporal and spatial mixing layers, denoting the origin and evolution of streamwise vortices, are discussed.

Large-Eddy Simulation

In this section, the LES formalism for incompressible flows in spectral space is considered. The detailed description of this formalism may be found elsewhere (Lesieur and Métais, 1996; Lesieur, 1997). For compressible flows, the LES formalism may be found in Comte et.al. (1994), Ducros et.al. (1996) and Silvestrini (1996), among others.

Full Periodic Problem

Let (k,t) and (k,t) be the spatial Fourier transform of the velocity u(x,t) and the scalar q (x,t) fields of an incompressible flow. Assuming periodicity in the three spatial directions, and using pseudo-spectral methods (Canuto et.al., 1988), the conservation equations of mass, momentum and scalar read in spectral space as :

where k=(k_x,k_y,k_z) stands for the wave number vector, F denotes the Fourier transform operator, w is the vorticity vector and P=p/r +u^.u/2, is the generalized pressure. The pressure P is eliminated in eq. (2) by projection on the plane orthogonal to k, in order to respect the incompressibility condition.

The LES formalism introduces the filtering operation:

where (k) is the Fourier transform of the filter function G(k). Here is used the cut-off filter in Fourier space, defined as:

where k_c=p/D is the cutoff wave number associated to the grid mesh D . Introducing the operation (5) in Eq. (1), (2) and (3), the LES equations read:

The terms on the r.h.s of Eqs. (6) and (7) denote the supergrid-scale and sub-grid scale transfers due to nonlinear terms involved in Navier-Stokes equations in Fourier space. The supergrid-scale transfers need no modelling since they can be explicitly calculated in the large-eddy simulation as:

where P is the projector on the plane normal to the wave number vector k. The unknown subgrid-scale transfers t _çkç>kc(k,t) and t^q _çkç>kc(k,t), should be modelled. Following Kraichnan ideas (Kraichnan, 1976), it was proposed to model these transfers with the aid of the spectral eddy viscosity and diffusivity (Chollet and Lesieur, 1981), as:

where models to calculate n_t(k, k_c,, t) and k _t(k, k_c,, t) should be introduced.

Partial Periodic Problem

Now we want to take into account the streamwise developing character of the mixing layer and therefore we need to change the temporal problem to a spatial developing problem, where no periodic conditions may be assumed in the streamwise direction.

Let now û(x,k_2D ,t) and (x,k_2D.,t) be the spatial bidimensional Fourier transform of the velocity and passive scalar fields, where we assume that x is the streamwise no periodic direction and k_2D=(k_y , k_z) is the wave vector defined in yz plane. With this decomposition the mass conservation equation read now as:

and the momentum and passive scalar conservation equations as:

where Ñ =(¶ /¶ x, ik_y , ik_z) and Ù denotes the Fourier transform.

The main difference between the system of Eqs. (12,13,14) with reference to the system Eqs. (1,2,3), is that now the pressure can no longer be eliminated by projection and then we should solve the Poisson equation that arrives when we take Ñ of Eq. 13, that is:

Since we have now two directions where the variables are described in spectral space and one in physical space, two filters should be introduced: the cutoff in spectral space (the 2D counterpart of Eq. (5)) and a top-hat filter for the physical direction, which may be defined as:

The filtering operation applied to Eqs. (13,14, 15) give us the LES equations to be used in the spatial developing mixing layer case.

Subgrid-Scale Models

Here will be discussed some subgrid-scale models developed at the Grenoble turbulence school (for a complete description of these models see Lesieur & Métais, 1996). Others models developed elsewhere may be found in Smagorinsky (1963), Germano et al. (1991) and Ghosal et al. (1995). Comparison of some of these subgrid-scale models in academic tests cases were reported in Comte et al. (1994) and Comte et al. (1995).

In Spectral Space

Assuming a k^-5/3 inertial range at wave numbers greater than k_c, it was proposed (Chollet & Lesieur, 1981) to renormalize the eddy viscosity with the aid of [E(k_c.,t)/k_c]^1/2, where E(k,t) is the three-dimensional kinetic-energy spectrum. More precisely, the eddy viscosity in spectral space writes

with

The constant 0.267 was obtained with the aid of the EDQNM (Eddy-Damped Quasi-Normal Markovian) non-local interactions theory (Lesieur, 1997), using leading-order expansions in powers of the small parameter k/k_c, and assuming that E(k) follows a Kolmogorov law extending above the cutoff. In eq. (20), K(k/k_c) displays a strong over- shoot (cusp-behaviour) in the vicinity of k/k_c=1 (Kraichnan, 1976). This is due to local or semi-local interactions in the neighborhood of k_c. If one goes back to physical space, the plateau part of the spectral eddy viscosity corresponds to a classical eddy-viscosity formulation, which assumes, in fact, a separation of scales between supergrid and sub-grid scales. This is of course wrong, and fixes the limits of the eddy-viscosity formulation. Therefore, the cusp part of the spectral eddy viscosity is important since it contains effects beyond the classical eddy-viscosity concept. The eddy-diffusivity was found to have, qualitatively, the same behaviour, with a corresponding turbulent Prandtl number P_r^t=n _t^¥ /k _t^¥ approximately constant and taken equal to 0.6 (Lesieur, 1997).

The major drawback of the eddy viscosity described by Eq. (17) is that it assumes a Kolmogorov spectrum at the cutoff. This condition is obviously not satisfied in transitional regions, or close to a wall, even at high Reynolds numbers. To avoid this problem, eddy coeficients may now be evaluated in a less restrictive context than previously. Assuming that the kinetic energy spectrum follows a power law E(k)µ k ^-m instead of a Kolmogorov law, it is found (Métais & Lesieur, 1992):

where the Kolmogorov constant C_K=1.4 and the associated turbulent Prandtl number

The model defined by equations (17), (21) and (22) was called the spectral-dynamic model (SDM). The model was used by Lamballais (1996) for LES of turbulent channel flow with excellent results. A full presentation of the model may be found in Silvestrini et al., 1998.

Note finally that eq. (21) is valid only for m£3 . For m>3, the choice was to set the eddy-viscosity equal to zero. From a practical viewpoint, this may be justified by considering that if the kinetic energy spectrum is steep enough, there is no energy pile-up at high wave numbers, so that no subgrid-scale modelling is actually necessary.

In Physical Space

To determine eddy viscosities in physical space, the kinetic energy at the smallest resolved scale D=p/k_c should be measured. One of these local spectra is F_2D(x,t), the second-order structure function of the resolved velocity field, defined as:

and related to the three dimensional kinetic energy spectrum in isotropic turbulence through Batchelor's formula (Batchelor, 1953):

In the case of infinite Kolmogorov spectra, energy-conservation arguments yield the structure-function model (SF model) (Métais & Lesieur, 1992), defined by

As it involves velocity increments instead of derivatives, the SF model has the advantage of being defined independently of the numerical scheme used. It is nevertheless not much better for transitional flows than the Smagorinsky model: low wave number velocity fluctuations corresponding to unstable modes yield n_t’s so large that affects the growth rate of weak instabilities, like Tollmien-Schlichting waves (Ducros, 1995).

One way of remedying this is to apply a high-pass filter onto the resolved velocity field before computing its structure function (Ducros et al. 1996). With a triply-iterated second-order finite-difference Laplacian filter denoted ~, one finds » E(k)/E(k)»40³(k/k_c)⁹ for all k, almost independently of the grid mesh and the velocity field. With the same formalism used for the structure-function model, this yields the filtered structure-function model (FSF model), defined by:

With this model it was possible to perform a LES of a spatially-growing boundary layer (at Mach 0.5) between Re_x=3.3 10⁵ and 1.14 10⁶, which widely encompasses the transition region (Ducros et al.1996).

Mixing Layers Simulations

In this section simulations of temporal and spatial growing mixing layers are presented.

The temporal approximation is obtained by taking a reference frame moving with the average velocity (U₁+U₂)/2 (where U₁ and U₂ are the velocity of the two parallel streams).

Temporally Growing Mixing Layer

A DNS and a LES of temporal mixing layers differing in the Reynolds number will be presented at first. In both cases, the temporal mixing layers are initiated by a hyperbolic-tangent velocity profile, U tanh 2y/d_ito which is superposed a small quasi-2D random perturbation. d_i denotes the initial vorticity thickness. Then we will analyze the influence of the initial conditions by showing coherent vortices formed in two LES which differ only in the initial perturbation added to the base profile.

For all the simulations the computational domain is cubic of side L_x=L_y=L_z=4l_a, where l _a=7d _i=2p /k_a is the wavelength of the most amplified streamwise mode predicted by the inviscid linear-stability theory (Michalke, 1964). Such a domain allows two successive pairings of Kelvin-Helmholtz (KH) vortices during a simulation. Periodic boundary conditions are imposed in the streamwise (x) and spanwise (z) directions, while free-slip boundary conditions are employed for y=±L_y /2, by means of pure sine or cosine expansions. The time derivative is approximated by a third-order low-storage Runge-Kutta scheme (Williamson, 1980). Aliasing errors (Canuto et al. 1988) are minimized by taking more collocation points in physical space (120³) than Fourier modes (96³).

The DNS with an initial Reynolds number of Re_di=Ud_i /n =100, is presented first. This simulation is called DNSQ2DT. Fig. 1 denotes the vorticity structures of the mixing layer by visualization of vorticity lines. The threshold value of the vorticity norm is w_i /3 (w_ibeing the initial maximal vorticity modulus equal to 2U/d_i). At t=35d_i /U, the first pairing of Kelvin-Helmholtz (KH) vortices is complete. At this time, two KH vortices can be seen, with stretching of vorticity lines in between. These streamwise vortices are called hairpins and are characterized (as we will see) by pairs of vortices of different signs. At t=70d_i /U, the end of the simulation, only one KH vortex remains. The side views show the stretching of vorticity lines at an angle near of 45° with respect to a horizontal plane. The origin of these streamwise vortices, which has been observed in laboratory experiments for a long time (Konrad 1976, Bernal & Roshko 1986), may be explained by the intense deformation rate imposed by the KH vortices in the stagnation zones. The KH vortices strain the vorticity lines, which are originally oriented in the spanwise direction, and align them in the streamwise direction.

The stretching of vorticity lines may be analyzed considering the vorticity equation for a perfect fluid:

and assuming that the vorticity in the stagnation region is weak (êêw_j ê<<êêS_ijêê). Eq. (27) gives the main direction of straining of streamwise vorticity of 45° with respect to a horizontal plane.

Now, results from a LES using the spectral-dynamic model (DM) with an initial Reynolds number of Re_di=2000, are presented. This simulation is called LESQ2DT. In this LES, the spectral-dynamic model is used in its "standard " version, defined by equations (17), (20) and (21). The spectrum slope m is calculated at each time step (and at each sub-step of the Runge-Kutta method), from the three-dimensional kinetic energy spectrum, using a least-square method applied to wave numbers ranging between k_c/k<k<k_c.

The time evolution of the slope of the three-dimensional kinetic energy spectrum close to the cutoff (m), and of E(k_c) is presented on Fig. 2. Until t=10d_i /U no eddy-viscosity model is applied since m£3. Between t=10d_i /U (time of KH vortices roll-up) and 20d_i /U (beginning of the first pairing), the slope of the 3D kinetic energy spectrum takes a rather constant value of m=25. After that, and until the end of the simulation, m tends asymptotically to 2. At this time, the Reynolds number based on the local vorticity thickness d (defined as 2U/ï á w_zñ (y=o)ï ) is Re_d =24000.

In Fig. 3 can be seen the time evolution of the non-normalized three-dimensional kinetic energy spectrum, for the simulation LESQ2DT. Mixing-layer experiments at this Reynolds number do possess a very good k^-5/3 Kolmogorov law over a quite long range at large wave numbers. The figure shows a Kolmogorov law only over a short range, whereas the slope is steeper near k_c, (close to –2, in agreement with Fig.2). This is the main disagreement between the spectral dynamic model and experiments in mixing layers.

Fig. 4 shows the vorticity-modulus isosurfaces of the quasi 2D mixing layer at t=35 and 75d_i /U with a threshold value of w_i. The figure shows the moment of the two pairing of KH vortices and the intense stretching of streamwise vortices. Note that the isovalue has increased three times with reference to the previous DNS. At the end of the simulation, there is only one KH vortex. Note finally the presence of intense small-scale vortices.

For the two simulations presented, DNSQ2DT and LESQ2DT, we analyze now the evolution of the "diameter" of the streamwise vortices with viscosity. Fig. 5 shows isolines of vorticity modulus in the same transversal plane and at the same time, for the DNS at Re=100 and for LES at Re=2000. In the DNS we can measure d»2d_iwhile in the LES d»d_i. These values may be compared with an experimental one, d»13d_i, obtained at Re»1400 (Huang & Ho, 1990). The agreement is fairly good and the tendency is correct. This shows the strength of LES to reproduce features of turbulent flows.

Now, we present results of a LES (also with the SDM model) where the initial conditions are defined by the same base profile but the perturbation is a three-dimensional isotropic one. The initial Reynolds number is also Re_di=2000. The simulation was stopped t=60d_i /U and is called LES3DT.

Figure 6 shows the temporal evolution of m and E(k_c) for the whole simulation. The spectrum slope decreases initially from the high initial value (m»9). At t=10d_i/U, which correspond to the time of the vortex roll-up, we have m»3. It means that the eddy- viscosity was inactive (see equation 21) up to this instant, and that all the dissipation was due to molecular viscosity. Hence instabilities are allowed to grow without any influence of the eddy viscosity, which is certainly desirable. Between t=10 and 30d_i/U (moment of the first pairing), the slope m decreases from 3 to 2. After that, m remains very close to 2 up to the end of the simulation. The temporal evolution of E(k_c), which reaches its maximum at t=25d_i /U and then decreases slowly, might indicates that a "quasi- equilibrium state" characteristic of the self-similar regime was attained.

Statistics of the recorded velocity profiles were used to determine the temporal evolution of the local vorticity thickness, and compared with experimental data of spatially-growing mixing layers carried out by Bell and Mehta (1990). The l.h.s. of Figure 7 shows d(t). A fairly good linear growth is established very early at a rate of U^-1dd/dt=0.19. During the first pairing (t»30), the spreading slows down, and then it starts rising again at the same linear rate. In spite of the differences in the spatial growth of mixing layers reported in several works (Silvestrini, 1996), and also between the spatial and the temporal problem, the growth rate found here is very close the traditionally accepted mixing layer spatial growth of 0.18 reported in Brown and Roshko's experiments (Brown and Roshko, 1974).

The r.h.s. of Figure 7 and Figure 8 show, respectively, the mean streamwise velocity and velocity components variances at the end of the simulation (t=60d_i /U). The agreement between numerical and experimental data is good and seems to indicate that a self-similar state has been established at the end of the simulation, as far as mean and variances of velocity are concerned. To confirm this point, normalized three-dimensional kinetic energy spectra are presented on Figure 9 (the normalization is made by U and the local vorticity thickness d ). The good collapse of the different spectra for t=50,55 and 60d_i /U is another good indicator that a self-similar regime is attained. Note also that Bell and Metha have considered that a self-similar regime was established at a streamwise distance of about 250d_ifrom the splitter plate, with the velocity ratio l =(U₁ -U₂)/(U₁+U₂)=0.25. If this distance is transformed in a corresponding elapsed time for a temporal mixing layer, using the convection velocity U_c=(U₁+U₂)/2 and writing

the value found t_self=62.5d_i /U is very close to the time considered for present statistics.

Let us look now at the three-dimensional vortical structure. Figure 10 presents a perspective view of vorticity-modulus isosurfaces (threshold w_i), at t=14,26,40 and 60d_i /U. At t=14d_i /U, one can see a dislocated array of four rolling-up Kelvin-Helmholtz vortices, similar to the configuration found in previous DNS of Comte et al. (1992) and laboratory experiments of Chandrsuda et al. (1978), and called "helical pairing". Secondaries streamwise vortices are also stretched by the deformation field induced between the big vortices. At t=26d_i /U large structures pair. The subsequent pairing is more difficult to identify from the vorticity isosurfaces, mainly because of a rapid growth of small-scale structures. At the end of the simulation (t=60d_i /U), the vorticity field displays only the presence of intense small-scale vortices, with no obvious preponderant orientation. By contrast, the low-pressure field (see Figure 11) indicates the presence of one big quasi two-dimensional vortex, stretching thinner longitudinal vortices. Note however that the computational domain is too small at this instant, with regard to the vortex size.

Spatially Developing Mixing Layer

Here, results from a LES using the filtered structure-function (26) of a spatial growing mixing layer are presented. The simulation is called FSFQ2DS. The numerical code used, that solves Eqs. (13,15), combines pseudo-spectral methods in the spanwise and transverse directions with compact finite-difference of sixth order (Lele, 1992) in the streamwise direction. Free-slip conditions are still imposed upon the boundaries. Non-reflective outflow boundary conditions are approximated by a multi-dimensional extension of Orlansky's discretization scheme. The temporal integration is performed by means of a low-storage 3rd order Runge-Kutta scheme, with a fractional step procedure for the pressure-gradient correction.

The profile prescribed at the inlet is:

plus small-amplitude random perturbations. The velocity ratio is chosen as: R=(U₁ –U₂)=0.5. The domain’s dimension are L_x=112d_i , L_y=28d_i and L_z=14d_i, and the grid mesh is cubic with 384x96x48 collocation points. The last record from a previous DNS run a Re=100 was used to initialize the FSFQ2DS run (Silvestrini, 1996).

Fig. 12 shows an isosurface of vorticity modulus at the end of the simulation, with a threshold value of 2/3w_i. In the figure may be identified, from left to right, two KH vortices undergoing a first pairing, after that a cluster of three KH vortices undergoing also a first pairing, and, at the end, a billow made of 4 fundamental KH vortices, whose second pairing is in progress. Between the KH vortices, streamwise vortices are stretched by the same mechanism observed in temporal mixing layers.

Fig. 13 shows an isosurface of streamwise vorticity. The black and grey colors, denoting pairs of vortex of different sign, identify the "legs" of the hairpins. The figure enables us to analyze the evolution of the diameter of the streamwise vortices with the pairing of the KH vortices. After the pairing of the three KH vortices, streamwise vortices are less numerous but bigger than before the pairing. More quantitatively, if two transversal planes are fixed, before (x=58d_i), and after this pairing (x=88d_i), the diameter of the streamwise vortices can be measured. A loose estimation gives d_bef » 0.5d_iand d_aft » d _i, which may be transformed to the local vorticity thickness d_bef » 0.13d and d_aft» 0.14d with the aid of the streamwise evolution of the vorticity thickness (Silvestrini, 1996). These values seem to reinforce the idea that a pairing of KH vortices induces a merging of streamwise vortices, with their diameter being scaled with the local vorticity thickness d .

But this calculation has not reached self-similarity: kinetic-energy spectra in the down-stream region are in k^-5/2 and r.m.s. velocity fluctuations have a departure of about 20% from experiments. Therefore simulations in longer domain are necessary to understand the downstream evolution of coherent vortices to self-similar conditions. Results from a LES using the filtered structure function model but where the upstream conditions are now perturbed by an isotropic noise is full reported in Comte et al. (1998).

Conclusions

The LES formalism in spectral space for incompressible turbulence was reviewed. Some subgrid-scale models defined in spectral or physical space were discussed and their limitations discussed. LES in temporal and spatial mixing layers show the main characteristic of these flows: the formation and pairing of KH vortices with intense stretching of streamwise vortices in between. These streamwise vortices may also merge when the KH vortices pair. The diameter of these vortices seems to be scaled with the local vorticity thickness suggesting that both kind of merging are related. The strength of LES to reproduce the coherence of the vortex organization in a turbulent flows is to be remarked.

Acknowledgments

This work was developed at the MOST/LEGI team in Grenoble, France. The author is grateful to Prof M. Lesieur, E. Lamballais and P. Comte for useful discussions, and to P. Begou for computational assistance. Calculations were carried out at the IDRIS (Institut du Développement et des Ressources en Informatique Scientifique, Paris). The author is supported by a Research Fellowship from FAPERGS/RS.

Article presented at the 1st Brazilian School on Transition and Turbulence, Rio de Janeiro, September 21-25, 1998.

Technical Editor: Atila P. Silva Freire.

Publication Dates