Hi all,
TL;DR:
- Why isn’t the pressure boundary condition, with values from a previous solution, being applied on the border?
- For
near(x[0], 2.0), why isx[0] = 1.95being called? - How can I fix my code, i.e. apply the pressure boundary condition on the border?
I am trying to simulate fluid flow (Navier-Stokes without time dependency) in a long pipe by splitting it into several smaller simulations, as seen in the image below. I am using Dirichlet Boundary Conditions for velocity on the left and pressure on the right of each pipe.
In the first simulation, I assume a velocity and pressure boundary. The second simulation is connected to the first spatially, so I want to use the velocity profile on the right side of the first simulation as the profile on the left side of the second simulation, again setting pressure to a constant. This is working more-or-less as intended.
Figure 1: Splitting a long pipe into two sections, using adjacent edges from previous simulations as Dirichlet boundary conditions on future ones
In the third simulation, I want to use the pressure from the left side of the second simulation as the pressure profile on the right side of the third simulation. However, this part is not running correctly. It seems the pressure is still being set to 0, and the values being called to apply the boundary condition are not on the boundary.
My question has 3 parts.
- Why isn’t the pressure boundary condition being applied on the border?
- In simulation 3, why is
x[0] = 1.95being called? - How can I fix this to run as intended, i.e. apply the pressure boundary condition on the border?
Below is my best attempt at a MWE. I am running FEniCS 2018.2.0.dev0 on Ubuntu 18.04.1.
from fenics import *
class borderVelBC(UserExpression):
def getPrevSoln(self, prevSol):
self.prevSol = prevSol
def eval(self, values, x):
try:
val = self.prevSol(x[0], x[1])
except:
val = [0, 0]
values[0] = val[0]
values[1] = val[1]
def value_shape(self):
return (2,)
class borderPreBC(UserExpression):
def getPrevSoln(self, prevSol):
self.prevSol = prevSol
def eval(self, values, x):
try:
val = self.prevSol(x[0], x[1])
print(str(val) + " at " + str(x))
except:
print("attempt invalid at " + str(x))
val = -99
values = val
########## Simulation 1 ##########
lefBound = 0.0
rigBound = 2.0
botBound = 0.0
topBound = 1.0
mesh = RectangleMesh(Point(lefBound, botBound), Point(rigBound, topBound), 40, 20)
V = VectorElement("CG", mesh.ufl_cell(), 2)
Q = FiniteElement("CG", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, MixedElement((V, Q)))
noSlip = DirichletBC(W.sub(0), Constant((0.0, 0.0)), "near(x[1], 0.0) || near(x[1], 1.0)")
leftBC = DirichletBC(W.sub(0), Constant((0.01, 0.0)), "near(x[0], 0.0) && x[1] < 0.9 && x[1] > 0.1")
rightBC = DirichletBC(W.sub(1), Constant(0.0), "near(x[0], 2.0)")
bcs = [noSlip, leftBC, rightBC]
w = Function(W)
(u, p) = split(w)
(v, q) = TestFunctions(W)
F = inner(u, v)*dx + inner(grad(u), grad(v))*dx + inner(grad(p), v)*dx + inner(div(u), q)*dx
solve(F == 0, w, bcs = bcs)
(u_soln, p_soln) = w.split(deepcopy=True)
File("u1.pvd") << u_soln
File("p1.pvd") << p_soln
prevUSoln = u_soln
prevPSoln = p_soln
########## Simulation 2 ##########
lefBound = 2.0
rigBound = 4.0
botBound = 0.0
topBound = 1.0
mesh = RectangleMesh(Point(lefBound, botBound), Point(rigBound, topBound), 40, 20)
V = VectorElement("CG", mesh.ufl_cell(), 2)
Q = FiniteElement("CG", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, MixedElement((V, Q)))
noSlip = DirichletBC(W.sub(0), Constant((0.0, 0.0)), "near(x[1], 0.0) || near(x[1], 1.0)")
rightBC = DirichletBC(W.sub(1), Constant(0.0), "near(x[0], 4.0)")
bcFromSim = borderVelBC()
bcFromSim.getPrevSoln(prevUSoln)
leftBC = DirichletBC(W.sub(0), bcFromSim, "near(x[0], 2.0)")
bcs = [noSlip, leftBC, rightBC]
w = Function(W)
(u, p) = split(w)
(v, q) = TestFunctions(W)
F = inner(u, v)*dx + inner(grad(u), grad(v))*dx + inner(grad(p), v)*dx + inner(div(u), q)*dx
solve(F == 0, w, bcs = bcs)
(u_soln, p_soln) = w.split(deepcopy=True)
File("u2.pvd") << u_soln
File("p2.pvd") << p_soln
prevUSoln = u_soln
prevPSoln = p_soln
########## Sim. 1 Correction ##########
lefBound = 0.0
rigBound = 2.0
botBound = 0.0
topBound = 1.0
mesh = RectangleMesh(Point(lefBound, botBound), Point(rigBound, topBound), 40, 20)
V = VectorElement("CG", mesh.ufl_cell(), 2)
Q = FiniteElement("CG", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, MixedElement((V, Q)))
noSlip = DirichletBC(W.sub(0), Constant((0.0, 0.0)), "near(x[1], 0.0) || near(x[1], 1.0)")
leftBC = DirichletBC(W.sub(0), Constant((0.01, 0.0)), "near(x[0], 0.0) && x[1] < 0.9 && x[1] > 0.1")
bcFromSim = borderPreBC()
bcFromSim.getPrevSoln(prevPSoln)
rightBC = DirichletBC(W.sub(1), bcFromSim, "near(x[0], 2.0)")
bcs = [noSlip, leftBC, rightBC]
w = Function(W)
(u, p) = split(w)
(v, q) = TestFunctions(W)
F = inner(u, v)*dx + inner(grad(u), grad(v))*dx + inner(grad(p), v)*dx + inner(div(u), q)*dx
solve(F == 0, w, bcs = bcs)
(u_soln, p_soln) = w.split(deepcopy=True)
File("u1correction.pvd") << u_soln
File("p1correction.pvd") << p_soln
