No-Slip Boundary Condition

[parfenyev]

Hello everyone!

I'm studing an example with 2d turbulence and I would like to modify it to study a domain with no-slip boundaries. So, I have changed the model:

grid = RectilinearGrid(size=(128, 128), extent=(2π, 2π),
                              topology=(Bounded, Bounded, Flat))

u_bcs = FieldBoundaryConditions(south = ValueBoundaryCondition(0.0),
                                north = ValueBoundaryCondition(0.0));

v_bcs = FieldBoundaryConditions(east = ValueBoundaryCondition(0.0),
                                west = ValueBoundaryCondition(0.0));

model = NonhydrostaticModel(timestepper = :RungeKutta3,
                              advection = UpwindBiasedFifthOrder(),
                                   grid = grid,
                    boundary_conditions = (u=u_bcs, v=v_bcs),
                               buoyancy = nothing,
                                tracers = nothing,
                                closure = ScalarDiffusivity(ν=1e-5)
                           )

After the simulation is complete (the rest of the code has not changed) I see that, e.g., v-velocity is non-zero at two boundaries:

interior(v_timeseries[10], :, :, 1)

128×129 view(::Array{Float64, 4}, :, :, 1, 10) with eltype Float64:
 0.0   0.00891897   0.0438791   …   0.0144025    0.00409626  0.0
 0.0   0.0293911    0.0507547       0.0219927    0.00980964  0.0
 0.0   0.0250097    0.0341092       0.0359548    0.00680388  0.0
 0.0   0.0117144    0.0101972       0.0441318    0.0130078   0.0
 0.0  -0.00261697  -0.00331914      0.0386366    0.0172396   0.0
 0.0  -0.0150908   -0.0139833   …   0.0201599    0.0125234   0.0
 0.0  -0.0254312   -0.032287       -0.00531311   0.00124912  0.0
 0.0  -0.0241077   -0.0249426      -0.0264362   -0.00884336  0.0
 0.0  -0.0137805   -0.00727661     -0.0380925   -0.0174468   0.0
 0.0  -0.00191326   0.0165824      -0.0429263   -0.0255871   0.0
 0.0   0.021664     0.0471763   …  -0.0368307   -0.0275779   0.0
 0.0   0.0477113    0.0693482      -0.0102541   -0.00603961  0.0
 0.0   0.0482472    0.0720176       0.02151      0.0174318   0.0
 ⋮                              ⋱                            
 0.0  -0.0425603   -0.0817042      -0.00649985  -0.00137887  0.0
 0.0  -0.048267    -0.091053        0.011604     0.00588816  0.0
 0.0  -0.0502379   -0.0903577       0.0232738    0.0119298   0.0
 0.0  -0.0415356   -0.072126        0.0330497    0.017806    0.0
 0.0  -0.0233295   -0.0439425   …   0.03966      0.0200196   0.0
 0.0  -0.00403236  -0.0107361       0.0390461    0.017641    0.0
 0.0   0.0195767    0.02647         0.0291041    0.0135361   0.0
 0.0   0.0354432    0.0554998       0.0169403    0.0120396   0.0
 0.0   0.0445907    0.0726527       0.00732206   0.0149755   0.0
 0.0   0.0406634    0.0702592   …   0.00780794   0.0214572   0.0
 0.0  -0.0356932    0.0191723       0.0184108    0.0267699   0.0
 0.0   0.0171275   -0.0098884       0.00349168  -0.0209539   0.0

How can I set no-slip (zero velocity) boundary conditions?

[glwagner]

Hi @parfenyev ! Thanks for opening a discussion!

In short: what you're seeing is expected!

Oceananigans uses a finite volume method, which means that the numbers you see represent cell-averaged values for the field. Finite volumes solve equations in the "weak form", which means they time-step a volume-integrated form of the equations. They differ from "strong form" methods like finite difference or typical pseudospectral / global spectral methods in which the discretized field represents a continuous field evaluated at particular points. (One reason why heatmap is in fact a really nice way to visualize finite volume fields!)

In this case, you're printing values for y-velocity component v. At i=1, those values represent the average value of the y-velocity component in the cell that's adjacent to the wall. The y-velocity is indeed non-zero, because there's still flow close the wall. The value of the y-velocity on the wall is not represented in the code; instead the boundary conditions are enforced by adding a momentum flux contribution (which slows down the flow) to those boundary-adjacent cells, where the flow is non-zero.

So while you can't verify that the no-slip condition is enforced by looking for "0" somewhere in your field, you can verify that the effect of the no-slip condition is felt (ie, momentum should be extracted). One thing you can try is to visualize two cases: one with free-slip conditions (the default), and one with no-slip conditions (what you've written), and see if you can detect the difference. You also likely want to increase the diffusivity from 1e-5. For a resolution 128^2 with no-slip / viscous boundary layers, maybe something like 1e-3 or 1e-2 is more reasonable?

You can also try the initial condition

# streamfunction psi(x, y) = sin(x) * sin(y)

ui(x, y, z) =    sin(x) * cos(y)
vi(x, y, z) = - cos(x) * sin(y)

set!(model, u=ui, v=vi)

It'd be great to see the results too if you want to post them in this discussion!

[parfenyev]

Hi @glwagner ! Thank you for the answer.

To check the no-slip conditions, I simulated the flow in the channel in laminar regime. Below is my code:

using Oceananigans
​
grid = RectilinearGrid(size=(4, 128), extent=(2π, 1),
                              topology=(Periodic, Bounded, Flat))
​
# no-slip BC
u_bcs = FieldBoundaryConditions(south = ValueBoundaryCondition(0.0),
                                north = ValueBoundaryCondition(0.0));
​
# external force
F = 1.0
u_forcing(x, y, z, t) = F
v_forcing(x, y, z, t) = 0
​
# viscosity
nu = 0.1
​
model = NonhydrostaticModel(timestepper = :RungeKutta3,
                              advection = UpwindBiasedFifthOrder(),
                                   grid = grid,
                    boundary_conditions = (u=u_bcs,),
                                forcing = (u=u_forcing, v=v_forcing),
                               buoyancy = nothing,
                                tracers = nothing,
                                closure = ScalarDiffusivity(ν=nu))
​
set!(model, u=0, v=0)
simulation = Simulation(model, Δt=1e-4, stop_time=10.0)
​
u, v, w = model.velocities
simulation.output_writers[:fields] = JLD2OutputWriter(model, (; u, v),
                                                      schedule = TimeInterval(0.02),
                                                      filename = "two_dimensional_turbulence.jld2",
                                                      overwrite_existing = true)

run!(simulation)

using Plots

u_timeseries = FieldTimeSeries("two_dimensional_turbulence.jld2", "u")
xu, yu, zu = nodes(u_timeseries)

anim = @animate for i in 1:length(u_timeseries.times)

    ui = interior(u_timeseries[i], :, :, 1)

    u_lim = 1.25
    u_levels = range(-u_lim, stop=u_lim, length=50)

    kwargs = (xlabel="x", ylabel="y", aspectratio=3, linewidth=0, colorbar=true,
              xlims=(0, model.grid.Lx), ylims=(0, model.grid.Ly))

    u_plot = contourf(xu, yu, clamp.(ui', -u_lim, u_lim);
                       color = :balance,
                      levels = u_levels,
                       clims = (-u_lim, u_lim),
                      kwargs...)

    plot(u_plot, title=["Velocity"], layout=(1, 1), size=(600, 500))
end

mp4(anim, "video.mp4", fps=18)

and here is comparison with the theoretical Poiseuille profile:

# last frame and final velocity
last = length(u_timeseries.times)
vel = interior(u_timeseries[last], 1, :, 1)

# comparison with Poiseuille profile
plot(yu, vel, lw=2, xlim=(0,1), ylim=(0,1.3), framestyle=:box, label="numeric",
    gridstyle=:dash, xlabel = "y", ylabel = "velocity", thickness_scaling=1.5)
plot!(yu, yu->F/(2*nu)*yu*(1.0-yu), ls=:dash, lw=2, label="theory")

[glwagner]

theory and numerics seem to be in close agreement!