Documentation

The 2-D isotropic homogenous turbulence decay problem is a classical and benchmark Direct Numerical Simulation (DNS)problem. The geometry is a simple square domain with a dimension of [2π × 2π] and all four walls have periodic boundary conditions. The energy spectrum E(κ) for this problem follows the enstrophy cascade theory (Kraichnan and Montgomery, 1980) or the local theory of 2-D turbulence (Saffman, 1971) where the scaling power ranges between [κ⁻³,κ^−4.2] in the inertial range (Kida, 1985), κ is the wave number.

This article presents the Direct Numerical Simulation (DNS) implementation of the 2-D incompressible isotropic homogenous turbulence decay problem using OpenFoam (Jasak et al., 2007). The article includes some OpenFoam code snippets for demonstration. In addition, all OpenFoam files required to reproduce the results are attached. OpenFoam-10 is used for this simulation and can be downloaded from https://OpenFoam.org/release/10/.

pimpleFoam (Holzmann, 2016), an OpenFoam incompressible flow solver, is used to numerically solve the continuity Eq. (1) and the momentum Eq. (2).

(1)

(2)

A uniform square grid is generated for the [2π × 2π] square geometry with [1024 × 1024] cells using the blockMesh utility, where all boundaries have cyclic boundary conditions. The blockMesh file is located in the system folder.

The numerical schemes are assigned in the fvSchemes file in the system folder. Backward Euler is used to discretize the temporal derivative term, 3^rd order Gauss scheme is used for the gradient, divergence, and Laplacian operators. The faces’ boundary center variables are computed using 3^rd order interpolation.

The linear algebraic solvers for the pressure/velocity fields and the PIMPLE algorithm controls are defined in the system/fvSolution file. The first pressure field predictor is computed using the Generalized Geometric-Algebraic MultiGrid solver (GAMG) solver with GaussSeidel as a smoother. The final corrected pressure field is computed using the Diagonal Incomplete-Cholesky GaussSeidel (DICGaussSeidel) solver.

The PIMPLE solver’s outer and inner corrector number of iterations is 2. The nonorthogonality corrector iteration (nNonOrthogonalCorrectors) is set to 0 because the mesh is uniform, square, and orthogonal.

The code snippet below is from the parameters files in the system folder, this file is not an OpenFoam standard format but we created it to define all variables in a single place to reduce the file preparation time and errors. In this file, the simulation variables are defined such as the box geometry, mesh size, kinematic viscosityν, total kinetic energy of the initial field Ea, and the peak energy wave number of the initial field κ0. This file is called by other OpenFoam standard files to pull the variables from.

The boundary conditions for the velocity and pressure fields are shown in listings 2 and 3. The initial conditions for both velocity and pressure are set to zero for the entire domain (internalField), these initial conditions will be overwritten later using the boxTurb utility. All sides of the square domain have cyclic (periodic) boundary conditions, since it is a 2-D problem, the frontAndBack faces boundary condition are set to empty. The velocity field file is named U.orig, orig means original, after using the boxTurb utility to generate the initial velocity field, a new velocity file (U) will be created.

In the 2-D turbulence decay problem, the initial random velocity field is generated such that it satisfies the divergence-free constraint, . In OpenFoam, the boxTurb generates this random divergence-free velocity field. The generated velocity field follows the spectrum shown in Eq. (3). In equation Eq. (3), κ is the wave number, κ_p is the peak energy spectrum wave number, and Ea is a constant to scale the total energy contained in the domain.

Ea and κ_p values are assigned in the turbBoxDic file located in the constant directory as shown in listing 4. In the BoxTurbDict file, the system/parameters file is included to pull the Ea and the k0 values.

The turbBox utility originally creates 3-D velocity fields. For 2-D problems, it creates 3-D velocity fields, but after the first time step, the 3^rd velocity component vanishes.

(3)

Other simulation controls are defined in the system/controlDict file, such as the simulation time = 50, initial time step = 0.0001, maximum Courant number =0.7, and the solution writing intervals = 0.2.

The Allrun bash script shown in listing 5 is used to run all the OpenFoam commands automatically. First, it calls the application name from the controlDict file using the $(getApplication) command. then, the grid is generated using blockMesh, then the initial velocity field is generated using boxTurb command. The problem is solved in parallel, hence the decomposePar command is used to distribute the mesh cells to the different processor threads. The number of processors can be controlled from the system/decomposeParDict file. Finally, the runParallel $application command is used to start the pimpleFoam solver in parallel. The last two commands are for postprocessing, such as combining the solutions from the processor threads and computing the enstrophy field, respectively.

Fig. 1 shows the evolution of the vorticity contours for the decay turbulence. At the first time step, the vorticity contours are for the initial random divergence-free field with wave number . With time evolution, coherent large-scale vortices and small-scale vortices form. The size of the large-scale vortices increases with time, and the small-scale vortices dissipate due to the viscous dissipation. This can also be observed from the energy spectrum shown in Fig. 2, at later time steps the inertial range of the spectrum shift to the left (lower wave number or larger vortex size), and the energy spectrum magnitude decrease in the viscous dissipation range at high wave number values. The energy spectrum slope follows aligning with the local theory of 2-D turbulence (Saffman, 1971), (Kida, 1985).