Hello,
I am implementing Navier-Stokes equations in Fenics, following this example. The manifold \Omega is a rectangle, and \partial \Omega_{\rm IN} the left boundary of \Omega. Among the boundary conditions for the pressure increment, \phi, I impose
n^i \partial_i \phi = 0
on \partial \Omega_{\rm IN}. I implement these BCs by setting to zero the term
\langle n^i \partial_i \phi \, \nu_\phi \rangle_{\partial \Omega_{\rm IN}} = 0
in the variational functional, where n is the normal outwards vector and \nu_{\phi} the test function related to \phi.
However, if I then compute
\int _{\partial \Omega_{\rm IN}} n^i \partial_i \phi d S
I obtain a value which is significantly different from zero.
Here is a minimal working example (mesh files are here):
from __future__ import print_function
from fenics import *
from mshr import *
import meshio
import ufl as ufl
#CHANGE PARAMETERS HERE
tol = 1E-3
L = 2.2
h = 0.41
T = 0.001
num_steps = 4
#CHANGE PARAMETERS HERE
dt = T / num_steps # time step size
mu = 0.001
rho = 1.0
i, j, k = ufl.indices(3)
#create mesh
mesh=Mesh()
with XDMFFile("triangle_mesh.xdmf") as infile:
infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("line_mesh.xdmf") as infile:
infile.read(mvc, "name_to_read")
O = VectorFunctionSpace(mesh, 'P', 2, dim=2)
Q = FunctionSpace(mesh, 'P', 1)
v = Function(O)
v_ = Function(O)
sigma = Function(Q)
phi = Function(Q)
nu = TestFunction(O)
q = TestFunction(Q)
inflow = 'near(x[0], 0)'
outflow = 'near(x[0], 2.2)'
walls = 'near(x[1], 0) || near(x[1], 0.41)'
mf = dolfin.cpp.mesh.MeshFunctionSizet(mesh, mvc)
#test for surface elements
ds_l = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=2)
ds_r = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=3)
ds_t = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=4)
ds_b = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=5)
xdmffile_v = XDMFFile("v.xdmf")
xdmffile_sigma = XDMFFile( "sigma.xdmf")
v_n = Function(O)
v_ = Function(O)
sigma_n = Function(Q)
sigma_ = Function(Q)
inflow_profile_v = Expression(('4.0*1.5*x[1]*(0.41 - x[1]) / pow(0.41, 2)', '0'), degree=2, h=h)
# Define symmetric gradient
def epsilon(u):
# nabla_grad(u)_{i,j} = (u[j]).dx[i]
#sym(nabla_grad(u)) = nabla_grad(u)_{i,j} + nabla_grad(u)_{j,i}
# return sym(nabla_grad(u))
return as_tensor(0.5*(u[i].dx(j) + u[j].dx(i)), (i,j))
#the normal vector on the inflow and outflow
def n_inout():
x = ufl.SpatialCoordinate(mesh)
u = as_tensor([conditional(lt(x[0], L/2), -1.0, 1.0), 0.0] )
return as_tensor(u[k], (k))
# Define boundary conditions
bcv_inflow = DirichletBC(O, inflow_profile_v, inflow)
bcv_walls = DirichletBC(O, Constant((0, 0)), walls)
bc_phi_outflow = DirichletBC(Q, Constant(0), outflow)
bc_v = [bcv_inflow, bcv_walls]
bc_phi = [bc_phi_outflow]
V = 0.5 * (v_n + v_)
Deltat = Constant(dt)
mu = Constant(mu)
rho = Constant(rho)
# step 1
F1v = rho * ( (v_[i] - v_n[i]) / Deltat + v_n[j] *((v_n[i]).dx(j)) ) * nu[i] * dx \
+ ( 2.0 * mu * (epsilon(V))[i, j] * (nu[i]).dx(j) + sigma_n * (nu[i]).dx(i) ) * dx \
+ (- sigma_n * n_inout()[i] * nu[i] - mu * (V[j]).dx(i) * n_inout()[j] * nu[i]) * ds
F2 = ( (phi.dx(i)) * (q.dx(i)) + (rho / Deltat) * ((v_[i]).dx(i)) * q ) * dx
F3v = ((v_[i] - v[i]) - (Deltat / rho) * (phi.dx(i))) * nu[i] * dx
print("Starting time iteration ...", flush=True)
# Time-stepping
t = 0
for n in range(num_steps):
t += dt
solve(F1v == 0, v_, bc_v)
solve(F2 ==0, phi, bc_phi)
solve(F3v == 0, v)
# Update previous solution
v_n.assign(v)
sigma_n.assign(-phi + sigma_n)
# Write the solution to file
xdmffile_v.write(v_n, t)
xdmffile_sigma.write(sigma_n, t)
print("BC = ", assemble(((phi.dx(i)) * (n_inout())[i])**2 * ds_l))
print("\t%.2f %%" % (100.0 * (t / T)), flush=True)
print("... done.", flush=True)
The natural BC is imposed in F2 and the printout of the natural BC is in print("BC = ", assemble(((phi.dx(i)) * (n_inout())[i])**2 * ds_l)).
The output is
Starting time iteration ...
[...] BC = 5746230.869323557
25.00 %
[...]
BC = 3392726.0511582345
50.00 %
[...]
Why is this residual so large ? Thank you !