Problem C3.
5 Direct Numerical Simulation of the Taylor-Green Vortex at Re = 1600
Overview
This problem is aimed at testing the accuracy and the performance of high-order methods on the direct numerical simulation of a three-dimensional periodic and transitional ow dened by a simple initial condition: the Taylor-Green vortex. The initial ow eld is given by u = v w p = = = x y z cos cos , L L L x y z V0 cos sin cos , L L L 0, 2x 2y 0 V02 cos + cos p0 + 16 L L V0 sin
cos
2z L
+2
This ow transitions to turbulence, with the creation of small scales, followed by a decay phase similar to decaying homogeneous turbulence (yet here non isotropic), see gure 1.
Figure 1: Illustration of Taylor-Green vortex at t = 0 (left) and at tnal = 20 tc (right): iso-surfaces of the z -component of the dimensionless vorticity.
Governing Equations
The ow is governed by the 3D incompressible (i.e., = 0 ) Navier-Stokes equations with constant physical properties. Then, one also does not need to compute the temperature eld as the temperature eld plays no role in the uid dynamics. Alternatively, the ow is governed by the 3D compressible Navier-Stokes equations with constant physical properties and at low Mach number.
Flow Conditions
0 V0 L
The Reynolds number of the ow is here dened as Re = 1600.
and is equal to
In case one assumes a compressible ow: the uid is then a perfect gas with c = cp /cv = 1.4 and the Prandtl number is P r = p = 0.71, where cp and cv are the heat capacities at constant pressure and volume respectively, is the dynamic shear viscosity and is the heat conductivity. It is also assumed that the gas has zero bulk viscosity: v = 0. The Mach number used is small enough that the solutions obtained for the velocity and pressure elds are indeed very 0 close to those obtained assuming an incompressible ow: M0 = V c0 = 0.10, p0 where c0 is the speed of sound corresponding to the temperature T0 = R 0 . The initial temperature eld is taken uniform: T = T0 ; thus, the initial density eld is taken as = Rp T0 . The physical duration of the computation is based on the characteristic L convective time tc = V and is set to tnal = 20 tc . As the maximum of the 0 dissipation (and thus the smallest turbulent structures) occurs at t 8 tc , participants can also decide to only compute the ow up to t = 10 tc and report solely on those results.
Geometry
The ow is computed within a periodic square box dened as L x, y, z L.
Boundary Conditions
No boundary conditions required as the domain is periodic in the three directions.
Grids
The baseline grid shall contain enough (hexahedral) elements such that approximately 2563 DOFs (degrees of freedom) are obtained: e.g., 643 elements when using p = 4 order interpolants for the continuous Galerkin (CG) and/or discontinuous Galerkin (DG) methods. Participants are encouraged, as far as can be aorded, to perform a grid or order convergence study.
Mandatory results
The temporal evolution of the kinetic energy integrated on the domain : 1 0 vv d . 2
k = dE dt .
Each participant should provide the following outputs:
Ek =
The temporal evolution of the kinetic energy dissipation rate:
The typical evolution of the dissipation rate is illustrated in gure 2.
Figure 2: Evolution of the dimensionless energy dissipation rate as a function of the dimensionless time: results of pseudo-spectral code and of variants of a DG code. The temporal evolution of the enstrophy integrated on the domain : 1 0 d . 2
E=
This is indeed an important diagnostic as is also exactly equal to 2 E 0 for incompressible ow and approximately for compressible ow at low Mach number (cfr. section 8).
L | | = Figure 3: Iso-contours of the dimensionless vorticity norm, V 0 x t 1, 5, 10, 20, 30, on a subset of the periodic face L = at time tc = 8. Comparison between the results obtained using the pseudo-spectral code (black) and those obtained using a DG code with p = 3 and on a 963 mesh (red).
The vorticity norm on the periodic face illustration is given in gure 3.
x L
= at time
t tc
= 8. An
Important: All provided values should be properly non-dimensionalised: e.g., 3 V2 V2 V0 L 1 divide t by V = tc , Ek by V02 , by L = t0 , E by L0 2 = t2 , etc. 0 c
c
Suggested additional results
In incompressible ow ( = 0 ), the kinetic energy dissipation rate , as obtained from the Navier-Stokes equations, is: = dEk 1 =2 dt 0 S : S d
Furthermore participants are encouraged to provide additional information.
where S is the strain rate tensor. It is also easily veried that this is equal to =2 E. 0
1
If possible, the temporal evolution of the integral k reported in addition to dE dt and the enstrophy.
S : S d shall be
In compressible ow, the kinetic energy dissipation rate obtained from the Navier-Stokes equations is the sum of three contributions:
1
=2
1 0
S d : S d d
where S d is the deviatoric part of the strain rate tensor, v 1 2 ( v ) d 2 = 0 where v is the bulk viscosity and 1 3 = 0
p v d .
The second contribution is zero as the uid is taken with v = 0. The third contribution is small as compressibility eects are small due to the small Mach number. The main contribution is thus the rst one; given that the compressibility eects are small, this contribution can also be approximated using the enstrophy integral. If possible, the temporal evolution of k addition to dE dt and the enstrophy.
1
and
shall also be reported, in
participants are furthermore encouraged to provide numerical variants of kinetic energy dissipation rate (eg. including jump terms for DG methods) and compare those to the consistent values.
Reference data
The results will be compared to a reference incompressible ow solution. This solution has been obtained using a dealiased pseudo-spectral code (developed at Universit e catholique de Louvain, UCL) for which, spatially, neither numerical dissipation nor numerical dispersion errors occur; the time-integration is performed using a low-storage 3-steps Runge-Kutta scheme [2], with a dimensionless timestep of 1.0 103 . These results have been grid-converged on a 5123 grid (a grid convergence study for a spectral discretization has also been done by van Rees et al. in[1]); this means that all Fourier modes up to the 256th harmonic with respect to the domain length have been captured exactly (apart from the time integration error of the Runge-Kutta scheme). The reference solutions are to be found in the following les: spectral Re1600 512.gdiag provides the evolution of dimensionless valk ues of Ek , = dE dt and E . wn slice x0 08000.out provides the dimensionless vorticity norm on the x plane L = . One can use the python script read and plot w.py to extract the data and visualize the vorticity eld or use another language following the format described in the script.
References
[1] W. M. van Rees W. M., A. Leonard, D. I. Pullin and P. Koumoutsakos, A comparison of vortex and pseudo-spectral methods for the simulation of periodic vortical ows at high Reynolds numbers, J. Comput. Phys., 230(2011), 2794-2805. [2] J. H. Williamson, Low-storage Runge-Kutta schemes, J. Comput. Phys. 35(1980)
