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
u_n.assign(u_)
p_n.assign(p_)
#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)
print(u_n(x,y))