# How can point loads be applied to the end of the cantilever beam?

``````from fenics import *
from mshr import *

L = 2
H = 0.067
b = 1

EI = 5e6
Ny = 100
Nx = 100

beam = RectangleMesh(Point(0, 0), Point(L, H), Nx, Ny)
beam.init()

def eps(v):

E = Constant(200e9)
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)

f = 0
w = 4e3
P = 6e3
f_thernal = Constant((0, -f))
T_plane = Constant((0, -w))
T_point = Constant((0, -P))
V = VectorFunctionSpace(beam, "CG", 2)

class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x, H) and on_boundary
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x, 0) and on_boundary
class Point(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and abs(x-L) < 0.01 and abs(x-H) < 0.01

boundaries = MeshFunction("size_t", beam, 1)

left = Left()
pt = Point()
top = Top()

boundaries.set_all(0)
top.mark(boundaries, 1)
left.mark(boundaries, 2)
pt.mark(boundaries, 3)
ds = Measure("ds", domain=beam, subdomain_data=boundaries)

du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_))*dx

left_end_displacement = Constant((0.0, 0.0))
bc = [DirichletBC(V, left_end_displacement, boundaries, 2)]

l = dot(T_point, u_)*ds(3)

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

plot(1e3*u, mode="displacement")

print("Maximal deflection:", -u(L, H/2))
print("Beam theory deflection:", float(((P*L**3)/(3*EI))))
``````

I’m trying to apply a point load to the end of the cantilever beam. I kept searching for Google and tried to remove the pointwise of bc and redefine the point as a small area, but the result was not good.
I’m looking for a way to fix this. The code above works fine without errors, but the result is strange. I don’t think there’s any grammatical error.

Hello

The below illustrates how you can apply a point load on a cantilever beam:

``````from dolfin import *

L = 10
W = 1
mesh = RectangleMesh(Point(0, 0), Point(L, W), 20, 3)
V = VectorFunctionSpace(mesh, 'P', 2)
tol = 1e-6
class LeftEdge(SubDomain):
def inside(self, x, on_boundary):

return on_boundary and abs(x) < tol

class RightEdge(SubDomain):
def inside(self, x, on_boundary):

return on_boundary and abs(x - L) < tol

def eps(v):

E = Constant(200e9)
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)

left_edge = LeftEdge()
right_edge = RightEdge()

boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
left_edge.mark(boundaries, 2)

#Zero Dirichlet B.C on the left side
bc = DirichletBC(V, Constant((0, 0)), boundaries , 2)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)

def __init__(self, **kwargs):
super().__init__()

self.point = kwargs['point']
self.value = kwargs['value']
def eval(self, value, x):
if near (x, self.point) and near(x, self.point):
value = self.value
value = self.value
else:
value = 0
value = 0
def value_shape(self):
return (2,)
T = PointLoad(point=(10.,1.), value=(0,-1e9), degree=1)
#Solve
a = inner(sigma(u), eps(v))*dx
L =dot(T, v)*ds

# Compute solution
u = Function(V)

solve(a == L, u, bc)
File("disp.pvd") << u
``````
``````from dolfin import *

L = 2
W = 0.067
mesh = RectangleMesh(Point(0, 0), Point(L, W), 20, 3)
V = VectorFunctionSpace(mesh, 'P', 2)
tol = 1e-6
class LeftEdge(SubDomain):
def inside(self, x, on_boundary):

return on_boundary and abs(x) < tol

class RightEdge(SubDomain):
def inside(self, x, on_boundary):

return on_boundary and abs(x - L) < tol

def eps(v):

E = Constant(200e9)
EI = 5e6
P = 6e3
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)

left_edge = LeftEdge()
right_edge = RightEdge()

boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
left_edge.mark(boundaries, 2)

#Zero Dirichlet B.C on the left side
bc = DirichletBC(V, Constant((0, 0)), boundaries , 2)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)

def __init__(self, **kwargs):
super().__init__()

self.point = kwargs['point']
self.value = kwargs['value']
def eval(self, value, x):
if near (x, self.point) and near(x, self.point):
value = self.value
value = self.value
else:
value = 0
value = 0
def value_shape(self):
return (2,)
T = PointLoad(point=(2, 0.067), value=(0,-P), degree=1)
#Solve
a = inner(sigma(u), eps(v))*dx
l =dot(T, v)*ds

# Compute solution
u = Function(V, name="Displacement")
solve(a == l, u, bc)

plot(1e3*u, mode="displacement")

print("Maximal deflection:", -u(L, W/2))
print("Beam theory deflection:", float(((P*L**3)/(3*EI))))
``````

I checked your answer and calculated it by modifying it according to the physical properties I wanted. However, the results do not match the theoretical values. I would appreciate it if you could check if I made any mistakes.

Solving linear variational problem.
Object cannot be plotted directly, projecting to piecewise linears.
Maximal deflection: 6.496007594039327e-05
Beam theory deflection: 0.0032

L = 2
W = 0.067
EI = 5e6
P(point load) = 6e3

In the above implementation, `T` is interpreted as a traction which is interpolated into `V` and integrated over the external surface of the beam. The figure below shows the y-component of the traction as it is interpreted by your code. This is very different from applying a point load.

To implement a point load, you should use `PointSource`, e.g. as below

``````from dolfin import *

L = 2
W = 0.067
mesh = RectangleMesh(Point(0, 0), Point(L, W), 20, 3)
V = VectorFunctionSpace(mesh, 'P', 2)
tol = 1e-6
class LeftEdge(SubDomain):
def inside(self, x, on_boundary):

return on_boundary and abs(x) < tol

class RightEdge(SubDomain):
def inside(self, x, on_boundary):

return on_boundary and abs(x - L) < tol

def eps(v):

E = Constant(200e9)
EI = 5e6
P = 6e3
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)

left_edge = LeftEdge()
right_edge = RightEdge()

boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
left_edge.mark(boundaries, 2)

#Zero Dirichlet B.C on the left side
bc = DirichletBC(V, Constant((0, 0)), boundaries , 2)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)

T = PointSource(V.sub(1), Point(2, 0.067), -P)

#Solve
a = inner(sigma(u), eps(v))*dx
l = dot(Constant((0.0,0.0)), v)*ds

K = assemble(a)
b = assemble(l)
T.apply(b)
bc.apply(K)
bc.apply(b)

# Compute solution
u = Function(V, name="Displacement")
solve(K, u.vector(), b)

plot(1e3*u, mode="displacement")

print("Maximal deflection:", -u(L, W/2))
print("Beam theory deflection:", float(((P*L**3)/(3*EI))))
``````

producing

``````Maximal deflection: 0.0031860791474304178
Beam theory deflection: 0.0032
``````
3 Likes

I really appreciate. But I have another problem. My goal is to calculate the deflection amount of the cantilever beam with a uniform distribution load and a point load at the same time. However, it is difficult to apply the uniform distribution load in the code you gave me. Is there a way to apply point load and uniform distribution load at the same time?

I solved the interpretation of the uniform distribution load using the following code.

class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x, W) and on_boundary

w = 4e3

top = Top()

top.mark(boundaries, 1)

l = dot(Constant((0, -w)), v)*ds(1)

However, this code does not apply to the code you gave me with different solve format. I guess the solve format I used is solve (a==l, u, bc), but the solve format in the code you gave me is a matrix. I want to find a way to solve this problem.

To apply the distributed load, replace the following line in the code I provided:

``````l = dot(Constant((0.0,0.0)), v)*ds
``````

with your variational formulation, i.e.:

``````l = dot(Constant((0, -w)), v)*ds(1)
``````

Also, don’t forget to define a `Measure` for the surface integral after the boundaries have been marked, i.e.

``````ds = Measure('ds', mesh, subdomain_data=boundaries)
``````

A complete code is below:

``````from dolfin import *

L = 2
W = 0.067
mesh = RectangleMesh(Point(0, 0), Point(L, W), 20, 3)
V = VectorFunctionSpace(mesh, 'P', 2)
tol = 1e-6
class LeftEdge(SubDomain):
def inside(self, x, on_boundary):

return on_boundary and abs(x) < tol

class RightEdge(SubDomain):
def inside(self, x, on_boundary):

return on_boundary and abs(x - L) < tol

def eps(v):

E = Constant(200e9)
EI = 5e6
P = 6e3
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)

left_edge = LeftEdge()
right_edge = RightEdge()

boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
left_edge.mark(boundaries, 2)

class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x, W) and on_boundary
w = 4e3
top = Top()
top.mark(boundaries, 1)
ds = Measure('ds', mesh, subdomain_data=boundaries)

#Zero Dirichlet B.C on the left side
bc = DirichletBC(V, Constant((0, 0)), boundaries , 2)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)

T = PointSource(V.sub(1), Point(2, 0.067), -P)

# Variational formulations
a = inner(sigma(u), eps(v))*dx

# Solution #1: point load only
l = dot(Constant((0,0)), v)*ds
K = assemble(a)
b = assemble(l)
T.apply(b)
bc.apply(K)
bc.apply(b)
u_pt = Function(V, name="Displacement, point load")
solve(K, u_pt.vector(), b)
print("Maximal deflection, point load only:", -u_pt(L, W/2))

# Solution #2: distributed load only
l = dot(Constant((0, -w)), v)*ds(1)
K = assemble(a)
b = assemble(l)
bc.apply(K)
bc.apply(b)
u_dist = Function(V, name="Displacement, distributed load")
solve(K, u_dist.vector(), b)
print("Maximal deflection, distributed load:", -u_dist(L, W/2))

# Solution #3: distributed load plus point load
l = dot(Constant((0, -w)), v)*ds(1)
K = assemble(a)
b = assemble(l)
T.apply(b)
bc.apply(K)
bc.apply(b)
u = Function(V, name="Displacement, combined load")
solve(K, u.vector(), b)
print("Maximal deflection, combined load:", -u(L, W/2))

plot(1e3*u, mode="displacement")
print("Beam theory deflection:", float(((P*L**3)/(3*EI))))
``````
1 Like