Defining boundary conditions inside holes

General
1

Hey everyone.

I’m trying to use a self-made program to calculate displacements in a flat plate. This flat plate has a rectangular hole inside it. I want to define a boundary condition consisting on a constant force being applied only on 1 side of the rectangle (leaving the other 3 untouched). I’ve tried different ways but I don’t get it. The program works but the displacements are nearly 0 even with big forces. I feel like the problem is on the BC. Is there a way to do this properly?

Thanks

2

Without supplying any code along with this statement, it is really hard to tell what has gone wrong with your code.

I would suggest making a minimal reproducible example

3

Hey again.

Below I’m adding my actual code. I used Paraview to visualise the domain and everything looks okay. As I mentioned, the program runs properly but clearly states boundary conditions may not have been applied properly. As I also mentioned above, I raised my forces to 1e-1 and still displacements are almost 0. There are three holes, but I want to calculate the solution as if the force was only applied on the upper lip of the upper hole. There might be other programming issues as I’m new to this, but I may have not seen them. Thanks in advance.

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(-8.0,0.0), Point(8.0,0.0)])
domain = domain - Rectangle(Point(-0.5,5.),Point(0.5,6.79)) - Rectangle(Point(-0.5,11.79),Point(0.5,13.58)) - Rectangle(Point(-0.5,18.58),Point(0.5,20.37))
mesh = generate_mesh(domain,4)

vtkfile = File('holein/malla.pvd')
vtkfile <<  mesh

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

#caracterizamos material E en N/cm^2
E = Constant(3.3e5)
nu = Constant(0.45)
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,1e-3))

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

#contorno

class Left(SubDomain):
	def inside(self, x, on_boundary):
		return near (x[0], -8.) 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 near (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)
ds = Measure('ds', subdomain_data=facets)

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

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

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

Thanks