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

Kolmogorov flow (2D)

The Kolmogorov flow in a periodic box $\mathrm{\Omega }=[0,1{]}^2$ is initiated via the force field

$$f(x,y)=\left(\begin{array}{c}\sin (\pi ky)\\ 0\end{array}\right)$$

where k is the wavenumber where energy is injected.

Packages

We just need the IncompressibleNavierStokes and a Makie plotting package.

using CairoMakie
using IncompressibleNavierStokes

Setup

Define a uniform grid with a steady body force field.

n = 256
axis = range(0.0, 1.0, n + 1)
setup = Setup(;
    x = (axis, axis),
    Re = 2e3,
    bodyforce = (dim, x, y, t) -> (dim == 1) * 5 * sinpi(8 * y),
    issteadybodyforce = true,
);
ustart = random_field(setup, 0.0; A = 1e-2);
nothing #hide

Plot body force

Since the force is steady, it is just stored as a field.

heatmap(setup.bodyforce[:, :, 1])

Solve unsteady problem

state, outputs = solve_unsteady(;
    setup,
    ustart,
    tlims = (0.0, 2.0),
    Δt = 1e-3,
    processors = (
        rtp = realtimeplotter(; setup, nupdate = 100),
        ehist = realtimeplotter(;
            setup,
            plot = energy_history_plot,
            nupdate = 10,
            displayfig = false,
        ),
        espec = realtimeplotter(;
            setup,
            plot = energy_spectrum_plot,
            nupdate = 10,
            displayfig = false,
        ),
        log = timelogger(; nupdate = 100),
    ),
);
nothing #hide

Field plot

outputs.rtp

Energy history

outputs.ehist

Energy spectrum

outputs.espec

Copy-pasteable code

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

using GLMakie
using IncompressibleNavierStokes

n = 256
axis = range(0.0, 1.0, n + 1)
setup = Setup(;
    x = (axis, axis),
    Re = 2e3,
    bodyforce = (dim, x, y, t) -> (dim == 1) * 5 * sinpi(8 * y),
    issteadybodyforce = true,
);
ustart = random_field(setup, 0.0; A = 1e-2);

heatmap(setup.bodyforce[:, :, 1])

state, outputs = solve_unsteady(;
    setup,
    ustart,
    tlims = (0.0, 2.0),
    Δt = 1e-3,
    processors = (
        rtp = realtimeplotter(; setup, nupdate = 100),
        ehist = realtimeplotter(;
            setup,
            plot = energy_history_plot,
            nupdate = 10,
            displayfig = false,
        ),
        espec = realtimeplotter(;
            setup,
            plot = energy_spectrum_plot,
            nupdate = 10,
            displayfig = false,
        ),
        log = timelogger(; nupdate = 100),
    ),
);

outputs.rtp

outputs.ehist

outputs.espec

---

This page was generated using Literate.jl.