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

kamilova-maths · March 8, 2021, 7:25pm

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.

nate · March 8, 2021, 8:15pm

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

kamilova-maths · March 9, 2021, 8:39am

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

nate · March 9, 2021, 3:32pm

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

kamilova-maths · March 9, 2021, 4:16pm

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

nate · March 12, 2021, 2:17am

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)

kamilova-maths · March 12, 2021, 8:34am

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.

Topic		Replies	Views
Slip condition on flat boundary, Stokes with constant viscosity variational formulation	2	404	February 19, 2021
Mini element - external force variational formulation	2	594	February 7, 2020
How to impose symmetry conditions on domain centreline, Stokes equations? variational formulation	3	563	November 30, 2022
Why is the normal stress non-zero at boundary where I impose zero normal stress? variational formulation	3	567	March 19, 2024
Steady State Stokes Equations with Iterative Solver mathematics	2	1084	March 24, 2020

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

Related topics