NUMERICAL SIMULATION OF TAYLORGREEN VORTEX

DECAY IN LAMINAR AND TURBULENT REGIMES

 

I. Shirokov, T. Elizarova

Moscow State University, Moscow, Russian Federation

Keldysh Institute of Applied Mathematics, RAS, Moscow, Russian Federation

ivanshirokov@inbox.ru, telizar@mail.ru

 

 

Contents

 

Introduction

Statement of the problem

Mathematical model and numerical solution method

Numerical results for Re=1600

Numerical results for Re=280

Numerical results for Re=100

Comparison of the results for different Re and the energy spectrum

Conclusions

References

 

 

Abstract

 

We report the results of numerical modeling and visualization for the classical TaylorGreen vortex decay problem basing on smoothed gas dynamic equation system namely quasi- gas dynamic (QGD) equations for viscous compressible gas flow. QGD system can be obtained by temporal averaging of Navier-Stokes system and differs form it by strongly non-linear terms with the small coefficient that has the dimension of time. The QGD system determines the temporal evolution of the gas density, velocity, and pressure as functions of Eulerian coordinates and time. The temperature, pressure, and density are related by the ideal gas state equation. The TaylorGreen vortex decay is examined for nitrogen at normal temperature with Mach number equal to 0.1. Variations in the Reynolds number are produced by varying the flow density and pressure.

 

The QGD system is approximated by an explicit time-differencing scheme. All spatial derivatives are approximated by central differences to second-order accuracy, and the time derivatives are approximated with the first-order accuracy. The computation based on the explicit scheme corresponds to the temporal evolution of the gas dynamic flow. The boundary conditions are periodic in three directions, which means that a system of identical TaylorGreen vortices is contained in an unbounded domain. The periodic boundary conditions are implemented by introducing shadow cells adjacent to the physical boundary.

 

The computations were carried out on the multiprocessor computer system K-100 at the Keldysh Institute of Applied Mathematics of the Russian Academy of Science. The numerical algorithm was parallelized as based on a decomposition of the computational domain by planes x=const with the use of the MPI standard. Our code is portable across different platforms supporting the C language and the MPI standard.

 

It is shown that QGD equations provide a uniform adequate numerical simulation of both laminar and turbulent evolution of the vortex in a unified manner without varying the parameters of the algorithm. For high Reynolds numbers (Re=1600) numerical results show the transition to turbulence, with the creation of small scale structures accomplished by special evolution of the kinetic energy decay with Kolmogorov scaling of kinetic energy spectrum. For small Reynolds numbers (Re=280 and 100) vortex decay remains laminar. Comparison with reference data show, that for turbulent flow QGD method requires less grid points to achieve the same accuracy compared with different variants of high-order DNS methods. On identical spatial grids at Re=1600, the QGD algorithm was found to be more accurate than the LES method with the Smagorinsky subgrid-scale eddy viscosity model.

 

The evolution of the vortex flow obtained in the numerical simulation is visualized by depicting the z-component of the vorticity in a 3D domain. Data visualization demonstrates a development of the 3D flow structures and gives opportunity to compare them for laminar and turbulent flow regimes, including a laminar-turbulent transition.

 

Key words: laminar-turbulent transition, quasi-gasdynamic equations, finite-different schemes, TaylorGreen vortex.

 

 

Introduction

 

The evolution of a vortex flow defined at the initial time as the TaylorGreen vortex (1937) is an illustrative example of an initially single steady vortex that decays with the formation of a cascade of progressively smaller vortex structures, which in turn gradually decay due to dissipation in the flow. This problem has been well studied theoretically at small Reynolds numbers, when the flow remains laminar, and at high Reynolds numbers, when the initial laminar flow transforms into a turbulent eddy cascade. For this reason, the evolution of the TaylorGreen vortex has been used as a benchmark solution for testing and validation of numerical algorithms designed for turbulent flow simulation (see, for example, [15]). A detailed visualization of the numerical results provides more insight into such 3D unsteady processes and ensures their adequate treatment.

 

In this study, we examine the capabilities of the quasi-gasdynamic (QGD) equations as applied to laminar-turbulent transition modeling in the TaylorGreen flow. Methods for constructing the QGD equations and examples of their use can be found, for example, in [68] and more recent publications. Some numerical results concerning laminar-turbulent transition simulation based on the QGD equations were presented in [913]. The effect of QGD dissipation on the spectral characteristics of the turbulent flow developing in the decay of a 3D vortex was examined in [14]. Additionally, an algorithm for constructing kinetic energy spectra was described in detail and the accuracy of the numerical simulation of this problem in the case of mesh refinement and a varying relaxation parameter was discussed.

 

In this paper, the QGD equations are used to compute the decay of the TaylorGreen vortex in nitrogen at normal temperature. Variations in the Reynolds number are produced by varying the flow density and pressure. The numerical results are compared with reference data in laminar (Re=100 and Re=280) and turbulent (Re=1600) vortex decay regimes. Inspection of the numerical results shows that the decay of the initial vortex flow is accompanied by the generation of multiscale vortex structures and their nonlinear interaction. Subsequently, these vortices decay. Additionally, the differences between the vortex structures in laminar and turbulent flow regimes and the moment at which a developed turbulent flow is formed at Re=1600 can easily be seen. The numerical results presented below were obtained on the K-100 multiprocessor computer system [15].

 

It is well known that the numerical simulation of 3D unsteady flows on multiprocessor computer systems produces large amounts of data, which can hardly be analyzed without using a suitable visualization procedure [16]. In this study, the main goal of flow visualization is to depict vortex structures in a 3D flow field, to keep track of their temporal evolution and transformations, and to compare vortex structures in various flow regimes. The last goal is achieved with the help of TecPlot visualization software. An analogous numerical approach and similar software tools were used earlier to compute free-surface flow oscillations in a tank subject to nonuniform movements [17].

 

 

Statement of the Problem

 

Following [3], we consider a cubical domain in Cartesian coordinates, where m. The domain is filled with a gas (nitrogen). The state of the gas is described by the following variables: is the density, , , are the components of the macroscopic gas velocity, and is the pressure. The velocity components are collectively denoted by (similar notation is used for other vectors and tensors).

 

As is customary, the initial flow field for the TaylorGreen vortex is given by (see [1])

 

, , ,

. (1)

 

The initial temperature distribution is uniform in space: K. The initial density is determined by the ideal gas equation of state:

 

. (2)

 

The gas constant for nitrogen is R=297 J/(kgK). Note that the initial flow parameters , , and are also related by the equation of state . The initial Mach number is defined as

 

, (3)

 

where the speed of sound in nitrogen under the initial conditions is m/s and the specific heat ratio is [18]. Since the Mach number in the problem is low, the flow can be considered weakly compressible [3].

 

The Reynolds number is defined as

 

, (4)

 

where kg/(ms) is the dynamic viscosity of nitrogen at the temperature K [19]. Below, the flows computed at Re=100, Re=280, and Re=1600 are compared with the reference data from [1], [4], and [3], respectively. In [2] a similar problem of isotropic turbulence decay was solved in the inviscid case by applying five high-order accurate finite-difference algorithms in the LES approximation, namely, Jameson's multistage scheme, Roe-TVD-MUSCL scheme, and various ENO schemes.

 

The initial values are specified in the following order. The initial characteristic velocity is determined by (3); the density , by (4); and the pressure is given by . Then the initial distributions of velocity and pressure are specified using (1) and the initial density distribution is given by (2).

 

The boundary conditions are periodic in three directions, which means that a system of identical TaylorGreen vortices is contained in an unbounded domain.

 

 

Mathematical Model and Numerical Solution method

 

The gas flow is governed by the macroscopic system of QGD equations [8]. This system determines the temporal evolution of the gas density, velocity, and pressure as functions of Eulerian coordinates and time. The temperature, pressure, and density are related by the ideal gas equation of state (2). The total energy per unit volume E and the total specific enthalpy H are calculated using the formulas and , where is the specific heat ratio.

 

In Cartesian coordinates, the QGD system is written as [8, p. 94]

 

, (5)

, (6)

. (7)

 

The mass flux density is determined as follows [8, p. 97]:

 

, . (8)

 

The viscous stress tensor and the heat flux are given by [8, pp. 83, 96]

 

, (9)

, (10)

, . (11)

 

Here, is the Kronecker delta (for and for ) and is internal energy per unit mass. The thermal conductivity is given by the relation [8, p. 71]

 

 

, (12)

 

where Pr is the Prandtl number. In the problem under study, Pr=0.71 [3].

The dynamic viscosity involved in expressions (9)(11) for and is defined as a function of temperature [18]:

 

, (13)

 

where describes the intermolecular interaction in the gas. In this study, we use [16]. The volume viscosity is set to zero according to [3]: .

 

The relaxation parameter involved in (8)(11) is defined as

 

, (14)

 

where is the local speed of sound and h is the spatial mesh size. The terms containing the coefficient represent subgrid dissipation, which smoothes out fluctuations of gasdynamic variables on scales on the order of the mesh size. The coefficient can be regarded as a tuning parameter determining the degree of subgrid dissipation.

 

To solve problem (5)(14) with initial conditions (1) by a finite-difference method, in the domain we introduce a uniform grid in space and time: , , , , . Note that the number of time steps is not fixed prior to the computations. For all gasdynamic quantities depending on coordinates, we introduce grid functions, for example, for density, and likewise for the other variables. The dimensions of the grid functions are the same as those of the corresponding physical quantities.

 

The problem is approximated by an explicit time-differencing scheme. All spatial derivatives are approximated by central differences to second-order accuracy, and the time derivatives are approximated to first-order accuracy. The algorithm for constructing the difference scheme is the same as in [1114]. The computation based on the explicit scheme corresponds to the temporal evolution of the gasdynamic flow. The periodic boundary conditions are implemented by introducing shadow cells adjacent to the physical boundary. Then, for example, a grid with N=130 is denoted as 1283 (i.e., the shadow cells are not taken into account).

 

The initial free path of gas molecules is defined as , where [16]. In the mathematical model used, it is required that . This condition was satisfied in all the computations presented below.

 

The time step is determined by the Courant condition [8, p. 140]: , where is the Courant number. Note that is involved in expression (14) for the relaxation parameter .

 

As in [1114], the computations were performed in dimensional variables. To compare the numerical results with reference data, the former were nondimensionalized using the reference parameters , , and . Accordingly, the nondimensional time is given by , where s, and the nondimensional kinetic energy per unit volume is .

 

The evolution of the vortex flow obtained in the numerical simulation is visualized by depicting the z-component of the vorticity in a 3D domain. Vorticity visualization is conventional in the study of the TaylorGreen flow problem [13]. It provides an illustrative representation of the flow structures under study and makes it possible to compare numerical results with reference data.

The computations were carried out on the multiprocessor computer system K-100 at the Keldysh Institute of Applied Mathematics of the Russian Academy of Science. The numerical algorithm was parallelized as based on a decomposition of the computational domain by planes x=const with the use of the MPI standard. This technique has been successfully applied in [1114]. Note that our code is portable across different platforms supporting the C language and the MPI standard.

 

 

Numerical results for Re=1600

 

Figure 1 shows level surfaces (isosurfaces) of the -component of vorticity, i.e., the z-component of the curl of the velocity field:

 

.

 

Note that was computed in dimensionless form. The level surfaces are presented for a sequence of dimensionless times from t=0 to t=22.5. The red color corresponds to Vz=0.7, and the blue color, to Vz=-0.7. Here, the grid consists of 1283 points, the mesh size is m, and .

 

(a) (b)

 

(c) (d)

 

(e) (f)

 

Fig. 1.

 

It can be seen that the initially regular velocity distribution (1) (Fig. 1a) breaks up into smaller structures. At the times t=5 and 10 (Figs. 1b, 1c) the velocity distribution retains some features of the initial TaylorGreen vortex, but by t=15 (Fig. 1d) the flow becomes chaotic with nearly no relation to the initial vortex structure. This chaotic velocity distribution is typical of turbulent flows. This flow evolution agrees well with the TaylorGreen vortex evolution observed at Re>500 in [1] and at Re=1600 in [3, 5]. Note that the dimensionless time interval corresponds to the characteristic period of rotation of the initial vortex (1).

(a) (b)

 

Fig. 2,

 

Figure 2a shows the temporal evolution of the kinetic energy per unit volume averaged over the computational domain :

 

.

 

The kinetic energy and time are represented in dimensionless units. The solid curves depict our results obtained on a grid consisting of 643 points with m and on a grid of 1283 points with m, . The dashed curve shows the reference data from [3] produced by a dealiased pseudospectral method and various discontinuous Galerkin methods.

 

Figure 2b displays the kinetic energy dissipation rate as a function of time. As in Fig. 2a, the solid curves depict our results, while the dashed curve shows the reference data from [3]. The maximum value of the dissipation rate indicates the laminar-turbulent transition and the formation of the KolmogorovObukhov spectrum (Fig. 6). Figure 2 shows that the results obtained on the finer grid are much more similar to the reference data than those on the coarser grid.

 

According to the numerical results obtained for in [14], an increase in the coefficient from 0.1 to 0.9 leads to smoother plots of and to the maximum of shifted toward earlier times. A similar effect is observed with a decreasing number of spatial grid points.

 

DNS simulations of the TaylorGreen flow have been presented in a number of studies, for example, in [1, 35]. Specifically, the computations in [5] were preformed on a grid consisting of 2563 nodes by applying methods of high-order accuracy in time and space. The number of nodes used in these computations was estimated as the minimum possible one required for the vortex structures corresponding to a given Reynolds number to be well resolved by applying the DNS approach. The results produced by the QGD algorithm on a 1283 grid (Fig. 2) agree well with the data of [1, 5] obtained on a 2563 grid.

 

The LES approach allows one to reduce the number of computational nodes in the numerical algorithm. The flow simulation based on the LES method employing the Smagorinsky filtering model with 643 nodes is presented in [5]. According to these results, the maximum dissipation rate is understated and reaches values of ~ 0.006, while the reference value is ~ 0.012. Thus, the Smagorinsky LES model is too dissipative for the TaylorGreen vortex flow. Note that a dynamic variant of the Smagorinsky model is used in [5], in which case the magnitude of subgrid dissipation is automatically adapted to the resolved scales in the flow. The QGD algorithm on a 643 grid gives a maximum dissipation rate of ~ 0.01 (Fig. 2b).

 

 

Numerical results for Re=280

 

Figure 3 shows level surfaces of the -component of vorticity. The yellow color corresponds to Vz=0.2, and the blue color, to Vz=-0.2. Here, the grid consists of 1283 points, the mesh size is m, and .

(a) (b)

 

(c) (d)

 

(e) (f)

 

Fig. 3.

 

Inspection of Fig. 3 reveals that the velocity field at Re=280 preserves the general character of the initial regular distribution, although relatively small structures are generated. In this case, the flow pattern becomes partially turbulent.

 

(a) (b)

Fig. 4.

 

Figure 4 displays the temporal evolution of (a) the kinetic energy and (b) its dissipation rate at Re=280. The solid and dash-dotted curves depict our results, while the dashed curves show the reference data from [4] produced by the Fergus pseudospectral code and a discontinuous Galerkin method with 963 degrees of freedom. The number of computational points in the QGD algorithm is 643 (with m, dash-dotted curve) and 1283 (with m, solid curve); in both cases, . As in the case of Re=1600, the results on the finer grid are closer to the reference data than those on the coarser grid, although the differences at Re=280 are rather small.

 

 

Numerical results for Re=100

 

Figure 5 presents level surfaces of the -component of vorticity computed at Re=100. Here, the grid consists of 1283 points, the mesh size is m, and .

 

(a) (b)

(c) (d)

(e) (f)

Fig. 5.

 

(a) (b)

Fig. 6.

 

Figure 6b shows the temporal evolution of the kinetic energy and its dissipation rate computed on three grids consisting of 323 points (with h=10-3 m), 643 points (with m), and 1283 points (with m), . As in the case of Re=280, the numerical results obtained on the finest grid are closer to the reference data than those produced on the coarser grids. The reference profile (dashed curve) was given in [1]. Note that the laminar TaylorGreen flow at Re=100 can be fairly well simulated even on the very coarse grid with 323 nodes.

 

The maximum value of for Re=100 occurs at t=4.5. At this time, the flow still retains some of the initial anisotropy of the TaylorGreen vortex (Fig. 5b). Later, the flow pattern also remains anisotropic. Thus, the whole decay of the single vortex at Re=100 remains laminar, so that no chaotic field of small vortices typical of Re=1600 is formed (Figs. 1d1f). This pattern of the TaylorGreen vortex decay agrees with the data presented in [1] for laminar flows at Re<500.

 

 

Comparison of the Results for Different Re and the Energy Spectrum

 

Figure 7 shows the profiles of the dissipation rate as a function of time for three Reynolds numbers: 100, 280, and 1600. Additionally, Fig. 7 presents the reference data from [1, 4, 3]. As the Reynolds number increases, the maximum of shifts to the right. The results obtained at Re=100 and 1600 are in quantitative agreement with the corresponding curves plotted in Fig. 7 in [1], which were computed on a grid with 2563 points by applying the DNS approach.

 

Note that the QGD model used does not require any modification or parameter tuning in the transition from laminar to turbulent flows. All variants are computed in a unified manner.

 

Figure 8 shows the spectral density of kinetic energy at Re=1600 calculated using the same technique as in [14]. Additionally, Fig. 8 presents a line with a Kolmogorov slope of -5/3. The spectrum was calculated at the time t=8.5, when reaches its maximum value (Fig. 2b). The energy spectrum profile is similar to that for Re=1600 plotted in Fig. 4 in [4] and to the profile for Re=1500 presented in Fig. 8.3 in [5]. The spectra in Fig. 8 were computed on grids with 643 and 1283 points. It can be seen that the energy profiles obtained on both grids are fairly similar, which suggests that the turbulent smoothing is physically adequate. As the number of grid points increases, the spectral density curve becomes smoother and approaches the -5/3 line.

Fig. 7. Fig. 8.

 

Figure 9 shows streamlines and contours of the vorticity Vz in the plane z=0.008 m for (a) Re=1600 and (b) Re=100 at t=20. A developed vortex flow pattern is observed at Re1600, while at Re=100 no small vortex structures are seen. Moreover, the value of Vz at Re=1600 is much higher than that at Re=100. Thus, the difference between the turbulent and laminar cases is easily seen. The symmetry of the computed solution in Fig. 9 is the same as in initial conditions (1), which suggests that the numerical algorithm is correct and stable.

(a) (b)

Fig. 9.

 

 

Conclusions

 

The DNS approach was used to simulate TaylorGreen vortex decay by applying the QGD algorithm. For all considered Reynolds numbers, the numerical results obtained in the laminar and turbulent cases were found to be in good qualitative and quantitative agreement with reference data from several studies, where the TaylorGreen vortex was simulated using various high-order accurate numerical methods (finite-difference, spectral, Galerkin) on finer grids in the framework of the DNS approach. On identical spatial grids at Re=1600, the QGD algorithm was found to be more accurate than the LES method with the Smagorinsky subgrid-scale eddy viscosity model.

 

Based on the QGD algorithm, the DNS simulation of TaylorGreen vortex decay in both laminar and turbulent cases can be implemented in a unified manner without varying the parameters of the algorithm. The results presented suggest that the QGD algorithm is promising as applied to the numerical simulation of a laminar-turbulent transition and turbulent gasdynamic flows.

 

The detailed flow visualization presented in this paper makes it possible to keep track of the transformations of the flow regimes under study, including a laminar-turbulent transition.

 

 

Acknowledgments

 

This work was supported by the Russian Foundation for Basic Research, project no. 13-01-00703-a.

 

 

References

 

[1] M. Brachet, D. Meiron, S. Orszag, B. Nickel, R. Morf, U. Frisch. Small-scale structure of the Taylor-Green vortex. J. Fluid Mech., vol. 130, 1983, pp. 411452.

 

[2] E. Garnier, M. Mossi, P. Sagaut, P. Comte, M. Deville. On the use of shock-capturing schemes for large-eddy simulation. J. of Computational Physics, vol. 153, 1999, pp. 273311.

 

[3] Institute of Aerodynamics and Flow Technology. Available at: http://www.as.dlr.de (German Aerospace Center)

 

[4] JB Chapelier, M. De La Llave Plata, F. Renac, E. Martin. Final abstract for ONERA Taylor-Green DG participation. 1st International Workshop On High-Order CFD Methods. January 7-8, 2012 at the 50th AIAA Aerospace Sciences Meeting, Nashville, Tennessee.

 

[5] D. Fauconnier. Development of a Dynamic Finite Difference Method for Large-Eddy Simulation. PhD thesis. Ghent University. ISBN 978-90-8578-235-3. NUR 978, 928. Wettelijk depot: D/2008/10.500/55.

 

[6] B.N. Chetverushkin. Kinetic Schemes and Quasi-Gas Dynamic System of Equations. CIMNE, Barselona, 2008.

 

[7] Yu.V. Sheretov. Dinamika sploshnykh sred pri prostranstvenno-vremennom osrednenii [Continuum Dynamics under Spatiotemporal Averaging]. SPC Regular and 560 Chaotic Dynamics, Moscow-Izhevsk, 2009, in Russian.

 

[8] T.G. Elizarova, Quasi-Gas Dynamic Equations, Springer, Dordrecht, 2009. IBSN 978-3-642-0029-5.

 

[9] T. G. Elizarova and P. N. Nikol'skii. Chislennoe modelirovanie laminarno-turbulentnogo perekhoda v techenii za obratnym ustupom [Numerical simulation of the laminarturbulent transition in the flow over a backward-facing step]. Vestnik Moskovskogo Universiteta, Ser. 3, Fiz. Astron., no 4, 2007, pp. 1417.

 

[10] T.G. Elizarova, V.V. Seregin. Filtered simulation method for turbulent heat and mass transfer in gas dynamic flows. Proceedings of the 6th International Symposium on turbulence, heat and mass transfer, Rome, Italy, 1418, September 2009. Edited by K. Hanjalie, Y. Nagano, S. Jakirilic. Begell House Inc., p. 383386.

 

[11] I. A. Shirokov, T.G. Elizarova, I. Fedioun, J.-C. Lengrand. Numerical simulation of the laminar-turbulent boundary-layer transition on a hypersonic forebody using the quasi-gas dynamic equations. 4th European Conference for Aerospace Sciences (EUCASS) 410 July 2011, St-Petersburg.

 

[12] T. G. Elizarova, I. A. Shirokov. Chislennoe modelirovanie kolebatelnogo techeniya v okrestnosti giperzvukovogo letatelnogo apparata [Numerical simulation of a non-stationary flow in the vicinity 590 of a hypersonic vehicle]. Matematicheskoe modelirovanie [Math. Model Comput. Simul] vol. 4, no 4, 2012, pp. 410418.

 

[13] I. A. Shirokov, T. G. Elizarova. Pryamoe chislennoe modelirovanie laminarno-turbulentnogo perekhoda v sloe vyazkogo szhimaemogo gaza [Direct simulation of laminar-turbulent transition in a viscous compressible gas layer]. Prikladnaya matematika i informatika [Comput. Math. Model], vol. 25, no 1, 2014, pp. 2748.

 

[14] T.G. Elizarova, I.A. Shirokov. Testirovanie KGDalgoritma na primere zadachi o raspade odnorodnoy izotropnoy turbulentnosti [Testing of the QGD algorithm by a homogeneous isotropic gas-dynamic turbulence decay]. Preprint, Keldysh Institute of 615 Applied Mathematics of Russian Academy of Science, no 35, 2013, p. 19. Available at http://library.keldysh.ru/preprint.asp?id=2013-35.

 

[15] -100. Available at: http://www.kiam.ru

 

[16] A.E. Bondarev, V.A. Galaktionov, V.M. Chechetkin. Scientific visualization in computational fluid dynamics. Scientific Visualization, vol. 2, no 4, 2010, pp. 126. Available at: http://sv-journal.com/2010-4/01.php

 

[17] T. Elizarova, D. Saburin. Mathematical modeling and visualization of the sloshing in the ice-breakers tank after impact interaction with ice barrier. Scientific Visualization, vol. 5, no 4, 2013, pp. 118135. Available at: http://sv-journal.com/2013-4/07.php

 

[18] Bird G.A. Molecular Gas Dynamics and the Direct Simulation of Gas Flows. Clarendon press, Oxford, 1998.

 

[19] Table of Physical Constants: Handbook, Ed. by I. K. Kikoin (Atomizdat, Moscow, 1976) [in Russian].