How to Define/Fix Boundary Conditions at Internal Corners

variational formulation
1

I am solving a linear elasticity problem on a quarter circle, and there are mixed BC applied to the mesh. Dirichlet on the left and bottom and Neumann on the top and right. However, when plotting the corners at where there is a circle cut out of a square the stress is incorrect, I am assuming this is due to units there not being defined. I have attached the necessary code and images for reference, I would appreciate any help please.

As seen in this image at (3,0) and (0,3) the stress contour is acting unusually, to know if the problem is correct the top of the quarter circle should be a value of 3 and as there are more mesh refinements the value at the point gets larger and larger.

cornerstress.PNG
cornerstress.PNG1301×1165 77 KB 

# Bottom Dirichlet BC
class bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0) and on_boundary

# Left Dirichlet BC   
class left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0) and on_boundary

# Right NBC   
class right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], l) and on_boundary # where l is length of square

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

V = VectorFunctionSpace(mesh, "CG", 1)
G = BoundarySource(mesh, degree=2)

# Define Boundary Conditions (fix)
bc = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
bc.set_all(0)
right().mark(bc, 1)
top().mark(bc, 2)
ds = fe.ds(subdomain_data=bc)
bcs = [DirichletBC(V, G, bottom()), DirichletBC(V, G, left())]

# Define Variational Problem
u = TrialFunction(V)
d = u.geometric_dimension()
v = TestFunction(V)
f = Constant((0, 0))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(G, v)*ds(1) + dot(G, v)*ds(2)
2

I would think that you would want to replace the line

bcs = [DirichletBC(V, G, bottom()), DirichletBC(V, G, left())]

with something like

bcs = [DirichletBC(V.sub(1),Constant(0),bottom()), DirichletBC(V.sub(0),Constant(0),left())]

if the bottom and left are intended to be symmetry planes. What this is doing is constraining the x-component of displacement (V.sub(0)) on the left and the y-component of displacement (V.sub(1)) on the bottom to be zero.

2 Likes