Why is the normal stress non-zero at boundary where I impose zero normal stress?

variational formulation
1

Hello,
I am solving Stokes flow for a rectangular domain. The top of the domain is split into two halves, where the boundary condition on the right half should be zero normal stress. Even though my solution looks ok when I plot it, the normal stress I calculate through that boundary during post-processing is not zero. How can this be?
This is my minimal working example

from dolfin import *
from mshr import *
import matplotlib.pyplot as plt
import numpy as np
# Define mesh and geometry - We solve for half of the domain we need, and impose symmetry

mesh = RectangleMesh(Point(0, 0), Point(1, 1), 100, 100)
n = FacetNormal(mesh)

# Define Taylor--Hood function space W
V = VectorElement("CG", triangle, 2)
Q = FiniteElement("CG", triangle, 1)
W = FunctionSpace(mesh, MixedElement([V, Q]))

# Define Function and TestFunction(s)
w = Function(W)
(u, p) = split(w)
(v, q) = split(TestFunction(W))

# Define the viscosity and bcs

mu = Constant(1.0)
# mu = Expression('exp(-a*pow(x[0],2))', degree=2, a=10)

u_in = Constant(-2.0)
u_c = Constant(-1.0)

inflow = 'near(x[1], 1.0) && x[0]<=0.5'
wall = 'near(x[0], 1)'
centre = 'near(x[0], 0.0)'
outflow = 'near(x[1], 0.0)'
bcu_inflow = DirichletBC(W.sub(0), (0.0, u_in), inflow)
bcu_wall = DirichletBC(W.sub(0), (0.0, u_c), wall)
bcu_outflow = DirichletBC(W.sub(0), (0.0, u_c), outflow)
bcu_symmetry = DirichletBC(W.sub(0).sub(0), Constant(0.0), centre)
bcs = [bcu_inflow, bcu_wall, bcu_outflow, bcu_symmetry]
# Define stress tensor

x = SpatialCoordinate(mesh)


def epsilon(v):
    return sym(as_tensor([[v[0].dx(0), v[0].dx(1), 0],
                          [v[1].dx(0), v[1].dx(1), 0],
                          [0, 0, 0]]))


# stress tensor
def sigma(v, p):
    return 2*mu*epsilon(v)-Id(p)


def Id(p):
    return as_tensor([[p, 0, 0],
                      [0, p, 0],
                     [0, 0, p]])


def cond(v):
    return sym(as_tensor([[v[0].dx(0), v[0].dx(1)],
                          [v[1].dx(0), v[1].dx(1)]]))


def sigmabc(v, p):
    return 2*mu*cond(v) - p*Identity(2)


# Define the variational form
f = Constant((0, -1))

colors = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
colors.set_all(0)  # default to zero
# We match the colours to the defined sketch in the Fenics chapter
CompiledSubDomain("near(x[0], 0.0)").mark(colors, 4)
CompiledSubDomain("near(x[1], 1.0) && x[0]<=0.5").mark(colors, 0)
CompiledSubDomain("near(x[1], 1.0) && x[0]>=0.5").mark(colors, 1)
CompiledSubDomain("near(x[0], 1.0)").mark(colors, 2)  # wall
CompiledSubDomain("near(x[1], 0.0)").mark(colors, 3)  # outflow

# Create the measure
ds = Measure("ds", domain=mesh, subdomain_data=colors)

a1 = (inner(sigma(u, p), epsilon(v))) * dx
a2 = (- div(u) * q - dot(f, v)) * dx

b1 = - dot(dot(sigmabc(u, p), v), n) * ds(0)
b3 = - dot(dot(sigmabc(u, p), v), n) * ds(2)
b4 = - dot(dot(sigmabc(u, p), v), n) * ds(3)
F = a1 + a2

# Solve problem
solve(F == 0, w, bcs)

# Plot solutions
(u, p) = w.split()
# We check stress at the end

stress = (dot(as_vector([sigmabc(u, p)[0, 0], sigmabc(u, p)[1, 1]]), n) ) * ds(1) # stress calculation at ds(1) - this normal stress should be zero!
flux = dot(u, n) * dot(u, n) * ds(1)

total_stress_new = assemble(stress)
total_flux_new = assemble(flux)
print("Total stress on boundary ds1 is", total_stress_new)
print("Total flux on boundary ds1 is", total_flux_new)

Does anyone know why my normal stress is not zero? Or maybe a better way to impose zero normal stress at ds(1)?

2

Furthermore, I can calculate the stress explicitly (as I am probably doing it wrong the way it is in the code), as

stress = (sigmabc(u, p)[0, 0]*n[0] + sigmabc(u, p)[1, 1]*n[1])*ds(1)

which is also non-zero. But

stress = (sigmabc(u, p)[0, 0]*n[0] )*ds(1)

is zero. Now I’m really confused… The normal vector at the boundary ds(1) should be as_vector(n[0], 0), it’s just a vector pointing up, right? Then why does adding the component multiplied by n[1] make it non-zero? I know that for this particular example it doesn’t matter and I can just write it times n[0] as above, but my goal is to solve on a domain where the normal to the boundary does not align with the components of the velocity vectors (e.g. an inclined boundary on a Polygon domain), and so I can’t just make the normal vector (1,0). Any advice would be greatly appreciated. I have been stuck on this for weeks and have been unable to make progress on my actual problem, with the Polygon boundary.

3

Have you ever solved? I have the same problem!

4

The normal vector at ds(1) should be (0, 1), not (1, 0), as you are looking at the top boundary

So it is expected that the first term is 0, and that the second term gives a non-zero contribution