Dear,
I am trying to carry out a project of a cantilever in the multicomponent or multibody style. There would be two domains or two meshes connected to the center by circular elements, such as screws or welds. For the time being I am not considering the contact between the meshes or domains, only that the regions of the holes transmit the loads. See the image below.
I have a code below, but I think it is not as I would like. I believe that the main errors would be to consider how to actually insert two meshes or two overlapping domains and how to make the boundary conditions of the holes between one domain and another.
Can anybody help me? @bleyerj
from fenics import *
from mshr import *
import matplotlib.pyplot as plt
L1 = 4.0
H = 1.0
tol = 1E-14
# Load
F = 200
ds_size = 0.1
class Left(SubDomain):
def inside(self, x, on_boundary):
return (x[0] <= 2.5 + tol) #and on_boundary
class Right(SubDomain):
def inside(self, x, on_boundary):
return (x[0] >= 1.5 - tol) #and on_boundary
#Define geometry for background
x1 = 2.0
y1 = 0.25
x2 = 2.0
y2 = 0.75
r = 0.1
domain = Rectangle(Point(0, 0), Point(L1, H)) \
- Circle(Point(x1, y1), r) - Circle(Point(x2, y2), r)
mesh = Mesh(generate_mesh(domain, 200))
# Initialize mesh function for interior domains
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
domains.set_all(0)
Left().mark(domains, 1)
Right().mark(domains, 2)
# Define new measures associated with the interior domains
dx = Measure("dx", domain=mesh, subdomain_data=domains)
# Boundary Condition
class Left1(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0) and on_boundary
class Circle1(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near((x[0]-x1)**2+(x[1]-y1)**2, r**2, 1e-3)
class Circle2(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near((x[0]-x2)**2+(x[1]-y2)**2, r**2, 1e-3)
# Load
class LoadSurface(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], L1) and between(x[1], ((H-ds_size)/2.0, (H+ds_size)/2.0))
loadSurface = LoadSurface()
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
loadSurface.mark(boundaries, 1)
ds = ds(subdomain_data=boundaries)
t = Expression(('q_x', 'q_y'), degree = 2, q_x = 0.0, q_y = -F/ds_size)
## Constitutive relation
def eps(v):
return sym(grad(v))
E = Constant(1e5)
nu = Constant(0.3)
model = "plane_stress"
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
if model == "plane_stress":
lmbda = 2*mu*lmbda/(lmbda+2*mu)
def sigma(v):
return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)
## Variational formulation
## Resolution
V = VectorFunctionSpace(mesh, 'Lagrange', degree=2)
du = TrialFunction(V)
u_ = TestFunction(V)
left1 = DirichletBC(V, Constant((0.0, 0.0)), Left1())
uc = Function(V)
bc1 = DirichletBC(V.sub(0), Constant(0.0), Circle1())
bc1.apply(uc.vector())
bc2 = DirichletBC(V.sub(0), Constant(0.0), Circle2())
bc2.apply(uc.vector())
bcs = [left1, bc1, bc2]
a = inner(sigma(du), eps(u_))*dx(1) + inner(sigma(du), eps(u_))*dx(2)
L = dot(t, u_)*ds(1)
u = Function(V, name="Displacement")
solve(a == L, u, bcs, solver_parameters={'linear_solver': 'bicgstab',
'preconditioner': 'hypre_amg'})
plot(u, mode="displacement")
plt.title("displacement")
plt.show()
File("u.pvd") << u
print("Maximal deflection:", - u(L1, H/2.)[1])
print("Beam theory deflection:", float((4*F/ds_size*L1**3)/(E*H**3)))
Vsig = TensorFunctionSpace(mesh, "DG", degree=0)
sig = Function(Vsig, name="Stress")
sig.assign(project(sigma(u), Vsig))
print("Stress at (0,H):", sig(0, H))