# How can I solve Stokes flow with an external force function correctly?

Hi all,
I am solving Stokes flow on a rectangular domain, with constant viscosity. I believe that my combination of boundary conditions is incorrect, as my external force is scaled out of the solution, i.e. I get the same thing no matter what f is. I have been reading about how to circumvent this problem by imposing alternative conditions at the bottom, but I can’t seem to find the right combination.

The top of my domain is split into two halves, ds0 (left) and ds1 (right). I want to impose an inlet velocity on ds0, and zero stress on ds1. Then an outlet velocity at the bottom (ds3). Is there a way to impose something equivalent to this that does not scale out my external force f? I attach the code below.

from dolfin import *
from mshr import *
import matplotlib.pyplot as plt
import numpy as np
# Define mesh and geometry
mesh = RectangleMesh(Point(0, 0), Point(0.2, 1), 60, 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 = Constant(1.0)
u_in = Constant(-2.0)
u_c = Constant(-1.0)

inflow = 'near(x[1], 1.0) && x[0]<=0.1'
weird = 'near(x[1], 1.0) && x[0]>=0.1'
wall = 'near(x[0], 0.2)'
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_wall, bcu_inflow, bcu_symmetry, bcu_outflow]

# Define stress tensor
x = SpatialCoordinate(mesh)

def epsilon(v):
return (as_tensor([[2*v[0].dx(0), v[0].dx(1) + v[1].dx(0), 0],
[v[1].dx(0) + v[0].dx(1), 2*v[1].dx(1), 0],
[0, 0, 0]]))

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

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

# Define the variational form
f = Constant((0, -1)) # Nothing happens when I change this
colors = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
colors.set_all(0)  # default to zero

CompiledSubDomain("near(x[0], 0.0)").mark(colors, 4)
CompiledSubDomain("near(x[1], 1.0) && x[0]<=0.1").mark(colors, 0)
CompiledSubDomain("near(x[1], 1.0) && x[0]>=0.1").mark(colors, 1)
CompiledSubDomain("near(x[0], 0.2)").mark(colors, 2)  # wall

CompiledSubDomain("near(x[1], 0.0)").mark(colors, 3)  # outflow

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

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

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

# Plot solutions
(u, p) = w.split()
#File("Results/velocitySym.pvd") << u


I would really appreciate it if anyone had any ideas. I am very confused.

Can you make sure your code example is runnable? You have a hanging return ... line which is indented without a function definition.

Hi, I’ve just edited it now. Apologies, I wanted to delete the comments to make it shorter and ended up accidentally deleting def epsilon(v):

1 Like

I’m pretty sure you’re just running into the problem you had last time regarding incompressible flow. The results you get look correct to me. Physical intuition would show that what comes in at the top left just hits your prescribed block of flow in the lower portion of the domain, and exits at the top right. Imposing a greater body force will do nothing because your fluid cannot compress itself within your boundary constraints.

If I change your boundary conditions to free slip, and make the outflow boundary a true natual zero stress Neumann condition, imposing various body forces has an impact since the fluid now has “somewhere to go”.

f = (0, 0)^\top, f = (0, -100)^\top , f = (0, 100)^\top

1 Like

Hello,
What you say makes perfect sense, I am just going in circles. I see that now.

However, I was also trying to do the slip condition as a comparison to the zero stress condition, but I was unable to get the results you’re showing… on which boundary are you imposing the slip condition?

[The following is a description of what I tried, but clearly I’m imposing free slip on a different boundary than you - feel free to ignore it ]

I have tried removing bcu_outflow and replacing with

right = 'near(x[1], 1.0) && x[0]>=0.1' # right half of top boundary
bcu_slip = DirichletBC(W.sub(0).sub(1), Constant(0.0), right) # slip condition


and I get the image

which makes sense to me as the slip condition is essentially putting a lid on the top right boundary. It’s also invariant to changes in f, as like you said, the flow has nowhere to go. How are you getting recirculation zones? It seems like you are applying different inlet velocities rather than imposing a slip condition on that top right.

This is my full slip condition code (pretty much the same as above, just replacing that boundary condition)

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(0.2, 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)
u_in = Constant(-2.0) # inlet velocity
u_c = Constant(-1.0) # wall velocity, and outlet velocity

inflow = 'near(x[1], 1.0) && x[0]<=0.1' # left half of top boundary
wall = 'near(x[0], 0.2)' # right boundary
centre = 'near(x[0], 0.0)' # left boundary
outflow = 'near(x[1], 0.0)' # bottom boundary
right = 'near(x[1], 1.0) && x[0]>=0.1' # right half of top boundary
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)
bcu_slip = DirichletBC(W.sub(0).sub(1), Constant(0.0), right) # slip condition
bcs = [bcu_inflow, bcu_wall, bcu_symmetry, bcu_slip]

x = SpatialCoordinate(mesh)

# Define stress tensor
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)) # forcing term

colors = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
colors.set_all(0)  # default to zero
CompiledSubDomain("near(x[0], 0.0)").mark(colors, 4)
CompiledSubDomain("near(x[1], 1.0) && x[0]<=0.1").mark(colors, 0) # top left
CompiledSubDomain("near(x[1], 1.0) && x[0]>=0.1").mark(colors, 1) # top right
CompiledSubDomain("near(x[0], 0.2)").mark(colors, 2)  # wall
CompiledSubDomain("near(x[1], 0.0)").mark(colors, 3)  # outflow

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

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

F = a1 + a2

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

# Plot solutions
(u, p) = w.split()

File("Results/velocityIsothermal.pvd") << u
c = plot(u, title='Velocity Isothermal')
plt.colorbar(c)
plt.show()


I believe this is the code I used to produce the plots. Don’t forget that I just meddled with the boundary conditions to create pretty pictures. You should make sure the code matches the underlying physical model you’re attempting to discretise. In your code it looks like you mixed up the dyadic (x, y) pair when you enforce \vec{u} \cdot \vec{n} = 0 on the right boundary.

from dolfin import *
import matplotlib.pyplot as plt
import numpy as np
# Define mesh and geometry
mesh = RectangleMesh(Point(0, 0), Point(0.2, 1), 60, 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 = Constant(1)
u_in = Constant(-2.0)
u_c = Constant(-1.0)

inflow = 'near(x[1], 1.0) && x[0]<=0.1'
wall = 'near(x[0], 0.2)'
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).sub(0), 0.0, wall)
bcu_symmetry = DirichletBC(W.sub(0).sub(0), Constant(0.0), centre)
bcs = [bcu_wall, bcu_inflow, bcu_symmetry]

# Define stress tensor
x = SpatialCoordinate(mesh)

def epsilon(v):
return (as_tensor([[2*v[0].dx(0), v[0].dx(1) + v[1].dx(0), 0],
[v[1].dx(0) + v[0].dx(1), 2*v[1].dx(1), 0],
[0, 0, 0]]))

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

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

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

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

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


Thank you so much for the code and the explanation.
My condition at the wall is that the velocity should be -1 downwards, rather than zero normal velocity, which creates a boundary layer near the wall as the velocity adjusts to -1. I understand the difference now and I understand why my solution is unaffected by f.

1 Like