Darcy Flow in 3D

1

Hi,

I’m trying to implement the darcy flow in 3D. The analogue in 2D works fine, however in 3D it yields weird results. I would assume an evenly distributed inflow in contrast to what can be seen in the image. Using different boundary conditions (velocity Dirichlet-BC for the inflow) works in 3D as well, but I need the BCs to be as they are. So my question is, whether it is possible to realize that.

Here is the minimum working example:

from dolfin import *

class Bottom(SubDomain):
    def inside(self,x, on_boundary):
        return near(x[2],0) and on_boundary

class Top(SubDomain):
    def inside(self,x, on_boundary):
        return near(x[2],1) and on_boundary

def wall(x, on_boundary):
        return (near(x[0],0) or near(x[0],1) or near(x[1],0) or near(x[1],1)) and on_boundary

### MESH & FUNCTION SPACE ###

mesh    = UnitCubeMesh(10,10,10)

Qe      = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
Be      = FiniteElement("Bubble",   mesh.ufl_cell(), 4)
Ve      = VectorElement(NodalEnrichedElement(Qe, Be))
element = MixedElement(Ve,Qe)
W       = FunctionSpace(mesh, element)

n       = FacetNormal(mesh)
u,p     = TrialFunctions(W)
v,q     = TestFunctions(W)

### SUB-DOMAINS ###

bottom = Bottom()
top    = Top()

sub_domains = MeshFunction('size_t', mesh, 1)
sub_domains.set_all(0)

bottom.mark(sub_domains, 1)
top.mark(sub_domains, 2)

ds = Measure('ds', domain=mesh, subdomain_data=sub_domains)

### BOUNDARY CONDITION ###

bc = DirichletBC(W.sub(0), Constant((0,0,0)), wall)

### VARIATIONAL FORMULATION ###

a =  (dot(u,v) - div(v)*p - div(u)*q)*dx
L = -Constant(1)*dot(n,v)*ds(1) + Constant(2)*dot(n,v)*ds(2)

w = Function(W)

solve(a == L, w, bc, solver_parameters={'linear_solver': 'mumps'})

file = File("flow.pvd")
file << w.split()[0]

file = File("pressure.pvd")
file << w.split()[1]

file = File("subdomains.pvd")
file << sub_domains

The resulting flow looks like this:

flow
flow1271×755 134 KB

2

In your line

you should use

sub_domains = MeshFunction('size_t', mesh, 2)

instead. This is because your boundaries of a 3D mesh are Faces, i.e., 2D objects, whereas for your 2D implementation the boundaries are 1D elements, i.e., lines.
More generally, you could use:

sub_domains = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)

as MeshFunction for boundaries.

3 Likes