Example 7: Rayleigh-Bénard Instability

missing

This test is similar to Rayleigh-Taylor but even more amazing. The top side of the domain is hot and the bottom is cold. The hot fluid expands and rises up, then cools again and sinks. The resulting dynamical pattern is called a convective cell. It is somehow related to atmosphere and weather. You can read more on wikipedia.

module heat
include("../src/LagrangianVoronoi.jl")
using .LagrangianVoronoi


const g = 9.8 #gravitation acceleration
const mu = 1e-4 #dynamic viscosity
const H = 0.1 #height
const W = 0.1
const gamma = 1.4 #adiabatic index
const rho0 = 10.0

const cV = 10.0  #sp. heat capacity
const R = (gamma - 1.0)*cV
const thermal_k = mu*cV*gamma/0.71 #thermal conductivity

const P0 = 1e4
const c0 = sqrt(gamma*P0/rho0)
const Tu = P0/(rho0*R) #temperature of upper boundary (cooler)
const Td = 1000.0 #30.0 #Tu*(1.0 + contrast) #temperature of lower boundary (heater)

const export_path = "results/heat"
const dr = H/100
const v_char = 2.0
const dt = 0.1*dr/v_char
const nframes = 400
const t_end = 1.0

function print_info()
    thermal_D = thermal_k/(rho0*gamma*cV) #thermal diffusivity
    thermal_a = 1.0/Tu
    nu = mu/rho0
    @show Pr = nu/thermal_D # Prandtl number
    @show Ra = thermal_a*g*(Td-Tu)*(H^3)/(nu*thermal_D) # Rayleigh number
end

Define the initial state (exponential atmosphere) and the boundary condition for temperature.

function exp_atmo!(p::VoronoiPolygon)
    p.T = Tu
    p.P = P0*exp(-g*p.x[2]/(R*Tu))
    p.e = cV*p.T
    p.rho = p.P/(R*p.T)
    p.mass = p.rho*area(p)
    p.k = thermal_k
    p.cV = cV
end

function T_bc(::RealVector, bdary::Int)::Float64
    if (bdary == BDARY_DOWN)
        return Td
    elseif (bdary == BDARY_UP)
        return Tu
    end
    return NaN # NaN indicates an adiabatic wall
end

The rest is fairly standard, except we also need to include Fourier diffusion in our time-marching scheme.

mutable struct Simulation <: SimulationWorkspace
    grid::GridNSF
    solver::PressureSolver{PolygonNSF}
    E::Float64
    S::Float64
    E_kinetic::Float64
    Simulation() = begin
        xlims = (0.0, W)
        ylims = (0.0, H)
        domain = Rectangle(xlims = xlims, ylims = ylims)
        grid = GridNSF(domain, dr, xperiodic = true)
        populate_lloyd!(grid, ic! = exp_atmo!)
        return new(grid, PressureSolver(grid), 0.0, 0.0, 0.0)
    end
end

function step!(sim::Simulation, t::Float64)
    move!(sim.grid, dt)
    gravity_step!(sim.grid, -g*VECY, dt)
    ideal_eos!(sim.grid, gamma)
    find_pressure!(sim.solver, dt)
    pressure_step!(sim.grid, dt)
    ideal_temperature!(sim.grid) # find the temperature according to ideal gas laws
    fourier_step!(sim.grid, dt)  # update the temperature by Fourier diffusion
    heat_from_bdary!(sim.grid, dt, T_bc, gamma) # add heat flux from boundaries
    find_D!(sim.grid)
    viscous_step!(sim.grid, dt)
    find_dv!(sim.grid, dt)
    relaxation_step!(sim.grid, dt)
    return
end

function postproc!(sim::Simulation, t::Float64)
    sim.S = 0.0
    sim.E = 0.0
    sim.E_kinetic = 0.0
    P0 = rho0*c0^2/gamma
    for p in sim.grid.polygons
        sim.S += p.mass*p.cV*(log(abs(p.P/P0)) - gamma*log(abs(p.rho/rho0)))
        sim.E += p.mass*p.e
        sim.E_kinetic += 0.5*p.mass*norm_squared(p.v)
    end
    percent = round(100*t/t_end, digits = 5)
    println("t = $t ($(percent)%)")
    @show sim.E
    @show sim.S
    @show sim.E_kinetic
    println()
end

function main()
    print_info()
    sim = Simulation()
    run!(sim, dt, t_end, step!; path = export_path,
        vtp_vars = (:v, :P, :T, :rho),
        postproc! = postproc!,
        nframes = nframes
    )
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end

end

This page was generated using Literate.jl.