IncompressibleNavierStokes.jl

Appearance

Note: Output is not generated for this example (to save resources on GitHub).

Plane jets - 2D ​

Plane jets example, as presented in [4]. Note that the original formulation is in 3D.

We start by loading packages. A Makie plotting backend is needed for plotting. GLMakie creates an interactive window (useful for real-time plotting), but does not work when building this example on GitHub. CairoMakie makes high-quality static vector-graphics plots.

julia
using FFTW
using CairoMakie
using IncompressibleNavierStokes
using LaTeXStrings

Output directory

julia
outdir = joinpath(@__DIR__, "output", "PlaneJets2D")

Floating point type

julia
T = Float64

Backend

julia
backend = CPU()
# using CUDA; backend = CUDABackend()

Reynolds number

julia
Re = T(6_000)

Test cases (A, B, C, D; in order)

julia
# V() = sqrt(T(467.4))
V() = T(21.619435700313733)

U_A(y) = V() / 2 * (tanh((y + T(0.5)) / T(0.1)) - tanh((y - T(0.5)) / T(0.1)))

U_B(y) =
    V() / 2 * (tanh((y + 1 + T(0.5)) / T(0.1)) - tanh((y + 1 - T(0.5)) / T(0.1))) +
    V() / 2 * (tanh((y - 1 + T(0.5)) / T(0.1)) - tanh((y - 1 - T(0.5)) / T(0.1)))

U_C(y) =
    V() / 2 * (
        tanh(((y + T(1.0)) / 1 + T(0.5)) / T(0.1)) -
        tanh(((y + T(1.0)) / 1 - T(0.5)) / T(0.1))
    ) +
    V() / 4 * (
        tanh(((y - T(1.5)) / 2 + T(0.5)) / T(0.2)) -
        tanh(((y - T(1.5)) / 2 - T(0.5)) / T(0.2))
    )

U_D(y) =
    V() / 2 * (
        tanh(((y + T(1.0)) / 1 + T(0.5)) / T(0.1)) -
        tanh(((y + T(1.0)) / 1 - T(0.5)) / T(0.1))
    ) -
    V() / 4 * (
        tanh(((y - T(1.5)) / 2 + T(0.5)) / T(0.2)) -
        tanh(((y - T(1.5)) / 2 - T(0.5)) / T(0.2))
    )

# U(y) = U_A(y)
# U(y) = U_B(y)
U(y) = U_C(y)
# U(y) = U_D(y)

Random noise to stimulate turbulence

julia
U(x, y) = (1 + T(0.1) * (rand(T) - T(0.5))) * U(y)

# boundary_conditions = (
#     (PeriodicBC(), PeriodicBC()),
#     (PressureBC(), PressureBC())
# )

A 2D grid is a Cartesian product of two vectors

julia
n = 64
# n = 128
# n = 256
x = LinRange(T(0), T(16), 4n + 1), LinRange(-T(10), T(10), 5n + 1)
plotgrid(x...)

Build setup and assemble operators

julia
setup = Setup(x, Re, backend);
# setup = Setup(; x, Re, boundary_conditions, backend);

Initial conditions

julia
ustart = velocityfield(setup, (dim, x, y) -> dim == 1 ? U(x, y) : zero(x));
nothing #hide

Real time plot: Streamwise average and spectrum

julia
function meanplot(state; setup)
    (; xp, Iu, Ip, Nu, N) = setup.grid

    umean = lift(state) do (; u, p, t)
        reshape(sum(u[1][Iu[1]]; dims = 1), :) ./ Nu[1][1] ./ V()
    end

    K = Nu[1][2] ÷ 2
    k = 1:(K-1)

    # Find energy spectrum where y = 0
    n₀ = findmin(abs, xp[2])[2]
    E₀ = lift(state) do (; u, p, t)
        u_y = u[1][:, n₀]
        abs.(fft(u_y .^ 2))[k.+1]
    end
    y₀ = xp[2][n₀]

    # Find energy spectrum where y = 1
    n₁ = findmin(y -> abs(y - 1), xp[2])[2]
    E₁ = lift(state) do (; u, p, t)
        u_y = u[1][:, n₁]
        abs.(fft(u_y .^ 2))[k.+1]
    end
    y₁ = xp[2][n₁]

    fig = Figure()
    ax = Axis(
        fig[1, 1];
        title = "Mean streamwise flow",
        xlabel = "y",
        ylabel = L"\langle u \rangle / U_0",
    )
    lines!(ax, xp[2][2:end-1], umean)
    ax = Axis(
        fig[1, 2];
        title = "Streamwise energy spectrum",
        xscale = log10,
        yscale = log10,
        xlabel = L"k_x",
        ylabel = L"\hat{U}_{cl} / U_0",
    )
    # ylims!(ax, (10^(0.0), 10^4.0))
    ksub = k[10:end]
    # lines!(ax, ksub, 1000 .* ksub .^ (-5 / 3); label = L"k^{-5/3}")
    lines!(ax, ksub, 1e7 .* ksub .^ -3; label = L"k^{-3}")
    scatter!(ax, k, E₀; label = "y = $y₀")
    scatter!(ax, k, E₁; label = "y = $y₁")
    axislegend(ax; position = :lb)
    # on(_ -> autolimits!(ax), E₁)

    fig
end

Solve unsteady problem

julia
state, outputs = solve_unsteady(;
    setup,
    ustart,
    tlims = (T(0), T(1)),
    method = RKMethods.RK44P2(),
    Δt = 0.001,
    processors = (
        rtp = realtimeplotter(;
            setup,
            # plot = fieldplot,
            # plot = energy_history_plot,
            # plot = energy_spectrum_plot,
            plot = meanplot,
            nupdate = 1,
        ),
        # anim = animator(; setup, path = "$outdir/vorticity.mkv", nupdate = 4),
        # vtk = vtk_writer(; setup, nupdate = 10, dir = outdir, filename = "solution"),
        # field = fieldsaver(; setup, nupdate = 10),
        log = timelogger(; nupdate = 100),
    ),
);
nothing #hide

Post-process ​

We may visualize or export the computed fields (u, p)

julia
outputs.rtp

Export to VTK

julia
save_vtk(state; setup, filename = joinpath(outdir, "solution"))

Plot pressure

julia
fieldplot(state; setup, fieldname = :pressure)

Plot initial velocity

julia
fieldplot((; u = u₀, p = p₀, t = T(0)); setup, fieldname = :velocitynorm)

Plot final velocity

julia
fieldplot(state; setup, fieldname = :velocitynorm)

Plot vorticity

julia
fieldplot(state; setup, fieldname = :vorticity)

Copy-pasteable code ​

Below is the full code for this example stripped of comments and output.

julia
using FFTW
using GLMakie
using IncompressibleNavierStokes
using LaTeXStrings

outdir = joinpath(@__DIR__, "output", "PlaneJets2D")

T = Float64

backend = CPU()
# using CUDA; backend = CUDABackend()

Re = T(6_000)

# V() = sqrt(T(467.4))
V() = T(21.619435700313733)

U_A(y) = V() / 2 * (tanh((y + T(0.5)) / T(0.1)) - tanh((y - T(0.5)) / T(0.1)))

U_B(y) =
    V() / 2 * (tanh((y + 1 + T(0.5)) / T(0.1)) - tanh((y + 1 - T(0.5)) / T(0.1))) +
    V() / 2 * (tanh((y - 1 + T(0.5)) / T(0.1)) - tanh((y - 1 - T(0.5)) / T(0.1)))

U_C(y) =
    V() / 2 * (
        tanh(((y + T(1.0)) / 1 + T(0.5)) / T(0.1)) -
        tanh(((y + T(1.0)) / 1 - T(0.5)) / T(0.1))
    ) +
    V() / 4 * (
        tanh(((y - T(1.5)) / 2 + T(0.5)) / T(0.2)) -
        tanh(((y - T(1.5)) / 2 - T(0.5)) / T(0.2))
    )

U_D(y) =
    V() / 2 * (
        tanh(((y + T(1.0)) / 1 + T(0.5)) / T(0.1)) -
        tanh(((y + T(1.0)) / 1 - T(0.5)) / T(0.1))
    ) -
    V() / 4 * (
        tanh(((y - T(1.5)) / 2 + T(0.5)) / T(0.2)) -
        tanh(((y - T(1.5)) / 2 - T(0.5)) / T(0.2))
    )

# U(y) = U_A(y)
# U(y) = U_B(y)
U(y) = U_C(y)
# U(y) = U_D(y)

U(x, y) = (1 + T(0.1) * (rand(T) - T(0.5))) * U(y)

# boundary_conditions = (
#     (PeriodicBC(), PeriodicBC()),
#     (PressureBC(), PressureBC())
# )

n = 64
# n = 128
# n = 256
x = LinRange(T(0), T(16), 4n + 1), LinRange(-T(10), T(10), 5n + 1)
plotgrid(x...)

setup = Setup(x, Re, backend);
# setup = Setup(; x, Re, boundary_conditions, backend);

ustart = velocityfield(setup, (dim, x, y) -> dim == 1 ? U(x, y) : zero(x));

function meanplot(state; setup)
    (; xp, Iu, Ip, Nu, N) = setup.grid

    umean = lift(state) do (; u, p, t)
        reshape(sum(u[1][Iu[1]]; dims = 1), :) ./ Nu[1][1] ./ V()
    end

    K = Nu[1][2] ÷ 2
    k = 1:(K-1)

    # Find energy spectrum where y = 0
    n₀ = findmin(abs, xp[2])[2]
    E₀ = lift(state) do (; u, p, t)
        u_y = u[1][:, n₀]
        abs.(fft(u_y .^ 2))[k.+1]
    end
    y₀ = xp[2][n₀]

    # Find energy spectrum where y = 1
    n₁ = findmin(y -> abs(y - 1), xp[2])[2]
    E₁ = lift(state) do (; u, p, t)
        u_y = u[1][:, n₁]
        abs.(fft(u_y .^ 2))[k.+1]
    end
    y₁ = xp[2][n₁]

    fig = Figure()
    ax = Axis(
        fig[1, 1];
        title = "Mean streamwise flow",
        xlabel = "y",
        ylabel = L"\langle u \rangle / U_0",
    )
    lines!(ax, xp[2][2:end-1], umean)
    ax = Axis(
        fig[1, 2];
        title = "Streamwise energy spectrum",
        xscale = log10,
        yscale = log10,
        xlabel = L"k_x",
        ylabel = L"\hat{U}_{cl} / U_0",
    )
    # ylims!(ax, (10^(0.0), 10^4.0))
    ksub = k[10:end]
    # lines!(ax, ksub, 1000 .* ksub .^ (-5 / 3); label = L"k^{-5/3}")
    lines!(ax, ksub, 1e7 .* ksub .^ -3; label = L"k^{-3}")
    scatter!(ax, k, E₀; label = "y = $y₀")
    scatter!(ax, k, E₁; label = "y = $y₁")
    axislegend(ax; position = :lb)
    # on(_ -> autolimits!(ax), E₁)

    fig
end

state, outputs = solve_unsteady(;
    setup,
    ustart,
    tlims = (T(0), T(1)),
    method = RKMethods.RK44P2(),
    Δt = 0.001,
    processors = (
        rtp = realtimeplotter(;
            setup,
            # plot = fieldplot,
            # plot = energy_history_plot,
            # plot = energy_spectrum_plot,
            plot = meanplot,
            nupdate = 1,
        ),
        # anim = animator(; setup, path = "$outdir/vorticity.mkv", nupdate = 4),
        # vtk = vtk_writer(; setup, nupdate = 10, dir = outdir, filename = "solution"),
        # field = fieldsaver(; setup, nupdate = 10),
        log = timelogger(; nupdate = 100),
    ),
);

outputs.rtp

save_vtk(state; setup, filename = joinpath(outdir, "solution"))

fieldplot(state; setup, fieldname = :pressure)

fieldplot((; u = u₀, p = p₀, t = T(0)); setup, fieldname = :velocitynorm)

fieldplot(state; setup, fieldname = :velocitynorm)

fieldplot(state; setup, fieldname = :vorticity)

This page was generated using Literate.jl.