FlatExtrapolationOpenBoundaryCondition in 2D natural convection simulation

[Dixi9]

Hi, I'm new to Oceananigans and currently running a 2D simulation of natural convection in a channel. The top boundary has a constant heat flux, the bottom is adiabatic, and I want open boundaries on the left and right. For this, I use FlatExtrapolationOpenBoundaryCondition(nothing) for west and east velocity BCs and GradientBoundaryCondition(0.0) for east buoyancy BCs. Initially, the temperature field diffuses as expected, but after some time (just few seconds), the velocity field suddenly surges across the domain (like wind goes through from west boundary), destroying the existing flow structure and collapsing it into a laminar pattern. Is this a known issue when using open boundaries, or am I missing something in the setup? I’ve tried many different initial conditions, and also smaller heat fluxes, but the issue still happens.
Thank you in advance for your help and suggestions!

Here is my code:

    u_bcs = FieldBoundaryConditions(
    top=ValueBoundaryCondition(0.0),
    bottom=ValueBoundaryCondition(0.0),
    west=FlatExtrapolationOpenBoundaryCondition(nothing),
    east=FlatExtrapolationOpenBoundaryCondition(nothing)
)

w_bcs = FieldBoundaryConditions(
    west=ValueBoundaryCondition(0.0),
    east=GradientBoundaryCondition(0.0)
)

@inline function top_flux(x, t)
    if 2.5 ≤ x ≤ 7.5
        return 1
    else
        return 0.0
    end
end

b_top_bcs = GradientBoundaryCondition(top_flux)

buoyancy_bc = FieldBoundaryConditions(
    top=b_top_bcs,
    bottom=GradientBoundaryCondition(0.0),
    west=ValueBoundaryCondition(0),
    east=GradientBoundaryCondition(0)
)

g̃ = (-1.0, 0.0, 0.0)

model = NonhydrostaticModel(; grid,
    boundary_conditions=(u=u_bcs, w=w_bcs, b=buoyancy_bc),
    # advection=UpwindBiased(order=5),
    advection=WENO(),
    timestepper=:RungeKutta3,
    tracers=:b,
    forcing=(u = u_forcing,),
    buoyancy=BuoyancyForce(BuoyancyTracer(), gravity_unit_vector=g̃),
    closure=ScalarDiffusivity(; ν, κ)
)

U-Velocity.mp4

Temperature.mp4

The issue occurs at Ra = 5e10 even with uniform heating at the top.

[Dixi9]

I’m trying to apply zero-gradient velocity boundary conditions on both west and east boundaries with no background velocity. But if I replace it with FlatExtrapolationOpenBoundaryCondition(0.1,relaxation_timescale=1) it seems to work. Is FEOBC(nothing) not suitable for both velocity boundaries now?

[jagoosw]

Hi @Dixi9,

Thanks for finding this issue. Even though FEOBC tries to make zero gradient across the boundary the interaction with the pressure solver does seem to cause it to suddenly create a veloicty when it goes non-zero on the boundary.

You can set a relaxation timescale without setting a boundary velocity like:

FlatExtrapolationOpenBoundaryCondition(; relaxation_timescale = 0.01)

It might also be worth making the domain wider relative to the forcing region to make sure that the obc isn't effecting the interesting dynamics.

[glwagner]

Should you set val=0 then? @jagoosw

[jagoosw]

Yeah if its set to zero then PAOBC is really just relaxing to zero so that might work

[Dixi9]

Hi @glwagner @jagoosw ,
Thanks for the suggestions. I tried PAOBC(val=0), but it pins the east/west velocity to 0, whereas I need a zero-gradient outflow/inflow. FEOBC(nothing) seems closer, but when turbulence reaches the east boundary I get a large spike in u and the solution was incorrect.

[jagoosw]

Why is a zero gradient required? Its not numerically possible to set zero gradient on a wall normal velocity as there isn't sufficient information. For example even if you directly set the velocity point on the wall to the be identical to the first interior point, then the gradient at the cell center is zero but we still don't know what the gradient at the face is since we have no information about the next point out of the domain. FEOBC attempts to approximate zero gradient but as you've noted it doesn't work very well.

[glwagner]

Would it work to relax to, say, the area-averaged velocity over the open boundary?