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

mesh
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