Hello Everyone,
I am trying to conduct three point bending test. i want to apply contact boundary condition like the figure show. how to perform it in FEniCS？
I am trying to apply a point boundary condition. But the result is not right. I think that the boundary conditions are different from the experimental ones. Here is my attempt
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
#creat mesh
##---------L=4W------------------------
wt=10.01
lt=4.0*wt
wn=0.35
an=4.9
mesh = Mesh()
editor = MeshEditor()
editor.open(mesh,'triangle', 2, 2)
editor.init_vertices(16)
editor.init_cells(16)
editor.add_vertex(0, Point(0.0, 0.0))
editor.add_vertex(1, Point(wt, 0.0))
editor.add_vertex(2, Point(2.0*wt-wn/2.0, 0.0))
editor.add_vertex(3, Point(2.0*wt, an))
editor.add_vertex(4, Point(2.0*wt+wn/2.0, 0.0))
editor.add_vertex(5, Point(3.0*wt, 0.0))
editor.add_vertex(6, Point(4.0*wt, 0.0))
editor.add_vertex(7, Point(0.0, an))
editor.add_vertex(8, Point(wt, an))
editor.add_vertex(9, Point(3.0*wt, an))
editor.add_vertex(10, Point(4.0*wt, an))
editor.add_vertex(11, Point(0.0, wt))
editor.add_vertex(12, Point(wt, wt))
editor.add_vertex(13, Point(2.0*wt,wt))
editor.add_vertex(14, Point(3.0*wt, wt))
editor.add_vertex(15, Point(4.0*wt, wt))
editor.add_cell(0, [0, 1, 7])
editor.add_cell(1, [1, 8, 7])
editor.add_cell(2, [1, 2, 8])
editor.add_cell(3, [2, 3, 8])
editor.add_cell(4, [4, 5, 3])
editor.add_cell(5, [5, 9, 3])
editor.add_cell(6, [5, 6, 9])
editor.add_cell(7, [6, 10, 9])
editor.add_cell(8, [7, 8, 11])
editor.add_cell(9, [8, 12, 11])
editor.add_cell(10, [8, 3, 12])
editor.add_cell(11, [3, 13, 12])
editor.add_cell(12, [3, 9, 13])
editor.add_cell(13, [9, 14, 13])
editor.add_cell(14, [9, 10, 14])
editor.add_cell(15, [10, 15, 14])
editor.close()
mesh=refine(mesh)
##---------mesh end -------------------
#--------------------------------------
# change num_refine for mesh refinement
#--------------------------------------
def refine_elemnt(mesh):
makers = MeshFunction("bool", mesh, mesh.topology().dim())
makers.set_all(False)
for cell in cells(mesh):
if abs(cell.midpoint().x()-23.0)<3.0:
makers[cell.index()] = True
mesh=refine(mesh,makers)
return mesh
num_refine = 5
for i in range(num_refine):
mesh=refine(mesh)
for i in range(0):
mesh=refine_elemnt(mesh)
#Import mesh file
#mesh = Mesh('mesh.xml')
#mesh = RectangleMesh(Point(0.0, 0.0), Point(46.0, 10.01), 400, 80)
## Output: the infomation of mesh
print(mesh.hmax())
print(mesh.hmin())
print(mesh.num_cells())
print(mesh.num_vertices())
print(mesh)
file_mesh = File("meshwang.pvd")
file_mesh << mesh
# Define function space
W = VectorFunctionSpace(mesh, 'CG', 1)
u, v = TrialFunction(W), TestFunction(W)
# Introduce manually the material parameters
E0 =3.04e3 #N/mm/mm test:3.04e3
nu=0.36 # test:0.36
mu = E0/(2*(1.0+nu)) #N/mm/mm test: E0/(2*(1.0+nu)) hirshikesh: 8.0e3
lmbda = E0*nu/((1.0+nu)*(1.0-2.0*nu)) #N/mm/mm test:E0*nu/((1.0+nu)*(1.0-2.0*nu)) hirshikesh:12.0e3
# Constituive functions
def epsilon(u):
return sym(grad(u))
def sigma(u):
return 2.0*mu*epsilon(u)+lmbda*tr(epsilon(u))*Identity(len(u))
# Boundary conditions
#bot = CompiledSubDomain("near(x[1], 0.0) && on_boundary")
def top(x, on_boundary):
tol = DOLFIN_EPS
return (abs(x[1] - wt) < tol) and on_boundary
def pointl(x, on_boundary):
tol = DOLFIN_EPS
return (abs(x[0] - 0.0) < tol) and (abs(x[1] -0.0) < tol)
def pointr(x, on_boundary):
tol = DOLFIN_EPS
return (abs(x[0] - 4.0*wt) < tol) and (abs(x[1] -0.0) < tol)
def pointm(x, on_boundary):
tol = DOLFIN_EPS
return (abs(x[0] - 2.0*wt) <= tol) and (abs(x[1] -10.01) <= tol)
load = Expression("t", t = 0.0, degree=1)
bc_ptl = DirichletBC(W, Constant((0.0,0.0)), pointl,method="pointwise")
bc_ptr = DirichletBC(W.sub(1), Constant(0.0), pointr,method="pointwise")
bc_ptm = DirichletBC(W.sub(1), load, pointm,method="pointwise")
bc_u = [bc_ptl, bc_ptr,bc_ptm ]
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
AutoSubDomain(top).mark(boundaries,1)
ds = Measure("ds")(subdomain_data=boundaries)
n = FacetNormal(mesh)
#plot(mesh)
# Apply a point load equal to 100 N downward in Y direction at coordinate (23.0,10.01):
class Load(UserExpression):
def __init__(self, **kwargs):
super().__init__(degree=kwargs["degree"])
self.point = kwargs['point']
self.value = kwargs['value']
def eval(self, value, x):
if near(x[0], self.point[0]) and near(x[1], self.point[1]):
value[0] = self.value[0]
value[1] = self.value[1]
else:
value[0] = 0
value[1] = 0
def value_shape(self):
return (2,)
Tra = Load(point=(20.02, 10.02), value=(0.0, -100.0), degree=1)
# Variational form
a = inner(sigma(u), grad(v))*dx
f = Constant((0.0, 0.0))
L = dot(f, v)*ds(1)
t = 0
deltaT=0.01 # load up
load_dis =[] # load displacement
load_force =[] # load force
forcex = []
forcey = []
x_dofs = W.sub(0).dofmap().dofs()
y_dofs = W.sub(1).dofmap().dofs()
# Staggered scheme
while t<0.35: # real 0.007 test 0.38 mm
# Compute solution
u = Function(W)
solve(a == L, u, bc_u)
f_ext_unknown = assemble(inner(sigma(u), grad(v))*dx)
bc_ptm.apply(f_ext_unknown)
Fy = 0
for i in y_dofs:
Fy += f_ext_unknown[i]
forcey.append(-Fy)
load_dis.append(t)
t += deltaT
load.t=t
x1 = [0.0,0.03081,0.07841,0.11827,0.16532,0.20054,0.25952,0.27453,0.30665,0.34078,0.35061]
y1 = [0.0,8.53131,21.25996,32.32954,44.35104,53.43517,69.61796,73.61607,80.71474,89.06452,90.69639]
plt.figure(1)
plt.plot(load_dis, forcey,'b-',label='simulation load'.format(0))
plt.plot(x1, y1,'o',label='test load'.format(0))
plt.legend()
plt.show()
print ('Simulation completed')