Why do Dirichlet BCs seem not to be satisfied exactly?

Hello everyone!

I am interested in Navier-Stokes-kind simulations with Fenics, but I seem not to understand how do the Dirichlet boundary conditions work. Could anyone please guide me to the right understanding?

I attach an example, entirely based on Fenics tutorial program ft07_navier_stokes_channel.py. I only fixed deprecated attribute “.array()”, deleted picture drawing and added some text output to show what I mean.

So, there are no-slip conditions on the walls at y=0 and y=1. But when I visualize solution with paraview or evaluate function in Fenics (like in the text output I added in the end of the program) I get essentially non-zero values, especially on first time steps. Why is that so? It looks like some extrapolation, so maybe Fenics kind of “forget” that there are no basis functions on the boundary or smth like that?

I use Fenics 2019.1.0 on Ubuntu 18.04, have installed it from Ubuntu repository.

I don’t see, how I can add a file, so I will post the entire code here. Hope that is not too bad…

from __future__ import print_function
from fenics import *
import numpy as np

T = 10.0           # final time
num_steps = 500    # number of time steps
dt = T / num_steps # time step size
mu = 1             # kinematic viscosity
rho = 1            # density

# Create mesh and define function spaces
mesh = UnitSquareMesh(16, 16)
V = VectorFunctionSpace(mesh, 'P', 1)
Q = FunctionSpace(mesh, 'P', 1)

# Define boundaries
inflow  = 'near(x[0], 0)'
outflow = 'near(x[0], 1)'
walls   = 'near(x[1], 0) || near(x[1], 1)'

# Define boundary conditions
bcu_noslip  = DirichletBC(V, Constant((0, 0)), walls)
bcp_inflow  = DirichletBC(Q, Constant(8), inflow)
bcp_outflow = DirichletBC(Q, Constant(0), outflow)
bcu = [bcu_noslip]
bcp = [bcp_inflow, bcp_outflow]

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)

# Define functions for solutions at previous and current time steps
u_n = Function(V)
u_  = Function(V)
p_n = Function(Q)
p_  = Function(Q)

# Define expressions used in variational forms
U   = 0.5*(u_n + u)
n   = FacetNormal(mesh)
f   = Constant((0, 0))
k   = Constant(dt)
mu  = Constant(mu)
rho = Constant(rho)

# Define strain-rate tensor
def epsilon(u):
    return sym(nabla_grad(u))

# Define stress tensor
def sigma(u, p):
    return 2*mu*epsilon(u) - p*Identity(len(u))

# Define variational problem for step 1
F1 = rho*dot((u - u_n) / k, v)*dx + \
     rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \
   + inner(sigma(U, p_n), epsilon(v))*dx \
   + dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds \
   - dot(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)

# Define variational problem for step 2
a2 = dot(nabla_grad(p), nabla_grad(q))*dx
L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx

# Define variational problem for step 3
a3 = dot(u, v)*dx
L3 = dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx

# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)

# Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]

# Time-stepping
t = 0
for n in range(num_steps):
    print("---------------- step ",n," ---------------------------")

    # Update current time
    t += dt

    # Step 1: Tentative velocity step
    b1 = assemble(L1)
    [bc.apply(b1) for bc in bcu]
    solve(A1, u_.vector(), b1)

    # Step 2: Pressure correction step
    b2 = assemble(L2)
    [bc.apply(b2) for bc in bcp]
    solve(A2, p_.vector(), b2)

    # Step 3: Velocity correction step
    b3 = assemble(L3)
    solve(A3, u_.vector(), b3)

    # Compute error
    u_e = Expression(('4*x[1]*(1.0 - x[1])', '0'), degree=2)
    u_e = interpolate(u_e, V)
    error = np.abs(u_e.vector().get_local() - u_.vector().get_local()).max()
    print('t = %.2f: error = %.3g' % (t, error))
    print('max u:', u_.vector().get_local().max())

    # Update previous solution
    #check Dirichlet bc 
    for v in vertices(mesh):
        x = v.point().x()
        y = v.point().y()
        if (y == 0)or(y == 1):
            print(x, y)

You can see that no BC is applied to the velocity during the velocity correction, “Step 3”, which explains the nonzero values. In the main documented Navier–Stokes demo, you can see that BCs are applied during the velocity correction step:


I see that BCs are not applied in the tutorial that you referenced, but I don’t work with projection methods, so I can’t say whether or not there is some good reason to omit the BC during the correction step.

1 Like

I’ve been doing some work with the IPCS method, and at least as I’ve seen, one Tends to impose boundary conditions on the tentative velocity step and the pressure correction step. The velocity correction step is treated as a projection. I cant say i’ve seen this issue being debated, but there is a wide literature and discussion on the appropriate boundary conditions for projection scheme in general. I’ll see if i can dig up some papers over the weekend for reference.

1 Like

Thank you very much, kamensky and dokken, for the answers. That’s exactly what I needed in terms of “guidance to the right understanding”. Will investigate it further.

P.S. dokken, I will be very grateful for any links or papers, but - only if you will have time. Anyway Google is in my direct disposal :slight_smile:

Here are some papers that discusses boundary conditions for splitting schemes:

  • A. Vreman, The projection method for the incompressible Navier–Stokes equations:
    The pressure near a no-slip wall, Journal of Computational Physics 263 (2014) 353–374. doi:10.1016/j.jcp.2014.01.035.
  • P. M. Gresho, R. L. Sani, On pressure boundary conditions for the incompressible
    Navier-Stokes equations, International Journal for Numerical Methods in Fluids 7 (10)
    (1987) 1111–1145. doi:10.1002/fld.1650071008.
  • R. L. Sani, J. Shen, O. Pironneau, P. M. Gresho, Pressure boundary condition for
    the time-dependent incompressible Navier–Stokes equations, International Journal for
    Numerical Methods in Fluids 50 (6) (2006) 673–682. doi:10.1002/fld.1062

Hope this can help you.

1 Like

Once again - thank you very much, dokken, for your time!