Chapter 2

The TGV is an exact solution of the incompressible Navier-Stokes equations:

∂t + V ⋅∇V = -∇P + ν∇2V,
∇⋅V = 0,

where V = (u,v,w) is the velocity of the fluid in the x, y and z directions respectively, P is the pressure, ρ is the density and ν is the kinematic viscosity. Note that the absence of a forcing function (body force) will result in viscosity dissipating all the energy in the fluid, which will then eventually come to rest.

In order to simulate the flow, the Velocity Correction Scheme will be used to evaluate the Navier-Stokes equations. This scheme uses a projection in which the velocity matrix and pressure are decoupled. For detailed information on the method, you can consult the User-Guide.

Quasi-3D Approach

Instead of constructing a 3D mesh of the domain and applying a spectral/hp discretisation in all three directions, it is more computationally efficient to use a quasi-3D approach. This approach allows the use of a spectral/hp expansion on a 2D mesh and to extend the computation in the third direction (typically the z-direction) using a simpler purely spectral harmonic discretisation (classically a Fourier expansion), if that direction is geometrically homogeneous. The user needs only specify the number of modes and basis type for the discretisation and the length along the z-direction.

Mesh and Boundary Conditions

The mesh we will use for the TGV flow is periodic in nature. It is therefore defined as -πL ≤ x,y,z ≤ πL, where L is the inverse of the wavenumber of the minimum frequency that can be captured on our mesh, i.e the largest length scale of our flow, known as the integral scale. It represents the size of the large, energy-carrying eddies responsible for most of the momentum transport. Given that the smallest length scale we can capture on our mesh Δx = (2πL ∕N )1D in one direction, it is clear that decreasing the integral scale L (for a periodic mesh, by decreasing its size) or increasing the number of mesh points N, allows us to capture smaller and smaller flow scales.

Since the TGV flow transitions to turbulence at a certain time, we would hope to capture the smallest turbulent length scales with our mesh. The smallest length scale η in one direction, as defined by Kolgomorov is:

-η ≈ Re- 3∕4

Then the number of points required for a homogeneous mesh in 3D to capture the Kolgomorov scales is (with L = O(1)):

       (L )3   (     )3
N3D  =  --   =  Re3 ∕4   = Re9 ∕4

Thus, for flow simulations of the TGV, which are usually run at Re = 1600, the number of uniform mesh points required to capture the smallest turbulent length scales is 16,190,861. The discretisation used in the current tutorial of N3D = 643 = 262,144 points does not have enough resolution to resolve these scales, and neither does the case with N3D = 1283 = 2,097,152 points, left as an exercise. The case with N3D = 2563 = 16,777,216 points would be sufficient but to simulate the flow with such a fine mesh requires the aid of a larger parallel computer, given the much larger memory requirements of the computation. When the flow is underresolved some stabilisation is typically required and this is why we will use a technique known as spectral vanishing viscosity for the lower resolution runs.

As mentioned above, the mesh used in this tutorial is a uniform square mesh of length -π ≤ x,y ≤ π (L = 1) containing 64 square elements, each with a spectral/hp expansion basis of order 8. This 2D mesh is then expanded in the z-direction to compute the quasi-3D flow. The homogeneous z-direction will be discretised using 64 planes. Boundary and initial conditions are also applied. The latter will be described in the following section. As for the former, given the periodic nature of the flow, periodic boundary conditions are applied on all edges of the mesh. The 2D meshes 64 and 128, showing all quadrature points, are shown in Fig. 2.1:


Figure 2.1 642 and 1282 element meshes.

Initial Conditions

We are solving an initial value problem and in order to march the problem forward in time we need to prescribe initial conditions for all variables. For this case the initial conditions correspond to the TGV flow, given by:

u(x,y,z,t = 0) = V 0 sin (x )
 L-cos( y)
  L-cos(z )
v(x,y,z,t = 0) = -V 0 cos( x )
  Lsin (y )
 Lcos( z)
 w(x,y,z,t = 0) = 0,
P(x,y,z,t = 0) = 0.


Figure 2.2 u and v velocity components on the z = 0 plane at t = 0 using the 643 elements mesh.


Figure 2.3 z-vorticity on the surface z = 0 at t = 0 and vorticity iso-surfaces at t = 0 on the 643 volume domain.

The u and v velocity components of the vector field are shown in Fig 2.2. Note how these fields are formed by contiguous patches of velocity of same magnitude and opposed direction. The resulting flow, when both velocity fields are combined, is presented in Fig. 2.3 as vorticity contours. The vortical flow is formed by contiguous pairs of counter-rotating vortices, with strength Γ and -Γ.