## Abstract

The results of three-dimensional numerical simulations of turbulent flows obtained by various authors are reviewed. The paper considers the turbulent mixing (TM) process caused by the development of the main types of instabilities: those due to gravitation (with either a fixed or an alternating-sign acceleration), shift and shock waves. The problem of a buoyant jet is described as an example of the mixed-type problem. Comparison is made with experimental data on the TM zone width, profiles of density, velocity and turbulent energy and degree of homogeneity.

## 1. Introduction

There are several approaches to the simulation of the turbulent flows of a multi-component compressible medium using hydrocodes—programs that approximate the Euler or Navier–Stokes equations.

The first and, maybe, the most reliable way from the viewpoint of accuracy is the direct simulation of such flows using three-dimensional hydrocodes without any turbulence models. Such a simulation method is called direct numerical simulation (DNS). We use the term DNS independently of the equations solved in a code, Euler or Navier–Stokes, though a number of authors use this term only for the Navier–Stokes equations. As a rule, in such calculations the initial density, velocity or interface perturbations are set. Then the perturbations are rather small and non-regular, and usually random-number generators are used for setting such perturbations. Such calculations were initiated not long ago, because till recently there have not been either computers with sufficient performance and memory or appropriate three-dimensional hydrocodes. This approach has significant restrictions related to the Reynolds number, which in this case depends on the number of computational cells. Therefore, it is suitable only for relatively small Reynolds numbers, though, as computer capabilities grow, the corresponding Reynolds numbers are growing with them.

However, the researchers had to hand two-dimensional hydrocodes, and so the calculations were performed in two-dimensional approximations. Although turbulent flow is essentially a three-dimensional phenomenon, this approach gave good results on the description of the turbulent flow development law. It was discovered that large eddies containing the major part of the turbulent energy were described rather adequately by two-dimensional codes. The first two-dimensional calculations of gravitational turbulent mixing (TM) were performed in [1], and two-dimensional calculations of shear TM were conducted in [2].

Nevertheless, we should admit that even nowadays modern computers have insufficient performance for direct three-dimensional simulation of turbulent flows for a rather wide range of problem types. They can be used only as a support tool for a deeper understanding of the processes occurring during TM.

From this standpoint, the second approach, namely the large eddy simulation (LES) method, based on direct simulation of large-scale eddies and the use of subgrid models for the simulation of eddies with scales less than the computational grid size, is more promising. However, for a compressible medium, such an approach has not yet been significantly developed; there are some works also showing that the effects of subgrid viscosity are small in comparison with the scheme viscosity of the difference scheme [3,4].

There are LES calculations performed without any subgrid models—the so-called implicit LES (ILES). Note that some authors often do not distinguish the ILES and DNS methods. It was the ILES method that was used in the 1990s when the appropriate codes were developed. The first calculations of gravitational mixing in a three-dimensional setting were performed in [5].

The low-cost approach to turbulence simulation, in terms of computational resources and availability, is the third one—the development of turbulence models (empirical or semi-empirical) and their implementation into hydrocodes (so-called RANS models—Reynolds averaged Navier–Stokes).

Attempts to describe turbulence using semi-empirical models have a rather long history and are described in detail in different publications. This mainly concerns the shear turbulence, occurring as the result of liquid layers’ relative sliding. Some complex models that include tens of parameters and allow the description of fine experimental results are known. However, because of their complexity, they are difficult to realize. Applications of simple models with a few parameters are strongly limited. In real-life conditions, one has to choose between the general character and complexity, on the one hand, and simplicity and realizability, on the other.

This paper describes the investigations performed by its authors using DNS of TM in the main types of elementary unstable flows, such as gravitational, shock wave and shear ones. Apparently, then, we should bear in mind that neither of the instabilities is met separately in the uncombined state; as a rule, there is a combination of their action, and therefore their separation is rather a conditional measure.

This review does not consider the issues of code development for direct simulation; it is devoted to the physical aspects of turbulent flows. All the investigations discussed below have been performed using the EGAK code for two-dimensional simulation and the TREK code for three-dimensional simulation. All the information on the use of the codes can be found in [6,7].

The review does not claim completeness and all references to the works by other researchers are given to indicate priorities, or for comparison with our own results. There are papers in which some of the problems considered here are discussed in more detail. For example, results of the numerical simulation of Rayleigh–Taylor turbulence are discussed in papers [8–10].

## 2. Gravitational mixing

The most comprehensive study using DNS was performed for gravitational TM, that is, mixing due to the Rayleigh–Taylor instability at constant gas–gas (liquid–liquid) interface acceleration from a heavier gas towards the lighter one.

### (a) Mixing without account for viscosity

This problem has been considered in a number of papers and has abundant information, both experimental and computational. A detailed description of this problem made by the authors of this paper is presented in [11–15].

#### (i) Calculation set-up

At the initial moment, two half-spaces divided by the plane *z*=*z*_{c}=0 are filled with liquids with the densities *ρ*_{1} and *ρ*_{2}=*nρ*_{1} (*n*>1). The initial geometry is shown in figure 1. The gravitational acceleration is directed from the heavy material towards the lighter one. The complete problem set-up can be found in [15]; the computation grid is 200×200×400.

#### (ii) Calculation results

The bitmap images of the flow typical for all the problems considered in this paper are shown later. In the given problem, the issue of transition to the self-similar mode and the constant of the turbulent mixing zone (TMZ) width self-similar growth are of primary interest. The self-similarity of this problem is expressed, in particular, by the change-over to the linear dependence of the TMZ width,
1.1where *τ*=*t*/*t*_{L}, , *L*_{t} is the full TMZ width in the *z* direction after averaging horizontal planes and *A* is the Atwood number. The d*F*/d*τ* inclination angle determines the value of the coefficient *α*_{a}=(d*F*/d*τ*)^{2} in the formula for the TMZ width at the self-similar phase,
1.2Figure 2 shows the dependence *F*(*τ*). For *n*=8.5–29, the values of *α*_{a} obtained in different experiments (see paper [16]) are in the range 0.15–0.2, slightly growing with *n*. The graphs corresponding to these values are shown in figure 2. As we can see in figure 2, there are two sections of the flow, which differ from each other by the inclination angle d*F*/d*τ*. In segment one of the graph this angle is larger than the inclination of the experimental dependence; in the second segment, vice versa, it is smaller. Such results were numerically obtained in [17]. The detailed flow analysis made in [14,15] shows that, in the first segment, the self-similarity of the flow is not present, though the inclination of d*F*/d*τ* is close to a constant. The self-similarity sets on only at the second phase, and therefore the value of the self-similar constant *α* should be measured in this segment; it turns out to be less than the known experimental values. The authors suppose that this is due to the fact that, during the processing of these data, the initial self-similar phase was not excluded. Note that the non-self-similarity of the initial segment can also be connected with the presence of interface local perturbations (LPs) at the initial moment (see §1*b*).

Now, let us illustrate the above using the experimental data of [18]. Figure 3 shows the time dependence of the function of the coordinate *z*_{2} of the heavy liquid penetration into the lighter one. As we can see, for *n*=40 at the final phase when the self-similar mode is achieved, the calculation leads to an inclination angle of d*F*_{2}/d*τ* that is close to the one observed in the experiment—figure 4 shows the minimum and the maximum values for *n*=36.5. Note that, using the dependence *F*_{2}(*t*) obtained over the value *α*_{2}=0.078 (*α*_{2}=(d*F*_{2}/d*τ*)^{2}) accepted in [18], the obtained inclination angle is larger than the observed one. This value has been, evidently, obtained without the exclusion of the initial segment where the self-similarity is not present. Note that there is the same dependence for other density differences.

Now, we are going to analyse the flow from the standpoint of achieving the self-similar mode. Figure 5 shows the reduced TMZ width
1.3following from the performed calculations. Here, the quantity *t*_{0} corresponds to the intersection of the extrapolated self-similar segment of the curve *F*(*t*), evaluated under the data of figure 2, with the abscissa axis.

In the calculations with *n*=10, 20, 40, the quantity *α* only after some time arrives at approximately constant values (which corresponds to the self-similar mode); *α*_{a}≈0.11, 0.15, 0.16, respectively. For this problem, the self-similar mode also expresses itself in reaching the stationary value of *E*(*t*)≡*k*_{m}/*L*_{t}*g*, where is the maximum value, over the TMZ width, of the averaged turbulent energy. The dependence *E*(*t*) is shown in figure 5. We can see that, in the calculations at the initial phase, the values of *E* are rather high (at this phase, the TMZ occupies a small number of the calculation cells, which leads to incorrect values of *E*) and increase along with *n*. At a later stage in the variants with *n*=3, *n*=10 and *n*=20, the quantity *E* reaches the approximately constant value *E*=*E*_{a}.

For a more grounded conclusion whether there is a transition to the self-similar mode, we should also use other turbulent quantities: in paper [15], the dependences of the maximum density and turbulent mass flow pulsations in TMZ, which also reach constant values, though after some time. In addition, the turbulent energy spectral analysis suggests that, in the calculations at the self-similar phase, the Kolmogorov law is well described.

Thus, the analysis suggests that, at the initial stage of the flow, the self-similar mode does not set on, and this segment should be excluded from the consideration of the flow self-similar characteristics. The same should be done with the experiments, where uncontrolled initial perturbations, which can influence the studied flows for a rather long time, are always present.

### (b) Turbulent mixing effect on the evolution of the interface local perturbations

In the problem considered above, the character of the initial perturbations provided one-dimensional (on average) flow evolution along the whole interface. The other extreme case of this interface evolution is realized if the initial perturbation has a local character. The papers [19–21] study the LP evolution on the accelerated interface. In [21], the self-similar solution of the LP evolution at constant acceleration (for the two-dimensional case) is obtained. It is shown that the perturbation configuration remains the same, and the perturbation amplitude grows proportionally to *gt*^{2}. The depth of the light material penetration into the heavy one for *A*∼1 is determined by the law *R*_{2}≅*βgt*^{2}/2, where *β*≅0.27 is a constant, which determines the growth of the ‘two-dimensional channel’ vertex. The growth constant obtained in the numerical calculations differs from its analytical value, and equals *β*≅0.44. The TMZ growth in the self-similar mode goes on under the same law as the LP growth, that is, the depth of the light material penetration in the heavy one is described by the quadratic dependence as well, where, however, the value of the self-similar constant is considerably lower.

Recently, several papers have appeared where the issue of the TMZ effect on the LP development is studied. Experimental study of this problem is described in [22], the numerical one in [23,24]. A numerical study of this problem was made in two-dimensional plane geometry simulating the long channel development located at the contact plane, and in a two-dimensional cylindrical set-up simulating the growth of the semi-spherical LP on the intersection of the axis with the contact plane. Evidently, the two-dimensional set-up of the perturbation development differs from the real three-dimensional case; therefore, the obtained results have qualitative character and require validation by three-dimensional calculations.

Figure 6 shows the system configuration at different moments in the calculations with LP as a long channel for *R*_{10}=5 mm (*R*_{10}/*R*_{turbo}=10) and turbulent zone. Here, *R*_{10} is the radius of the LP; *R*_{turbo} is the amplitude of background random perturbations on the interface. The rate of bubble growth with the initial radius 5 mm on the linear segment of the dependence *R*(*S*) is *χ*=d*R*/d*S*=0.3. The rate of the TMZ width growth in this problem is the same as in the problem without LP. Thus, the TMZ does not significantly affect the LP growth rate.

However, for smaller initial LP amplitudes, the results are different. The rate of the LP amplitude growth decreased as the relation *R*_{10}/*R*_{turbo} decreased.

The above is true for a three-dimensional semi-spherical perturbation. As in the case of a cylindrical perturbation, in this case, the TMZ growth rate tends to a constant. In the calculations with three-dimensional LP (as in the experiments), the LP growth decelerates depending on *R*_{10}/*R*_{turbo}. Figure 7 shows the dependences of the LP growth constant on *R*_{10}/*R*_{turbo}. As we can see for the two-dimensional perturbation, *β*=*f*(*R*_{10}/*R*_{turbo}) grows smoothly within the interval *R*_{10}/*R*_{turbo}≈0 to 7, and then rather sharply increases up to the extreme value *β*_{max}=0.3, where LP develops independently of the TMZ. The same picture is observed for the three-dimensional perturbation; however, the extreme value in this case is *β*_{max}≈0.6. Figure 7 shows the experimental dependence for the three-dimensional perturbation.

Thus, for both types of perturbations, different modes of perturbation growth are realized depending on the relation of the initial LP amplitudes and background perturbations on the interface.

### (c) Account for molecular viscosity

It has been shown above that, in three-dimensional calculations made without molecular viscosity, the initial non-self-similar phase of the TM is present, as well as in the experiments. The width of the corresponding section depends on the input data, on the used difference scheme and on the number of computation cells in the calculation, which indicates the considerable influence of the scheme viscosity on the solution. In view of this, the issues of the flow dependence on the molecular viscosity at the initial phase and of the comparative value of the molecular and the scheme viscosities in the calculations are of interest. In papers [25,26], the results of the three-dimensional calculations of the problem from §1*a* are given, taking account of the molecular viscosity. The experiments with viscous liquids are described in [18,27].

#### (i) Calculation results

Figure 8 shows the dynamics of the TM development in the form of isosurfaces of volume fraction at several moments in the calculations with the molecular viscosity coefficients *ν*=5×10^{−6} and *ν*=5×10^{−3}. Figure 8 demonstrates the effect of the molecular viscosity on the mixing process, that is, as the viscosity grows, the small-scale spectrum in the TMZ is presented less clearly.

Figure 9 shows the dependence *F*(*τ*). As we can see, at rather small values of the viscosity *ν*≤5×10^{−5} (the corresponding Reynolds number determined as *Re*_{ν}=*g*^{2}*t*^{3}/*ν* is *Re*_{ν}∼(0.5–4)×10^{6} for *t*=3–6), the calculation results with *ν*=5×10^{−5} and *ν*=5×10^{−6} are close to each other and are close to the calculation made without viscosity. This means that the implementation of the molecular viscosity *ν*≤5×10^{−5} into three-dimensional calculations provides an effect comparable with the scheme viscosity or less; and only the molecular viscosity *ν*>5×10^{−5} exceeds the scheme viscosity in terms of its effect on *F*(*τ*). For *ν*≤5×10^{−5}, the self-similar phase of the mixing is achieved. Figure 9 also shows the approximation with *α*_{a}=0.055; we can see that, at the self-similar phase, it is very close to the calculations made both with and without accounting for the viscosity. On the other hand, at rather high viscosity *ν*≥5×10^{−4} (*Re*_{ν}∼(0.5–4)×10^{5} for *t*=3–6), the results obviously differ from the preceding variants: the higher the viscosity, the longer is the initial (non-self-similar) phase with a large d*F*/d*τ* inclination angle. In fact, there is no transition to the self-similar mode in the calculations with high viscosity. The point is that rather high Reynolds numbers correspond to this phase, whereas in the calculations with higher viscosities such numbers are just not realized. In paper [26], these issues are considered in detail.

### (d) Alternating-sign acceleration

The TM problem on a plane interface between two incompressible liquids under the effect of alternating-sign acceleration which at the initial phase creates an unstable situation (this phase is considered above) and at the following phase a stable one is considered. This problem has been experimentally studied and described in [28,29]; and its numerical investigation is described in [30–35].

The experiments have shown that, in the case of immiscible liquids, at the stable flow phase the liquids are separated from each other; that is, after the first phase when the liquids are mixed up, they return to their initial state. However, the calculations of [30] performed using DNS in the one-liquid approximation and of [31] performed using the *k*–*ε* turbulence model failed to describe the separation process. Only the two-liquid approach used in [32–35] gave a positive result.

The results of the two-liquid calculations from the mentioned paper are shown in figure 10 as bitmap images in one horizontal cross section of the light material volume fraction for different times; figure 11 shows the dependence *F*(*τ*) obtained from the calculations. Figure 11 suggests that in the first phase the behaviour of *F*(*τ*) in all variants is nearly the same. At the second phase, in the calculations with one material, we can see only an insignificant decrease of *F*(*τ*), and after that it slowly grows. At the same time, in the calculations with two materials, the function *F*(*τ*) decreases linearly, so that its inclinations are nearly the same in both variants at different grids. The inclination angle is *α*^{(−)}≡d*F*/d*τ*≈−0.084—the corresponding linear approximation is shown in figure 11.

### (e) Homogeneous mixing degree and probability density function

The degree of mixing homogeneity in turbulent flows is very important for the determination of reaction rates—both chemical and nuclear. The same is true for the material concentration probability density function (PDF; usually for a heavy material).

In papers [36,37], the results of the above-discussed calculations performed in the two-liquid approximation, as well as the additional calculations in the one-liquid approximation, are considered from the standpoint of obtaining the homogeneity degree and the PDF. Note that such calculations are also presented in other papers [8,38–42]; table 1 shows the known calculated and experimental data on the homogeneity degree, *θ*, of the flow under consideration.

Table 1 suggests that there is, first, a considerable dispersion of the calculated data in terms of the homogeneity degree, and, second, considerable disagreement of the calculated and experimental data. Some authors [36,37] tried to perceive the nature of these differences. For this purpose, the connection of *θ* with the calculation method (with the tracking of the interface and without it) as well as with the method of data processing was established. The same connection was discovered for the PDF. The calculation results are compared with the experimental ones in [38,39].

The time dependence of *θ* in the point where the density quadratic pulsation reaches its maximum is shown in figure 12 for the two-fluid calculations with different molecular viscosities. In all the variants, initially, *θ*=1, because at the beginning the turbulence is not developed. Only as long as the turbulence develops does the value of *θ* decrease down to the approximately constant value *θ*≈*θ*_{c}. As we can see in all the variants, except two variants with the largest molecular viscosity, these values are close to each other and equal *θ*_{c}≈0.25–0.3.

Now, we shall consider the results of the one-liquid calculation for the integral over the TMZ value of the homogeneity degree, *Θ* [37]. At the initial phase, *Θ*=1 (figure 13), and then it decreases; however, as it reaches the minimum value, *Θ*_{min}≈0.2, it starts growing and achieves the value *Θ*_{a}≈0.8, typical for the self-similar phase of the calculations without interface tracing (table 1). In similar calculations [8], the increase of *Θ* from *Θ*_{min}≈0.43 to *Θ*_{a}≈0.8 was obtained.

As we can see in figure 12, in the calculations the molecular viscosity growth leads to the decrease of the homogeneous mixing degree—down to *θ*_{c}≈0.15 in the self-similar phase with *ν*=5×10^{−3}. This can be explained by the fact that a rather high value of the molecular viscosity leads to turbulence suppression at small scales of the order of the cell size. As shown in [37], it is the small-scale turbulence that determines the homogeneity degree, and its suppression leads to the decrease of its value.

The same effect should be provided by the molecular viscosity in the experiments. In fact, it is shown in [39] that an increasing gravity, that is, an increasing Reynolds number at the same path lengths, leads to an increased degree of homogeneous mixing.

Papers [36,37] also consider the behaviour of the PDF. It is shown that the PDF distribution, as well as the homogeneity degree, strongly depends on the way the calculation is set up—the one-liquid or the two-liquid one. The first set-up corresponds to physically mixing gases, and the second one to immiscible liquids. It is also shown that, when the results of the two-liquid calculations are spatially averaged, the mentioned quantities behave as in the one-liquid calculation.

## 3. Shear turbulent mixing

The trivial shear flow—a one-dimensional non-stationary plane mixing layer (with shear velocity *ω*_{0})—is considered in this section. Under certain conditions, many more complicated flows are similar to it: stationary plane mixing layer, the initial section of plane and round jets, initial mixing phase in a cylindrical vortex, etc. Paper [2] discusses the investigation of this problem based on the direct two-dimensional simulation; papers [43–46] do the same based on the three-dimensional simulation. The TMZ width, the velocity profiles, the turbulent energy and the mean quadratic pulsation of the transverse and longitudinal velocity components are compared with the measurements presented in [47–49]. The spectral analysis of the velocity pulsations in the TMZ is made, and the PDF is obtained. Some results of numerical studies are shown below; the calculation set-ups are specified in the mentioned papers.

Figure 14 shows the bitmap images of the velocity longitudinal component at two moments of time for one of the calculations made in the transverse direction. The time dependence of the TMZ, , is shown in figure 15, where , *τ*=*tω*_{0}/*L*_{x}. The corresponding straight line, the linear approximation with (the closest one in this variant to the three-dimensional simulation [45,46]) is shown in figure 15. The value –0.1 is true for several variants. Figure 16 shows the profiles of the scaled velocity and turbulent energy ().

## 4. Shock wave mixing

### (a) Shock wave interaction with the TMZ

Here, the turbulence evolution problem, occurring on the gas–gas interface in a shock tube at the interaction with shock waves, is considered. Experimentally, the problem has been investigated in [50,51]; in [52,53], the DNS results are presented.

The calculated data presented below are taken from [52]. The calculation set-up is presented in the same paper. Figure 17 shows the graph of the TMZ boundary dependence on time. We can see that the calculated TMZ widths agree well with the experimental data. Figure 17 also shows the *R*–*t* shock wave diagram. Note that significant growth of the TMZ width takes place when the first shock wave passes through it (*t*_{r1}≈5.8); a more intense TMZ width growth is observed when the second shock wave passes through it (*t*_{r2}≈6.4).

In the experiments, the pulsations of the velocity longitudinal component were measured by five sensors. Figure 18 shows the time dependence of the square longitudinal velocity component pulsation, , obtained in the calculations for the point with the coordinate 13.9 cm compared with experimental data. The results are averaged over the transverse cross section of the shock tube. More details on this problem can be found in [52].

### (b) Turbulent mixing development behind the front of the shock wave with high Mach number

Here, the problem of TM development behind the shock wave with high shock wave Mach numbers is considered. This problem is interesting because of the fact that the shock wave front can merge with the TMZ. In [54], this problem is considered experimentally, whereas in [55], computationally. The passing of the shock wave with the Mach numbers (*Ma*_{SW}>5) through an air–SF_{6} interface is considered. The complete problem and calculation set-ups are presented in the specified papers. It was obtained in the calculations, as well as in the experiments, that, at the Mach number higher than approximately 7, the shock wave merges with the TMZ boundary. This results in shock wave instability; thus the shape of the interface and the shock wave front correlate with each other. As an example, figure 19 shows the isosurfaces of the SF_{6} volume fraction over the level 0.99 and pressure 20 GPa, corresponding to the shock wave location.

## 5. Buoyant jet

In paper [56], the buoyant jet is studied under the condition of vertical inflow of heated air from a container into the atmosphere. This problem is interesting due to two factors: first, there is a joint action of gravitational and shear mixing; second, there are a lot of experimental data on the turbulence characteristics for this problem. In [57], the numerical simulation of this experiment is performed using the *k*–*ε* model; in [58], the three-dimensional simulation is described. The calculated data agree well with the experiments over all the available turbulence parameters. Some results of the three-dimensional investigations are shown below; the experiment and calculation set-ups can be found in the specified papers.

On the whole, the jet development can be observed on the bitmap images (in the cross section *y*=0) of the velocity components for one of the calculations (figure 20). In figure 20, the transverse segments show the cross sections for which the characteristics of the averaged flow and turbulence were provided in the experiments and in the calculations. Figure 21 shows the radial distribution of the mean buoyancy and density fluctuation in the mentioned three cross sections.

## 6. Some applications

The TM simulation approaches described above were used by various authors to solve some application problems. For example, papers [59,60,61] describe the simulation of the burning process of hydrogen–air mixtures in closed rooms of various configurations. These simulations were performed in the interest of providing safety in nuclear power production facilities.

Another example is the problem of plasma interaction with a magnetic field. One such problem was studied during experiments on the generation of superstrong magnetic fields using magnetic explosion generators [62]. Results of such experiments are often determined by the development of small-scale disturbances on the inner surface of the cylindrical liner used to compress a cavity with a magnetic field. The numerical simulation of the evolution of such disturbances due to the liner deceleration by the field is described in paper [63].

A similar problem (the problem of disturbances on a plasma surface due to its deceleration by a magnetic field) occurs in some space experiments with sources of plasma. Such a problem is described in paper [64].

## 7. Conclusion

The works reviewed above in this paper clearly demonstrate the capabilities of direct three-dimensional simulations of TM for different fluid flows. The problem on the agenda is to use the described approach for the investigation of more complex flows that are possible in designing various engineering structures.

## Footnotes

One contribution of 13 to a Theme Issue ‘Turbulent mixing and beyond: non-equilibrium processes from atomistic to astrophysical scales II’.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.