Rayleigh-Taylor instability in 2D

Two fluids with different temperatures start mixing.

using CairoMakie
using IncompressibleNavierStokes

Output directory for saving results

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

"/home/runner/work/IncompressibleNavierStokes.jl/IncompressibleNavierStokes.jl/docs/build/examples/generated/output/RayleighTaylor2D"

Precision

T = Float64

Float64

Grid

n = 50
x = tanh_grid(T(0), T(1), n, T(1.5)), tanh_grid(T(0), T(2), 2n, T(1.5))
plotgrid(x...; figure = (; size = (300, 600)))

Setup

setup = Setup(;
    x,
    boundary_conditions = (;
        u = ((DirichletBC(), DirichletBC()), (DirichletBC(), DirichletBC())),
        temp = ((SymmetricBC(), SymmetricBC()), (SymmetricBC(), SymmetricBC())),
    ),
);

Initial conditions

start = (;
    u = velocityfield(setup, (dim, x, y) -> zero(x)),
    temp = temperaturefield(setup, (x, y) -> one(x) * (1 + sinpi(x) / 50 > y)),
)

(u = [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], temp = [1.0 1.0 … 0.0 0.0; 1.0 1.0 … 0.0 0.0; … ; 1.0 1.0 … 0.0 0.0; 1.0 1.0 … 0.0 0.0])

Solve equation

state, outputs = solve_unsteady(;
    force! = boussinesq!, # Solve the Boussinesq equations
    setup,
    start,
    tlims = (T(0), T(10)),
    params = (;
        viscosity = T(1e-3),
        gravity = T(1.0),
        gdir = 2, # Gravity acts in the y-direction
        conductivity = T(1e-3),
        dodissipation = true,
    ),
    processors = (;
        rtp = realtimeplotter(;
            setup,
            nupdate = 20,
            fieldname = :temperature,
            size = (400, 600),
        ),
        log = timelogger(; nupdate = 200),
    ),
);

[ Info: t = 1.70454	Δt = 0.0085	umax = 0.24	itertime = 0.0095
[ Info: t = 3.40909	Δt = 0.0085	umax = 0.8	itertime = 0.0019
[ Info: t = 5.11363	Δt = 0.0085	umax = 0.82	itertime = 0.0019
[ Info: t = 6.81818	Δt = 0.0085	umax = 0.95	itertime = 0.0019
[ Info: t = 8.52272	Δt = 0.0085	umax = 1	itertime = 0.0019
[ Info: Finished after 1174 time steps and 3.8 seconds

Copy-pasteable code

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

using WGLMakie
using IncompressibleNavierStokes

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

T = Float64

n = 50
x = tanh_grid(T(0), T(1), n, T(1.5)), tanh_grid(T(0), T(2), 2n, T(1.5))
plotgrid(x...; figure = (; size = (300, 600)))

setup = Setup(;
    x,
    boundary_conditions = (;
        u = ((DirichletBC(), DirichletBC()), (DirichletBC(), DirichletBC())),
        temp = ((SymmetricBC(), SymmetricBC()), (SymmetricBC(), SymmetricBC())),
    ),
);

start = (;
    u = velocityfield(setup, (dim, x, y) -> zero(x)),
    temp = temperaturefield(setup, (x, y) -> one(x) * (1 + sinpi(x) / 50 > y)),
)

state, outputs = solve_unsteady(;
    force! = boussinesq!, # Solve the Boussinesq equations
    setup,
    start,
    tlims = (T(0), T(10)),
    params = (;
        viscosity = T(1e-3),
        gravity = T(1.0),
        gdir = 2, # Gravity acts in the y-direction
        conductivity = T(1e-3),
        dodissipation = true,
    ),
    processors = (;
        rtp = realtimeplotter(;
            setup,
            nupdate = 20,
            fieldname = :temperature,
            size = (400, 600),
        ),
        log = timelogger(; nupdate = 200),
    ),
);

---

This page was generated using Literate.jl.