Taylor-Green Vortex

The Taylor-Green vortex (TGV) is a canonical benchmark used to assess a solver’s ability to reproduce unsteady, spatially periodic flows with known reference data. The initial conditions take the form of smooth trigonometric functions, making the TGV an ideal test for equation-based initialization: the solver must evaluate \(\sin\) and \(\cos\) expressions at every lattice node at startup and reproduce them to floating-point precision before the time-marching begins.

Two complementary sub-cases are used here. Case 00a is a 2-D TGV, which admits a closed-form analytical solution valid for all time and therefore provides a stringent quantitative accuracy and convergence test. Case 00b is a 3-D TGV, for which no analytical solution exists once the flow has transitioned to turbulence; the dissipation rate is compared against the pseudo-spectral DNS of van Rees et al.[1].

2-D Taylor-Green Vortex (Case 00a)

The 2-D TGV is a doubly periodic flow on the square domain \([0,\, 2\pi] \times [0,\, 2\pi]\) initialised with a single Fourier mode. Because the nonlinear advection terms are exactly zero for this initial condition, the flow decays purely through viscous diffusion and the velocity field remains a rescaled copy of the initial condition at all times. This property yields an exact analytical solution that the LBM result can be compared against directly.

Initial conditions

The velocity and density fields at \(t = 0\) are:

(1)\[\begin{split}u_x(x, y, 0) &= U_0 \cos(k x) \sin(k y) \\ u_y(x, y, 0) &= -U_0 \sin(k x) \cos(k y) \\ \rho(x, y, 0) &= \rho_0\end{split}\]

where \(k = 2\pi / L\) is the wave number, \(L\) is the domain length, and \(U_0\) is the peak velocity amplitude. For a lattice of \(N\) nodes per side with lattice spacing \(\Delta x = 1\), the wave number in lattice units is \(k = 2\pi / N\).

These fields are prescribed via the models.initialization.equations block introduced in PR #443. Setting use_equation_init: true and providing symbolic expressions for \(u_x\), \(u_y\), and \(\rho\) causes the solver to evaluate them on the GPU before the first time step, placing the initial condition to machine precision without any approximation from interpolation or file I/O.

Analytical solution

The exact solution for \(t > 0\) is an exponential decay:

(2)\[\begin{split}u_x(x, y, t) &= U_0 \exp\!\left(-2\nu k^2 t\right) \cos(k x) \sin(k y) \\ u_y(x, y, t) &= -U_0 \exp\!\left(-2\nu k^2 t\right) \sin(k x) \cos(k y)\end{split}\]

The spatially integrated kinetic energy per unit area is:

(3)\[E(t) = \frac{U_0^2}{4} \exp\!\left(-4\nu k^2 t\right) = E_0 \exp\!\left(-4\nu k^2 t\right)\]

The characteristic decay time is \(t_\nu = (4\nu k^2)^{-1}\). For a well-resolved simulation the solver result should track (2) with a velocity error that decreases as \(O(\Delta x^2)\) when the grid is refined, consistent with the second-order accuracy of the LBM streaming step.

Simulation parameters

Four grid resolutions are run to assess convergence. All cases use the D2Q9 velocity set with the RR-BGK collision operator (the production 3rd-order Hermite operator).

Parameter	Value
Velocity set	D2Q9
Collision operator	RRBGK
Relaxation time \(\tau\)	0.501 - 0.508 (varies with grid size)
Peak velocity \(U_0\)	0.01 (lattice units)
Reference density \(\rho_0\)	1.0
Reynolds number	240 (\(Re = U_0 N / \nu\), constant across grids)
Domain	\([0,\, 2\pi] \times [0,\, 2\pi]\), periodic in both directions
Grid sizes \(N\)	8, 16, 32, 64 nodes per side
Boundary conditions	Fully periodic (periodic_dims: [true, true])
Initialization	Equation-based (use_equation_init: true)

The Mach number \(\mathrm{Ma} = U_0 / c_s = U_0 \sqrt{3} \approx 0.017\) is well below the weakly-compressible limit of 0.1, so compressibility errors are negligible.

Note

The convergence study runs four independent simulations (N = 8, 16, 32, 64). The expected slope on a log-log plot of the \(L_2\) velocity error versus \(\Delta x\) is 2, confirming second-order spatial accuracy.

3-D Taylor-Green Vortex (Case 00b)

The 3-D TGV starts from a smooth, low-amplitude initial condition that is linearly unstable. At moderate Reynolds numbers the flow develops a cascade of smaller vortical structures and eventually reaches a regime of decaying turbulence. Because an exact solution does not exist for \(t > 0\), validation is performed by comparing the kinetic energy dissipation rate \(\varepsilon(t) = -\mathrm{d}E/\mathrm{d}t\) against the spectral DNS results of van Rees et al.[1].

Initial conditions

The velocity and pressure fields at \(t = 0\) are the standard TGV initial condition on the periodic cube \([0,\, 2\pi]^3\):

(4)\[\begin{split}u_x(x, y, z, 0) &= V_0 \sin(k x) \cos(k y) \cos(k z) \\ u_y(x, y, z, 0) &= -V_0 \cos(k x) \sin(k y) \cos(k z) \\ u_z(x, y, z, 0) &= 0 \\ \rho(x, y, z, 0) &= \rho_0 + \frac{\rho_0 V_0^2}{16 c_s^2} \left[\cos(2kz) + 2\right] \left[\cos(2kx) + \cos(2ky)\right]\end{split}\]

where \(k = 2\pi / L\) and \(V_0\) is the peak velocity. The density expression in (4) is the weakly-compressible approximation of the incompressible pressure distribution. In the current simulation \(\rho_0 = 1\) is used uniformly (i.e. the density correction is omitted) because at \(\mathrm{Ma} \approx 0.07\) the maximum perturbation is of order \(10^{-3}\) and does not affect the velocity dynamics.

As with Case 00a, these fields are set via models.initialization.equations, evaluating the trigonometric expressions on the GPU at startup.

Validation approach

The volume-averaged kinetic energy is:

(5)\[E(t) = \frac{1}{N^3} \sum_{i,j,k} \frac{1}{2} \rho(i,j,k,t) \left(u_x^2 + u_y^2 + u_z^2\right)\]

The dissipation rate is estimated by finite difference:

(6)\[\varepsilon(t) = -\frac{\mathrm{d}E}{\mathrm{d}t} \approx -\frac{E(t + \Delta T) - E(t - \Delta T)}{2\,\Delta T}\]

where \(\Delta T\) is the interval between saved snapshots.

The DNS reference of van Rees et al.[1] was computed with a pseudospectral code at \(Re_\Gamma = 1/\nu = 1600\) on grids up to \(768^3\). The reference dissipation curve is digitised from Figure 8a and 8b of the paper (768³ PS results, the most converged resolution). The Reynolds number is defined with characteristic length \(L^* = 1/k_0 = 1\) (the inverse wavenumber), so in lattice units \(Re = V_0 N / (2\pi\nu)\). The peak dissipation time and magnitude are well established and serve as the primary quantitative targets.

Simulation parameters

Parameter	Value
Velocity set	D3Q27
Collision operator	RRBGK
Relaxation time \(\tau\)	0.503056
Peak velocity \(V_0\)	0.04 (lattice units)
Reference density \(\rho_0\)	1.0
Domain	\([0,\, 2\pi]^3\), periodic in all three directions
Grid size \(N\)	256 nodes per side
Boundary conditions	Fully periodic (periodic_dims: [true, true, true])
Initialization	Equation-based (use_equation_init: true)
Reference DNS	van Rees et al.[1], \(Re_\Gamma = 1600\)

Note

The Reynolds number follows the van Rees convention \(Re_\Gamma = 1/\nu\) with characteristic length \(L^* = 1\) (inverse wavenumber) and \(V_0 = 1\). In lattice units this becomes \(Re = V_0 N / (2\pi\nu)\), giving \(\nu \approx 0.001019\) and \(\tau = 3\nu + 0.5 \approx 0.503056\).

[1] (1,2,3,4)

Wim M. van Rees, Anthony Leonard, and D. I. Pullin. A comparison of vortex and pseudo-spectral methods for the simulation of periodic vortical flows at high Reynolds numbers. Journal of Computational Physics, 230(8):2794–2805, 2011. doi:10.1016/j.jcp.2010.11.031.