Taylor-Green Vortex Decay

The Taylor-Green vortex is a canonical test for viscous flow solvers. In the incompressible limit, the vortex decays exponentially due to viscosity, with a known exact solution. We solve it using the compressible Navier-Stokes equations at low Mach number.

Problem Setup

The exact solution for the velocity field is:

\[v_x = -U_0\cos(kx)\sin(ky)\,\mathrm{e}^{-2\nu k^2 t}, \qquad v_y = U_0\sin(kx)\cos(ky)\,\mathrm{e}^{-2\nu k^2 t},\]

where $k = 2\pi/L$, $\nu = \mu/\rho_0$ is the kinematic viscosity.

using FiniteVolumeMethod
using StaticArrays

gamma = 1.4
eos = IdealGasEOS(gamma)
mu = 0.01
Pr = 0.72
ns = NavierStokesEquations{2}(eos, mu = mu, Pr = Pr)

NavierStokesEquations{2, IdealGasEOS{Float64}}(EulerEquations{2, IdealGasEOS{Float64}}(IdealGasEOS{Float64}(1.4)), 0.01, 0.72)

Physical parameters (low Mach number: $U_0 \ll c_s$):

rho0 = 1.0
P0 = 100.0   ## high pressure ensures low Mach
U0 = 0.01
L = 1.0
k = 2 * pi / L
nu = mu / rho0

function ic_tgv(x, y)
    vx = -U0 * cos(k * x) * sin(k * y)
    vy = U0 * sin(k * x) * cos(k * y)
    P = P0 - rho0 * U0^2 / 4.0 * (cos(2 * k * x) + cos(2 * k * y))
    return SVector(rho0, vx, vy, P)
end

ic_tgv (generic function with 1 method)

Solving

N = 32
mesh = StructuredMesh2D(0.0, L, 0.0, L, N, N)
t_final = 0.5

prob = HyperbolicProblem2D(
    ns, mesh, HLLCSolver(), CellCenteredMUSCL(MinmodLimiter()),
    PeriodicHyperbolicBC(), PeriodicHyperbolicBC(),
    PeriodicHyperbolicBC(), PeriodicHyperbolicBC(),
    ic_tgv; final_time = t_final, cfl = 0.3
)
coords, U, t_end = solve_hyperbolic(prob)

32×32 Matrix{Tuple{Float64, Float64}}:
 (0.015625, 0.015625)  (0.015625, 0.046875)  …  (0.015625, 0.984375)
 (0.046875, 0.015625)  (0.046875, 0.046875)     (0.046875, 0.984375)
 (0.078125, 0.015625)  (0.078125, 0.046875)     (0.078125, 0.984375)
 (0.109375, 0.015625)  (0.109375, 0.046875)     (0.109375, 0.984375)
 (0.140625, 0.015625)  (0.140625, 0.046875)     (0.140625, 0.984375)
 ⋮                                           ⋱  
 (0.859375, 0.015625)  (0.859375, 0.046875)     (0.859375, 0.984375)
 (0.890625, 0.015625)  (0.890625, 0.046875)     (0.890625, 0.984375)
 (0.921875, 0.015625)  (0.921875, 0.046875)     (0.921875, 0.984375)
 (0.953125, 0.015625)  (0.953125, 0.046875)  …  (0.953125, 0.984375)
 (0.984375, 0.015625)  (0.984375, 0.046875)     (0.984375, 0.984375)

Comparison with Exact Solution

The exact solution decays exponentially:

decay = exp(-2 * nu * k^2 * t_end)
nx, ny = N, N

max_err_vx = 0.0
max_err_vy = 0.0
for i in 1:nx, j in 1:ny
    w = conserved_to_primitive(ns, U[i, j])
    x, y = coords[i, j]
    vx_exact = -U0 * cos(k * x) * sin(k * y) * decay
    vy_exact = U0 * sin(k * x) * cos(k * y) * decay
    global max_err_vx = max(max_err_vx, abs(w[2] - vx_exact))
    global max_err_vy = max(max_err_vy, abs(w[3] - vy_exact))
end

true

Visualisation

using CairoMakie

xc = [coords[i, 1][1] for i in 1:nx]
yc = [coords[1, j][2] for j in 1:ny]
vx_num = [conserved_to_primitive(ns, U[i, j])[2] for i in 1:nx, j in 1:ny]
vy_num = [conserved_to_primitive(ns, U[i, j])[3] for i in 1:nx, j in 1:ny]
vx_ex = [-U0 * cos(k * xc[i]) * sin(k * yc[j]) * decay for i in 1:nx, j in 1:ny]

fig = Figure(fontsize = 24, size = (1200, 500))
ax1 = Axis(
    fig[1, 1], xlabel = "x", ylabel = "y",
    title = L"v_x \text{ (numerical)}", aspect = DataAspect()
)
hm1 = heatmap!(ax1, xc, yc, vx_num, colormap = :RdBu)
Colorbar(fig[1, 2], hm1)

ax2 = Axis(
    fig[1, 3], xlabel = "x", ylabel = "y",
    title = L"v_x \text{ (exact)}", aspect = DataAspect()
)
hm2 = heatmap!(ax2, xc, yc, vx_ex, colormap = :RdBu)
Colorbar(fig[1, 4], hm2)
resize_to_layout!(fig)
fig

The numerical solution matches the exact viscous decay very well. The maximum $v_x$ error is $(round(max_err_vx, sigdigits=3)), which is small compared to the velocity amplitude $U_0 = $(U0)$.

Just the code

An uncommented version of this example is given below. You can view the source code for this file here.

using FiniteVolumeMethod
using StaticArrays

gamma = 1.4
eos = IdealGasEOS(gamma)
mu = 0.01
Pr = 0.72
ns = NavierStokesEquations{2}(eos, mu = mu, Pr = Pr)

rho0 = 1.0
P0 = 100.0   ## high pressure ensures low Mach
U0 = 0.01
L = 1.0
k = 2 * pi / L
nu = mu / rho0

function ic_tgv(x, y)
    vx = -U0 * cos(k * x) * sin(k * y)
    vy = U0 * sin(k * x) * cos(k * y)
    P = P0 - rho0 * U0^2 / 4.0 * (cos(2 * k * x) + cos(2 * k * y))
    return SVector(rho0, vx, vy, P)
end

N = 32
mesh = StructuredMesh2D(0.0, L, 0.0, L, N, N)
t_final = 0.5

prob = HyperbolicProblem2D(
    ns, mesh, HLLCSolver(), CellCenteredMUSCL(MinmodLimiter()),
    PeriodicHyperbolicBC(), PeriodicHyperbolicBC(),
    PeriodicHyperbolicBC(), PeriodicHyperbolicBC(),
    ic_tgv; final_time = t_final, cfl = 0.3
)
coords, U, t_end = solve_hyperbolic(prob)

decay = exp(-2 * nu * k^2 * t_end)
nx, ny = N, N

max_err_vx = 0.0
max_err_vy = 0.0
for i in 1:nx, j in 1:ny
    w = conserved_to_primitive(ns, U[i, j])
    x, y = coords[i, j]
    vx_exact = -U0 * cos(k * x) * sin(k * y) * decay
    vy_exact = U0 * sin(k * x) * cos(k * y) * decay
    global max_err_vx = max(max_err_vx, abs(w[2] - vx_exact))
    global max_err_vy = max(max_err_vy, abs(w[3] - vy_exact))
end


using CairoMakie

xc = [coords[i, 1][1] for i in 1:nx]
yc = [coords[1, j][2] for j in 1:ny]
vx_num = [conserved_to_primitive(ns, U[i, j])[2] for i in 1:nx, j in 1:ny]
vy_num = [conserved_to_primitive(ns, U[i, j])[3] for i in 1:nx, j in 1:ny]
vx_ex = [-U0 * cos(k * xc[i]) * sin(k * yc[j]) * decay for i in 1:nx, j in 1:ny]

fig = Figure(fontsize = 24, size = (1200, 500))
ax1 = Axis(
    fig[1, 1], xlabel = "x", ylabel = "y",
    title = L"v_x \text{ (numerical)}", aspect = DataAspect()
)
hm1 = heatmap!(ax1, xc, yc, vx_num, colormap = :RdBu)
Colorbar(fig[1, 2], hm1)

ax2 = Axis(
    fig[1, 3], xlabel = "x", ylabel = "y",
    title = L"v_x \text{ (exact)}", aspect = DataAspect()
)
hm2 = heatmap!(ax2, xc, yc, vx_ex, colormap = :RdBu)
Colorbar(fig[1, 4], hm2)
resize_to_layout!(fig)
fig

This page was generated using Literate.jl.