BC's for BDM1 x DG0

1

Hi there,

I’ve implemented a darcy flow problem, where I apply the pressure BC weakly, and the velocity BC strongly which results in a free slip flow. It works fine for geometries with perpendicular boundaries as shown in fig. 1 but results in unexpected velocity spikes for geometries where the pressure BC surface and it’s adjacent velocity BC surface are not perpendicular to each other as shown in fig. 2.

Is there a way to get rid of those velocity spikes?

Strong velocity BC is marked red. Weak pressure BC is marked green.

Fig 1:

u2
u21417×391 154 KB

Fig2:

u
u1435×404 204 KB

Here is the code:

from dolfin import *
from matplotlib import pyplot as plt
from mshr import *

class Boundary(SubDomain):
    def inside(self,x, on_boundary):
        return on_boundary

class Inflow(SubDomain):
    def inside(self,x, on_boundary):
        return (.3 <= x[0] <= .7) and x[1] < .001 and on_boundary

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

domain_points = [Point(0, 0), 
                 Point(1, 0),
                 Point(1, 1),
                 Point(0, 1),]

domain  = Polygon(domain_points)
mesh    = generate_mesh(domain, 50)

BDM1    = FiniteElement('BDM',triangle, 1)
DG0     = FiniteElement('DG',triangle, 0)
E       = MixedElement(BDM1, DG0)
W       = FunctionSpace(mesh,E)

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

sub_domains = MeshFunction('size_t', mesh, 1)
sub_domains.set_all(0)
Boundary().mark(sub_domains, 3)
Top().mark(sub_domains, 2)
Inflow().mark(sub_domains, 1)

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

a   = (dot(u,v) - div(u)*q - div(v)*p)*dx
L   = Constant(0)*dot(n,v)*ds(2) - Constant(1)*dot(n,v)*ds(1)
bc  = DirichletBC(W.sub(0), Constant((0,0)),sub_domains,3)

w   = Function(W)

solve(a == L, w,bc)

u,p = w.split()

plt.figure(dpi=400)
c = plot(u)
plt.colorbar(c)
plt.savefig("u.png")
2

Hello,

I think what you are seeing here is a singularity that can occur in solutions to the Poisson equation when the boundary condition is changed from a Dirichlet type to a Neumann type along a straight boundary.

Thanks,

Joe

1 Like