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

mesh
1 
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):
	return sym(grad(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[1], H) and on_boundary
class Left(SubDomain):
	def inside(self, x, on_boundary):
		return near(x[0], 0) and on_boundary
class Point(SubDomain):
	def inside(self, x, on_boundary):
		return on_boundary and abs(x[0]-L) < 0.01 and abs(x[1]-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)[1])
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.

2

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[0]) < tol

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

       return on_boundary and abs(x[0] - L) < tol

def eps(v):
    return sym(grad(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)

class PointLoad(UserExpression):

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

        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,)
#End Load
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
3 
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[0]) < tol

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

       return on_boundary and abs(x[0] - L) < tol

def eps(v):
    return sym(grad(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)

class PointLoad(UserExpression):

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

        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,)
#End Load
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)[1])
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

4

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.

image
image1080×336 24.6 KB

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[0]) < tol

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

       return on_boundary and abs(x[0] - L) < tol

def eps(v):
    return sym(grad(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)

#End Load
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)[1])
print("Beam theory deflection:", float(((P*L**3)/(3*EI))))

producing

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

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?

6

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[1], 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.

7

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[0]) < tol

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

       return on_boundary and abs(x[0] - L) < tol

def eps(v):
    return sym(grad(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)

# Top load
class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 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)

#End Load
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)[1])

# 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)[1])

# 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)[1])

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

Good morning
Your code has been a great help already.
I am struggling to adapt it for a 1D beam
Could you please suggest a version ?

9

Please make a minimal attempt to add a point load to a 1D beam. Asking someone to write a new code from scratch is a big ask.

10

What do you mean 1D beam ? If you are referring to 1D beam elements based on Euler-Bernoulli or Timoshenko beam theory this is very different from a 2D/3D solid formulation as above. In this case you should look at Elastic 3D beam structures — Numerical tours of continuum mechanics using FEniCS master documentation
A FEniCSx version will be coming soon

1 Like
12

Good afternoon

Actually I have been researching and working on the issue for quite some time now. I apologize if it seemed like I was asking for someone to write a new code from scratch, I thought it was just a few changes from the RectangleMesh example above.

However I’ve created a new thread here with my progress so far including the code to keep things organized and helpful for the community. I would greatly appreciate it if you could take a look and provide any suggestions or guidance.

Thank you for your time and understanding.