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

**URL:** <https://fenicsproject.discourse.group/t/how-can-i-solve-stokes-flow-with-an-external-force-function-correctly/5266>\
**Category:** variational formulation\
**Created:** [March 8, 2021, 7:25pm UTC](https://fenicsproject.discourse.group/t/how-can-i-solve-stokes-flow-with-an-external-force-function-correctly/5266 "2021-03-08T19:25:07Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![kamilova-maths](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/kamilova-maths/32/2272_2.png) [@kamilova-maths](https://fenicsproject.discourse.group/u/kamilova-maths)\
**Post date:** [March 8, 2021, 7:25pm UTC](https://fenicsproject.discourse.group/t/how-can-i-solve-stokes-flow-with-an-external-force-function-correctly/5266/1 "2021-03-08T19:25:08Z")

</div>

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.

```auto
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.

---

<div class="post-metadata">

**Author:** ![nate](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nate/32/17_2.png) [@nate](https://fenicsproject.discourse.group/u/nate)\
**Post date:** [March 8, 2021, 8:15pm UTC](https://fenicsproject.discourse.group/t/how-can-i-solve-stokes-flow-with-an-external-force-function-correctly/5266/2 "2021-03-08T20:15:26Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![kamilova-maths](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/kamilova-maths/32/2272_2.png) [@kamilova-maths](https://fenicsproject.discourse.group/u/kamilova-maths)\
**Post date:** [March 9, 2021, 8:39am UTC](https://fenicsproject.discourse.group/t/how-can-i-solve-stokes-flow-with-an-external-force-function-correctly/5266/3 "2021-03-09T08:39:05Z")

</div>

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):

---

<div class="post-metadata">

**Author:** ![nate](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nate/32/17_2.png) [@nate](https://fenicsproject.discourse.group/u/nate)\
**Post date:** [March 9, 2021, 3:32pm UTC](https://fenicsproject.discourse.group/t/how-can-i-solve-stokes-flow-with-an-external-force-function-correctly/5266/4 "2021-03-09T15:32:27Z")

</div>

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

 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/b/b575a2361e69bb7dfe0b7027d0e5b6f43feb772a.png) ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/a/a026869437273f3373d5185f9603a1194aee2b23.png) ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/f/fd8d4ea0a1c23b2402bd9c1d543701e5b2df7cc0.png)

---

<div class="post-metadata">

**Author:** ![kamilova-maths](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/kamilova-maths/32/2272_2.png) [@kamilova-maths](https://fenicsproject.discourse.group/u/kamilova-maths)\
**Post date:** [March 9, 2021, 4:16pm UTC](https://fenicsproject.discourse.group/t/how-can-i-solve-stokes-flow-with-an-external-force-function-correctly/5266/5 "2021-03-09T16:16:27Z")

</div>

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

```auto
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

 ![VelocityFreeSlip](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/d/db677ce514c9ea579fdd3c34304d0728b3d8bd96.png)

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)

```auto
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()

```

---

<div class="post-metadata">

**Author:** ![nate](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nate/32/17_2.png) [@nate](https://fenicsproject.discourse.group/u/nate)\
**Post date:** [March 12, 2021, 2:17am UTC](https://fenicsproject.discourse.group/t/how-can-i-solve-stokes-flow-with-an-external-force-function-correctly/5266/6 "2021-03-12T02:17:03Z")

</div>

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.

```python
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)

```

---

<div class="post-metadata">

**Author:** ![kamilova-maths](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/kamilova-maths/32/2272_2.png) [@kamilova-maths](https://fenicsproject.discourse.group/u/kamilova-maths)\
**Post date:** [March 12, 2021, 8:34am UTC](https://fenicsproject.discourse.group/t/how-can-i-solve-stokes-flow-with-an-external-force-function-correctly/5266/7 "2021-03-12T08:34:14Z")

</div>

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.