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.