2-D 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].

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}\begin{aligned} 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{aligned}\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}\begin{aligned} 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{aligned}\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.