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.

using FFTW
using CairoMakie
using IncompressibleNavierStokes
using LaTeXStrings

Output directory

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

Floating point type

T = Float64

Array type

ArrayType = Array
# using CUDA; ArrayType = CuArray
# using AMDGPU; ArrayType = ROCArray
# using oneAPI; ArrayType = oneArray
# using Metal; ArrayType = MtlArray

Reynolds number

Re = T(6_000)

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

# 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

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

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

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

Initial conditions

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

Real time plot: Streamwise average and spectrum

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

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)

outputs.rtp

Export to VTK

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

Plot pressure

fieldplot(state; setup, fieldname = :pressure)

Plot initial velocity

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

Plot final velocity

fieldplot(state; setup, fieldname = :velocitynorm)

Plot vorticity

fieldplot(state; setup, fieldname = :vorticity)

---

This page was generated using Literate.jl.