How to impose symmetry conditions on domain centreline, Stokes equations?

variational formulation
1

I am solving the Stokes flow equations a rectangular domain. As far as I can tell, my code works well. However, my domain is symmetric with respect to the vertical axis, and I want to impose symmetry instead of solving for both sides. Currently, my results show the exact same behaviour on the right and on the left, which is great, but in future I would like to look at a more complicated configuration and imposing symmetry would save time.

I want to impose zero tangential stress on the centerline for the flow problem.

The right and left walls of my domain move with velocity u_c, which is what drives the recirculation zones on the sides.

This is the code I am using:

from dolfin import *

#Define mesh and geometry
mesh = RectangleMesh(Point(-1, 0), Point(1, 1), 120, 60)

#Mesh that I want, half of my current domain
#mesh = RectangleMesh(Point(0, 0), Point(1, 1), 120, 60)
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 = Expression(‘exp(-a*pow(x[0],2))’, degree=2, a=5)
p0 = Expression(“1”, degree=1)
u_in = Constant(-2.0)
u_c = Constant(-1.0)
inflow = ‘near(x[1], 1.0) && x[0]<0.5 && x[0]>-0.5’
wall1 = ‘near(x[0], 1.0)’
wall2 = ‘near(x[0], -1.0)’
outflow = ‘near(x[1], 0.0)’
centre = ‘near(x[0], 0.0)’
bcu_inflow = DirichletBC(W.sub(0), (0.0, u_in), inflow)
bcu_wall1 = DirichletBC(W.sub(0), (0.0, u_c), wall1)
bcu_wall2 = DirichletBC(W.sub(0), (0.0, u_c), wall2)
#bcu_sym = Symmetry condition along centre, tangential stress = 0.
bcu_outflow = DirichletBC(W.sub(0), (0.0, u_c), outflow)
#bcs = [bcu_inflow, bcu_wall1, bcu_sym, bcu_outflow]
bcs = [bcu_inflow, bcu_wall1, bcu_wall2, bcu_outflow]

#Define the variational form
epsilon = sym(grad(u))
f = Constant((0, -9.81))
#F = (2muinner(epsilon, grad(v)) - div(u)*q - div(v)p)dx + p0dot(v,n)ds
F = (2muinner(epsilon, grad(v)) - div(u)*q - div(v)*p + dot(f,v))dx + p0dot(v,n)*ds

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

#Plot solutions
(u, p) = w.split()
File(“Stokes/velocity.pvd”) << u
File(“Stokes/pressure.pvd”) << p

Any help would be greatly appreciated!!

ExampleFullDomain
ExampleFullDomain1486×512 84 KB

2

I finally realised how to format my code correctly, I reattach it here

from dolfin import *

#Define mesh and geometry
mesh = RectangleMesh(Point(-1, 0), Point(1, 1), 120, 60)

#Mesh that I want, half of my current domain
#mesh = RectangleMesh(Point(0, 0), Point(1, 1), 120, 60)
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 = Expression(‘exp(-a*pow(x[0],2))’, degree=2, a=5)
p0 = Expression(“1”, degree=1)
u_in = Constant(-2.0)
u_c = Constant(-1.0)
inflow = ‘near(x[1], 1.0) && x[0]<0.5 && x[0]>-0.5’
wall1 = ‘near(x[0], 1.0)’
wall2 = ‘near(x[0], -1.0)’
outflow = ‘near(x[1], 0.0)’
centre = ‘near(x[0], 0.0)’
bcu_inflow = DirichletBC(W.sub(0), (0.0, u_in), inflow)
bcu_wall1 = DirichletBC(W.sub(0), (0.0, u_c), wall1)
bcu_wall2 = DirichletBC(W.sub(0), (0.0, u_c), wall2)
#bcu_sym = Symmetry condition along centre, tangential stress = 0.
bcu_outflow = DirichletBC(W.sub(0), (0.0, u_c), outflow)
#bcs = [bcu_inflow, bcu_wall1, bcu_sym, bcu_outflow]
bcs = [bcu_inflow, bcu_wall1, bcu_wall2, bcu_outflow]

#Define the variational form
epsilon = sym(grad(u))
f = Constant((0, -9.81))
#F = (2muinner(epsilon, grad(v)) - div(u)*q - div(v)p)dx + p0dot(v,n)ds
F = (2muinner(epsilon, grad(v)) - div(u)*q - div(v)*p + dot(f,v))dx + p0dot(v,n)*ds

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

#Plot solutions
(u, p) = w.split()
File(“Stokes/velocity.pvd”) << u
File(“Stokes/pressure.pvd”) << p

Sorry about that!

3

How did you get this term in F?

p0dot(v,n)ds
4

Hi,

you can simply impose a symmetry boundary condition for Stokes (or Navier-Stokes) flow as follows. The symmetry boundary condition states, that there must be no flux through the surface of symmetry and that the tangential component(s) are free. Therefore, the condition is given by

u \cdot n = 0

where n is the unit normal on the symmetry boundary.
As your symmetry boundary is given by the vertical line at x[0] = 0, the unit normal is given by n= [0,1]^T, so that the above condition is equivalent to u[0] = 0.

Therefore, in your code, you can impose this condition using

bcu_sym = DirichletBC(W.sub(0).sub(0), Constant(0.0), centre)

But the question that @Navier raised regarding your weak form, is of course still valid. Are you trying to impose a certain pressure level at the outlet?