How to impose point load at free end for cantilever beam (linear elasticity)?

Question

Followed example/suggestion to implement point load at specified point on the free end of cantilever beam as on Cantilever beam example (linear elasticity) . However, as mentioned in reply #2 on this link, I get u = 0 at all nodes. Can you please help me to locate the mistake/suggest correction?

For code see reply #2 of Cantilever beam example (linear elasticity) .

Comments

Please format your code properly by indenting it four spaces.

---

Provided a web-link for the code to avoid clutter/formatting issues.

---

Hello! I have written a code for 3D cantilever beam with specified point load at the free end. I am novice user of Fenics (1.0.0 on Windows) and would like to know if it is implemented correctly. Can any one check/help?

from dolfin import *
import numpy as np   

# Load mesh and define function space
d0 = 8.
d1 = 1.
d2 = 1.
n0 = int(d0*6)
n1 = int(d1*6)
n2 = int(d2*6)

# Mesh for cantilever beam
Omega = Box(0., 0., 0., d0, d1, d2, n0, n1, n2)

Omega.init()

V = VectorFunctionSpace(Omega, "CG", 1)

tol = 1.e-8 # tolerance 

# Left boundary 
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)

# Right boundary
class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 8.0)      

# Point (8,0.5,0.5) at which load is added
class point(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0],8.) and near(x[1],0.5,1e-2) and near(x[2],0.5,1e-2)

 #Boundary segments
left = Left()
pt = point()
right = Right()

# Initialize mesh function for boundary domains
boundaries = FacetFunction("uint", Omega)
boundaries.set_all(0)
left.mark(boundaries, 1)
pt.mark(boundaries, 2)      

 # Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant((0.0, 0.0, 00.0))

# Elasticity parameters
E = 210000.0
nu = 0.3

mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))

# Strain
def eps(v):
    return sym(grad(v))

# Stress
def sigma(v):
    return 2.0*mu*sym(grad(v)) + lmbda*tr(sym(grad(v)))*Identity(v.cell().d) 

# Strain energy 
def energy(v):
    return lmbda/2.0*(tr(eps(v)))^2 + mu*tr(eps(v)**2) 

# Define new measures associated with the interior domains and exterior boundaries
    ds = Measure("ds")[boundaries]

# Define variational problem
g_T= Constant((0.0, 0.0, -0.3))

a = inner(grad(v),sigma(u))*dx 

L= dot(v,g_T)*ds(2)

# Boundary condition
left_end_displacement = Constant((0.0, 0.0, 0.0))

bc = [DirichletBC(V, left_end_displacement, left)] 

# Compute solution
u = Function(V)
solve(a == L, u, bc);

# Save solution to VTK format
File("Output_elas3D_PointLoad.pvd", "compressed") << u
plot(u,interactive=True);

# Strain energy 
DG = FunctionSpace(Omega, 'CG', 1)
U = inner(sigma(u),eps(u))
E =project(U,DG)
E_AtVertex = E.vector().array()

#Element wise U
DG = FunctionSpace(Omega,"DG",0)
v = TestFunction(DG)
Vol =1   # CellVolume(Omega) gives error (??)
U_elementWise = assemble(U*v/Vol*dx).array() 

---

Hi, Did you solve the problem finally? I am stuck on the same thing.

As per below listed suggestions, when I apply pointwise Drichilet BC on the point, I get a displacement "resulting" from the boundary condition in the viper plot.

So, if the BC is (0,-7,0), the max displacement in the viper plot is near 7. So I feel that I am implying displacement instead of load.

Also, what ChrishZhang says about unchanging displacement w.r.t. elastic modulus but changing displacement w.r.t. point load is correct. Kristian says their is some kind of normalising function running in the background, is that correct?

If your code is verified and working now, pls. post it here to help end this thread. Thanks a lot!

Answer 1

Did you try

method="pointwise"

for

DirichletBC

?

Comments

This works for me in dolfin version 1.3.0

from dolfin import *
import numpy as np   

# Load mesh and define function space
d0 = 8.
d1 = 1.
d2 = 1.
n0 = int(d0*6)
n1 = int(d1*6)
n2 = int(d2*6)

# Mesh for cantilever beam
Omega = BoxMesh(0., 0., 0., d0, d1, d2, n0, n1, n2)

Omega.init()

V     = VectorFunctionSpace(Omega, "CG", 1)

tol = 1.e-8 # tolerance 

# Left boundary 
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)

# Right boundary
class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 8.0)      

# Point (8,0.5,0.5) at which load is added
class point(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0],8.) and near(x[1],0.5,0.) and near(x[2],0.5,0.)

#Boundary segments
left = Left()
pt = point()
right = Right()

# Initialize mesh function for boundary domains
boundaries = FacetFunction("uint", Omega)
boundaries.set_all(0)
left.mark(boundaries, 1)
pt.mark(boundaries, 2)      

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant((0.0, 0.0, 00.0))

# Elasticity parameters
E = 210000.0
nu = 0.3

mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))

# Strain
def eps(v):
    return sym(grad(v))

# Stress
def sigma(v):
    return 2.0*mu*sym(grad(v)) + lmbda*tr(sym(grad(v)))*Identity(v.cell().d) 

# Strain energy 
def energy(v):
    return lmbda/2.0*(tr(eps(v)))^2 + mu*tr(eps(v)**2) 

# Define new measures associated with the interior domains and exterior boundaries
ds = Measure("ds")[boundaries]

# Define variational problem
g_T= Constant((0.0, 0.0, -0.3))

a = inner(grad(v),sigma(u))*dx 

L= dot(v,g_T)*ds(2)

# Boundary condition
left_end_displacement = Constant((0.0, 0.0, 0.0))

bc = [DirichletBC(V, left_end_displacement, left),DirichletBC(V, g_T, pt, method="pointwise")] 

# Compute solution
u = Function(V)
solve(a == L, u, bc)

# Save solution to VTK format
File("Output_elas3D_PointLoad.pvd", "compressed") << u
plot(u,interactive=True)

# Strain energy 
#DG = FunctionSpace(Omega, 'CG', 1)
#U = inner(sigma(u),eps(u))
#E =project(U,DG)
#E_AtVertex = E.vector().array()
#
##Element wise U
#DG = FunctionSpace(Omega,"DG",0)
#v = TestFunction(DG)
#Vol =1   # CellVolume(Omega) gives error (??)
#U_elementWise = assemble(U*v/Vol*dx).array() 

---

I run your code.
When I modified the Young's modulus, it seems that the displacement field remained the same.
That really confused me.

---

Changing the Young's modulus rescales the displacements. The plotting function probably has some normalisation functionality that hides this.

---

not really, the displacemence has changed when the point load increase.

---

I never said otherwise.

---

Hi,
your answer is not correct. When you write:

DirichletBC(V, g_T, pt, method="pointwise")] 

you are imposing a fixed displacement g_T = -3.0, since these are Dirichlet conditions. It is not a point load. This is the reason why when ChrisZang modified the Young's modulus, he got the same displacement.

Also you wrote

L= dot(v,g_T)*ds(2)

so it seems you have a point load, but you can check that the contribution of L to the solution is zero. You can multiply L by any constant, you will get the same solution. This is because the boundary pt = point() is just a point, so when assembling the contribution of L is zero.

Finally, try to remove the Dirichlet condition

DirichletBC(V, g_T, pt, method="pointwise")] 

and your solution will be uniformly zero, which means you don't have any point load.

I mentionned in another question that I also have the same problem, i.e. I don't know how to apply a force only at a point. It seems one should use the function PointSource, but I could not make it work for vector functions. The only solution I found was to define a small Neumann boundary, but not reduced to a point, for instance in your case, if you redefine

class point(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and  abs(x[0]-8.0) < 0.2 and abs(x[1]-0.5) < 0.2  and abs(x[1]-0.5) < 0.2 

and you remove the Dirichlet condition

DirichletBC(V, g_T, pt, method="pointwise")] 

then it is working.

Still this solution is not very satisfying because the Neumann boundary is not a point, and I would be glad if someone can explain how to use PointSource with a vector-valued function.