Hi all. I’m an undergrad student modeling Antarctic ice flow around ice rumples as 2D Stokes flow around a circular obstacle. I read a few of the previous posts and tutorials on this, but unfortunately I’m still having trouble with defining my no-slip boundary conditions around the obstacle.
My mesh looks like this:
I want to establish no-slip boundary conditions around the obstacles. I did this using the residual method found here. The resulting velocity field looks like this:
However, I want to make my material non-Newtonian, so I decided to go back and use the normal linear variational solver that is used in most of the FEniCS tutorials.
First, I create a mesh, define subdomains, and export the mesh and subdomains as xml files:
mymesh = generate_mesh(Rectangle(Point(0.0,0.0), Point(L,H)) - Circle(Point(60,0.25*H),25) - Circle(Point(130, 0.75*H), 25), resolution)
File("meshes/mesh.xml.gz") << mymesh
mesh = Mesh("meshes/mesh.xml.gz")
# Sub domain for no-slip (mark whole boundary, inflow and outflow will overwrite)
class Noslip(SubDomain):
def inside(self, x, on_boundary):
return on_boundary
# Sub domain for inflow (right)
class Inflow(SubDomain):
def inside(self, x, on_boundary):
return x[0] > 1.0 - DOLFIN_EPS and on_boundary
# Sub domain for outflow (left)
class Outflow(SubDomain):
def inside(self, x, on_boundary):
return x[0] < DOLFIN_EPS and on_boundary
# Create mesh functions over the cell facets
sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
# Mark all facets as sub domain 3
sub_domains.set_all(3)
# Mark no-slip facets as sub domain 0, 0.0
noslip = Noslip()
noslip.mark(sub_domains, 0)
# Mark inflow as sub domain 1, 01
inflow = Inflow()
inflow.mark(sub_domains, 1)
# Mark outflow as sub domain 2, 0.2, True
outflow = Outflow()
outflow.mark(sub_domains, 2)
# Save sub domains to file
file = File("meshes/subdomains.xml.gz")
file << sub_domains
Then, I import the xml mesh files, define my variational problem with Taylor-Hood elements, define my Dirichlet boundary conditions, and then solve the problem.
mesh = Mesh("meshes/mesh.xml.gz")
sub_domains = MeshFunction("size_t", mesh, "meshes/subdomains.xml.gz")
# Define Taylor-Hood function space
VE = VectorElement("CG", mesh.ufl_cell(), 2)
QE = FiniteElement("CG", mesh.ufl_cell(), 1)
WE = VE * QE
W = FunctionSpace(mesh, WE)
V = FunctionSpace(mesh, VE)
Q = FunctionSpace(mesh, QE)
noslip = Constant((0, 0))
bc0 = DirichletBC(W.sub(0), noslip, sub_domains, 0) # subdomain marked 0 (noslip)
# Inflow boundary condition for velocity
# x0 = 1
inflow = Expression(("4*(x[1]*(YMAX-x[1]))/(YMAX*YMAX)", "0."), YMAX=YMAX, element = V.ufl_element())
bc1 = DirichletBC(W.sub(0), inflow, sub_domains, 1) # subdomain marked 1 (inflow)
# Boundary condition for pressure at outflow
# x0 = 0
zero = Constant(0)
bc2 = DirichletBC(W.sub(1), zero, sub_domains, 2) # subdomain marked 2 (outflow)
# Collect boundary conditions
bcs = [bc0, bc1, bc2]
eta = 1
I = Identity(2)
def epsilon(u): # strain rate tensor
return 0.5*(grad(u) + grad(u).T)
def sigma(u,p): # Cauchy stress tensor
return 2*eta*epsilon(u) - p*I # from https://fenicsproject.org/pub/tutorial/html/._ftut1008.html
# Define trial and test functions
u, p = TrialFunctions(W)
v, q = TestFunctions(W)
f = Constant((0, 0)) # set f, body force, to 0
a = inner(sigma(u,p), grad(v))*dx + q*div(u)*dx
L = inner(f, v)*dx
# Compute solution
w = Function(W)
solve(a == L, w, bcs)
The resulting velocity field looks like this:
As you can see, the velocity at the no-slip boundary on the obstacles is not 0, as it is in the previous velocity plot. My code runs without error, so I think it is able to find the appropriate facets to apply the boundary conditions to. Am I defining my no-slip boundary incorrectly? Any suggestions would be greatly appreciated!
