Boundary conditions not being applied

General
1

Hey there.

I’m playing with different geometries and loads being applied on specific parts of a domain. For this code, I’m using a triangle with holes, with a pressure P being applied on one lip of one hole (see pictures). As I observed a symmetry condition, I’m only working on one side of the triangle.

duda
duda907×790 5.74 KB

When it comes to boundary conditions and applying the pressure, I don’t get any ‘‘warning, error’’ message, but displacements and stresses are null, which is obviously not correct. The upper lip of the hole is on coordinates x=[0, 0.5] y=20.37.

Captura de pantalla 2023-04-11 a las 8.16.09
Captura de pantalla 2023-04-11 a las 8.16.092254×1268 72.7 KB

Here I add my full code, to see if anyone could help. I believe it has to do with the way I’m defining the boundary conditions or the subdomains, but I can’t help myself.

Thanks in advance.

from fenics import*
from dolfin import*
from mshr import*
from mshr import Polygon, generate_mesh

mesh=Mesh()

#longitudes en cm
domain = Polygon([Point(0.0,25.37), Point(0.0,0.0), Point(8.0,0.0)])
domain = domain - Rectangle(Point(0.,5.),Point(0.5,6.79)) - Rectangle(Point(0.,11.79),Point(0.5,13.58)) - Rectangle(Point(0.,18.58),Point(0.5,20.37))
mesh = generate_mesh(domain,4)
	
vtkfile = File('trabajo/domain.pvd')
vtkfile <<  mesh

def eps(v):
	return sym(grad(v))

#caracterizamos material E en N/cm^2
E = Constant(250)
nu = Constant(0.2)
model = "plane_stress"

G = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
if model == "plane_stress":
	lmbda = 2*G*lmbda/(lmbda+2*G)

def sigma(v):
	return lmbda*tr(eps(v))*Identity(2) + 2.0*G*eps(v)

#definir carga, cómo?
f=Constant((0,2e3))

V = VectorFunctionSpace (mesh, 'Lagrange', degree=2)
du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du),eps(u_))*dx
l = inner(f,u_)*ds(3)

#contorno

class Left(SubDomain):
	def inside(self, x, on_boundary):
		return near (x[0], 0) and on_boundary
class Right(SubDomain):
	def inside(self, x, on_boundary):
		return near (x[0], 8.) and on_boundary
class Bottom(SubDomain):
	def inside(self, x, on_boundary):
		return x[1]==0 and on_boundary
class Top(SubDomain):
	def inside(self, x, on_boundary):
		return near (x[0], 0.) and near (x[1], 20.37) and on_boundary

facets = MeshFunction("size_t", mesh, 1)
facets.set_all(0)
Left().mark(facets,1)
Right().mark(facets,2)
Bottom().mark(facets,3)
Top().mark(facets, 4)
ds = Measure('ds', subdomain_data=facets)

bc = [DirichletBC(V.sub(0),  Constant(0.), facets, 1)]

u=Function(V, name="Desplazamientos")
solve (a==l, u, bc)

vtkfile = File('trabajo/desplazamientos.pvd')
vtkfile << u

T = TensorFunctionSpace(mesh, 'CG', 1)
stress = Function(T,name="Stress")
stress.assign(project(sigma(u),T))

vtkfile = File('trabajo/tensiones.pvd')
vtkfile << stress
2

I guess the problem here is that you set near(x[0], 0.) here, which means you are trying to constrain the boundary condition to a single point (not the upper lip). If you want to mark all facets with that y coordinate, you need to use:

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

Hello

This was my first thought a couple of days ago. I tried again doing exactly the changes that you mentioned before, as I did days ago, and this is what I get:

Solving linear variational problem.
  *** Warning: Found no facets matching domain for boundary condition.
  *** Warning: Found no facets matching domain for boundary condition.

Still, displacements and stress magnitude are extremely close to 0.

Thanks again.

4

I cannot reproduce this with the following code:

from dolfin import*
from mshr import*

mesh=Mesh()

#longitudes en cm
domain = Polygon([Point(0.0,25.37), Point(0.0,0.0), Point(8.0,0.0)])
domain = domain - Rectangle(Point(0.,5.),Point(0.5,6.79)) - Rectangle(Point(0.,11.79),Point(0.5,13.58)) - Rectangle(Point(0.,18.58),Point(0.5,20.37))
mesh = generate_mesh(domain,4)
	
vtkfile = File('trabajo/domain.pvd')
vtkfile <<  mesh

def eps(v):
	return sym(grad(v))

#caracterizamos material E en N/cm^2
E = Constant(250)
nu = Constant(0.2)
model = "plane_stress"

G = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
if model == "plane_stress":
	lmbda = 2*G*lmbda/(lmbda+2*G)

def sigma(v):
	return lmbda*tr(eps(v))*Identity(2) + 2.0*G*eps(v)

#definir carga, cómo?
f=Constant((0,2e3))

V = VectorFunctionSpace (mesh, 'Lagrange', degree=2)
du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du),eps(u_))*dx
l = inner(f,u_)*ds(3)

#contorno

class Left(SubDomain):
	def inside(self, x, on_boundary):
		return near (x[0], 0) and on_boundary
class Right(SubDomain):
	def inside(self, x, on_boundary):
		return near (x[0], 8.) and on_boundary
class Bottom(SubDomain):
	def inside(self, x, on_boundary):
		return x[1]==0 and on_boundary
class Top(SubDomain):
	def inside(self, x, on_boundary):
		return near (x[1], 20.37) and on_boundary

facets = MeshFunction("size_t", mesh, 1)
facets.set_all(0)
Left().mark(facets,1)
Right().mark(facets,2)
Bottom().mark(facets,3)
Top().mark(facets, 4)
File("facets.pvd") << facets

yielding

image
image351×571 28.7 KB

5

I run your last code adding the next lines to obtain displacements and stress:

bc = [DirichletBC(V.sub(0),  Constant(0.), facets, 1)]

u=Function(V, name="Desplazamientos")
solve (a==l, u, bc)

vtkfile = File('trabajo/desplazamientos.pvd')
vtkfile << u

T = TensorFunctionSpace(mesh, 'CG', 1)
stress = Function(T,name="Stress")
stress.assign(project(sigma(u),T))

vtkfile = File('trabajo/tensiones.pvd')
vtkfile << stress

Displacements and stress are extremely close to 0. The force f seems irrelevant, as if it was not applied.

sdeda
sdeda1323×779 19.5 KB

This is my main problem. The force not being applied.