**NUMERICAL SIMULATION OF TAYLOR–GREEN 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**

** **

Mathematical model and numerical solution method

Comparison of the results for different Re and the energy spectrum

** **

** **

**Abstract**

We report the results of numerical modeling and visualization for the classical Taylor–Green 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 Taylor–Green 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 Taylor–Green 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, Taylor–Green vortex.*

The evolution of a vortex flow defined at the initial time as the Taylor–Green 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 Taylor–Green vortex has been used as a benchmark solution for testing and validation of numerical algorithms designed for turbulent flow simulation (see, for example, [1–5]). 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 Taylor–Green flow. Methods for constructing the QGD equations and examples of their use can be found, for example, in [6–8] and more recent publications. Some numerical results concerning laminar-turbulent transition simulation based on the QGD equations were presented in [9–13]. 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 Taylor–Green 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].

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 Taylor–Green 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/(kg·K). 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/(m·s) 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 Taylor–Green 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 [11–14]. 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 128^{3} (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 [11–14], 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 Taylor–Green flow problem [1–3]. 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 [11–14]. Note that our code is portable across different platforms
supporting the C language and the MPI standard.

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 128^{3}
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 Taylor–Green 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 Taylor–Green 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 64^{3} points with m and on a grid of 128^{3} 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 Kolmogorov–Obukhov 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 Taylor–Green flow have
been presented in a number of studies, for example, in [1, 3–5]. Specifically,
the computations in [5] were preformed on a grid consisting of 256^{3} 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 128^{3} grid (Fig. 2) agree well with the data of [1,
5] obtained on a 256^{3} 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 64^{3} 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
Taylor–Green 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 64^{3} grid gives
a maximum dissipation rate of ~ 0.01 (Fig. 2b).

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 128^{3} 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 96^{3} degrees of freedom. The number
of computational points in the QGD algorithm is 64^{3} (with m, dash-dotted curve) and 128^{3} (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.

Figure 5 presents
level surfaces of the -component of vorticity computed at Re=100. Here, the grid consists of 128^{3} 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 32^{3}
points (with *h*=10^{-3} m), 64^{3} points (with m), and 128^{3}
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 Taylor–Green flow at Re=100 can be fairly well
simulated even on the very coarse grid with 32^{3}
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 Taylor–Green 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. 1d–1f). This pattern of the Taylor–Green 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 256^{3}
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 64^{3} and 128^{3} 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****.**

The DNS approach was used to simulate Taylor–Green 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 Taylor–Green 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 Taylor–Green 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.

[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. 411–452.

[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. 273–311.

[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 laminar–turbulent transition in the flow over a backward-facing step]. Vestnik Moskovskogo Universiteta, Ser. 3, Fiz. Astron., no 4, 2007, pp. 14–17.

[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, 14–18, September 2009. Edited by K. Hanjalie, Y. Nagano, S. Jakirilic.
Begell House Inc., p. 383–386.

[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)* 4–10 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. 410–418.

[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. 27–48.

[14] T.G. Elizarova, I.A. Shirokov. Testirovanie KGD–algoritma 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. 1–26.
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-breaker’s tank after
impact interaction with ice barrier. *Scientific Visualization*, vol. 5, **no** 4, 2013, pp. 118–135. 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].