# Help solving a clasical cantiveler beam in comparition with theory

I’m solving a classical problem with clamped cantilever beam 3D and a force on its end. The beam has a circular section. I have this code:

``````import os
from typing import NamedTuple

import gmsh
import numpy as np
import ufl
from dolfinx import fem, io, mesh
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
from petsc4py import PETSc

# enforce 1 thread per core

comm = MPI.COMM_WORLD

# -------------------#
# Problem properties #
# -------------------#
class Material(NamedTuple):
name: str
E: PETSc.ScalarType
nu: PETSc.ScalarType
rho: PETSc.ScalarType

material = Material("Aluminum 6061", 68.9e9, 0.33, 2700)

# beam size
height = 150  # m
diameter = 6.5  # m

# --------------------#

# ------- #
# MESHING #
# ------- #
FIXED_TAG = 1000
TIP_SURFACE_TAG = 3000
VOLUME_TAG = 4000
ELEMENTS_ORDER = 2

def gmsh_tower(model: gmsh.model, name: str, h: float, d: float) -> gmsh.model:
model.setCurrent(name)

# Recombine tetrahedra to hexahedra
gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2)
gmsh.option.setNumber("Mesh.RecombineAll", 2)
gmsh.option.setNumber("Mesh.CharacteristicLengthFactor", 0.1)

circle = model.occ.addDisk(0, 0, 0, d / 2, d / 2)
extruded_geometry = model.occ.extrude([(2, circle)], 0, 0, h, numElements=[h], recombine=True)
model.occ.synchronize()

model.setPhysicalName(2, fixed_sf, "FIXED SURFACE")
boundary_entities = model.getEntities(2)
other_boundary_entities = []
for entity in boundary_entities:
if entity != circle:
other_boundary_entities.append(entity[1])
model.setPhysicalName(2, tip_sf, "TIP SURFACE")

model.mesh.generate(3)
model.mesh.setOrder(ELEMENTS_ORDER)
volume_entities = []
for entity in extruded_geometry:
if entity[0] == 3:
volume_entities.append(entity[1])
model.setPhysicalName(3, vol, "Mesh volume")
return model

def create_mesh(comm: MPI.Comm, model: gmsh.model, name: str, filename: str, mode: str):
msh, ct, ft = io.gmshio.model_to_mesh(model, comm, rank=0)
msh.name = name
ct.name = f"{msh.name}_cells"
ft.name = f"{msh.name}_facets"
with io.XDMFFile(msh.comm, filename, mode) as file:
msh.topology.create_connectivity(2, 3)
file.write_mesh(msh)
file.write_meshtags(
ct, msh.geometry, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry"
)
file.write_meshtags(
ft, msh.geometry, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry"
)
return (msh, ct, ft)

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)

# Create model
model_name = "Tower"
model = gmsh.model()
model = gmsh_tower(model=model, name=model_name, h=height, d=diameter)
model.setCurrent(model_name)
mesh_file = f"out_gmsh/mesh_rank_{MPI.COMM_WORLD.rank}.xdmf"
domain, cell_markers, facet_markers = create_mesh(
comm=MPI.COMM_SELF,
model=model,
name=model_name,
filename=mesh_file,
mode="w",
)

# -------------- #
# Function Space #
# -------------- #
dim = domain.geometry.dim
V = fem.functionspace(domain, ("Lagrange", ELEMENTS_ORDER, (dim,)))

mt = mesh.meshtags(domain, 2, load_facets, 1)
x = ufl.SpatialCoordinate(domain)

E = fem.Constant(domain, float(material.E))
nu = fem.Constant(domain, float(material.nu))
rho = fem.Constant(domain, float(material.rho))

# --------------------#
# Boundary conditions #
# --------------------#

# Fixed BC
fixed_facets = facet_markers.find(FIXED_TAG)
fixed_surface_dofs = fem.locate_dofs_topological(V, 2, fixed_facets)

u_bc = fem.Function(V)
fixed_bc = fem.dirichletbc(u_bc, fixed_surface_dofs)

# -------------------------#
# linear elastic equations #
# -------------------------#

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

def epsilon(u):

def sigma(u):
return lmbda * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)

du = ufl.TrialFunction(V)
u_ = ufl.TestFunction(V)
a = ufl.inner(sigma(du), epsilon(u_)) * ufl.dx
l = ufl.inner(f, u_) * ds(1)
u = fem.Function(V, name="Displacement")

bcs = [fixed_bc]

problem = LinearProblem(
a, l, u=u, bcs=bcs, petsc_options={"ksp_type": "cg", "pc_type": "gamg", "ksp_rtol": 1e-8}
)

uh = problem.solve()

R = diameter / 2
L = height

print("====================================================================")
print("                           DEFORMATION                              ")
print("====================================================================")
Ix = np.pi / 4 * R**4
print(f"         Analytic: {(Fx * L**3) / (3 * material.E * Ix)}")
print(f"         Computed: {u.x.array.max()}")

print("====================================================================")
print("                           REACTIONS                                ")
print("====================================================================")
Rx = -np.pi * R**2 * Fx
Ry = -np.pi * R**2 * Fy
Rz = -np.pi * R**2 * Fz
print(f"    Analytic: Rx = {Rx:.3f}, Ry = {Ry:.3f}, Rz = {Rz:.3f}")

# ====================================================================
Rx = fem.assemble_scalar(fem.form(-sigma(u)[0, 0] * ds))
Ry = fem.assemble_scalar(fem.form(-sigma(u)[0, 1] * ds))
Rz = fem.assemble_scalar(fem.form(-sigma(u)[0, 2] * ds))
print(f"      Direct: Rx = {Rx:.3f}, Ry = {Ry:.3f}, Rz = {Rz:.3f}")

# ====================================================================
residual = ufl.action(a, u) - l
v_reac = fem.Function(V)
virtual_work_form = fem.form(ufl.action(residual, v_reac))

def one(x):
values = np.zeros((1, x.shape[1]))
values[0] = 1.0
return values

u_bc.sub(0).interpolate(one)
fem.set_bc(v_reac.vector, bcs)
Rx = fem.assemble_scalar(virtual_work_form)

u_bc.vector.set(0.0)
v_reac.vector.set(0.0)
u_bc.sub(1).interpolate(one)
fem.set_bc(v_reac.vector, bcs)
Ry = fem.assemble_scalar(virtual_work_form)

u_bc.vector.set(0.0)
v_reac.vector.set(0.0)
u_bc.sub(2).interpolate(one)
fem.set_bc(v_reac.vector, bcs)
Rz = fem.assemble_scalar(virtual_work_form)

print(f"Virtual Work: Rx = {Rx:.3f}, Ry = {Ry:.3f}, Rz = {Rz:.3f}")

``````

with the solution:

``````====================================================================
DEFORMATION
====================================================================
Analytic: 1.8634166653601339e-06
Computed: 6.168997972153805e-05
====================================================================
REACTIONS
====================================================================
Analytic: Rx = -331.831, Ry = -0.000, Rz = -0.000
Direct: Rx = 0.022, Ry = 0.000, Rz = -335.679
Virtual Work: Rx = -331.815, Ry = -0.000, Rz = -0.000
``````

Soo seems that there is something wrong in the maximum displacement computation, as well as the “Direct” method to compute the reaction forces on the beam base, because its giving an approximate good results but in the wrong axis.

Also I want to compute the reaction bending moments, but I remove from here due to is really wrong results.

Just had a quick look so it needs more checking but it seems that you apply a load of intensity Fx distributed over the beam cross-section, the load in

``````l = ufl.inner(f, u_) * ds(1)
``````

1 Like

@bleyerj Thanks you are rigth that solve my deformation distortion. Thanks a lot.

Here I post the new version of my code. Also I add the Torques computations, but as I said before their differ too much from the analytics results.

``````
import os
from typing import NamedTuple

import gmsh
import numpy as np
import ufl
from dolfinx import fem, io, mesh
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
from petsc4py import PETSc

# enforce 1 thread per core

comm = MPI.COMM_WORLD

# -------------------#
# Problem properties #
# -------------------#
class Material(NamedTuple):
name: str
E: PETSc.ScalarType
nu: PETSc.ScalarType
rho: PETSc.ScalarType

material = Material("Aluminum 6061", 68.9e9, 0.33, 2700)

# beam size
height = 150  # m
diameter = 6.5  # m

# --------------------#

# ------- #
# MESHING #
# ------- #
FIXED_TAG = 1000
TIP_SURFACE_TAG = 3000
VOLUME_TAG = 4000
ELEMENTS_ORDER = 2

def gmsh_tower(model: gmsh.model, name: str, h: float, d: float) -> gmsh.model:
model.setCurrent(name)

# Recombine tetrahedra to hexahedra
gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2)
gmsh.option.setNumber("Mesh.RecombineAll", 2)
gmsh.option.setNumber("Mesh.CharacteristicLengthFactor", 0.1)

circle = model.occ.addDisk(0, 0, 0, d / 2, d / 2)
extruded_geometry = model.occ.extrude([(2, circle)], 0, 0, h, numElements=[h], recombine=True)
model.occ.synchronize()

model.setPhysicalName(2, fixed_sf, "FIXED SURFACE")
boundary_entities = model.getEntities(2)
other_boundary_entities = []
for entity in boundary_entities:
if entity != circle:
other_boundary_entities.append(entity[1])
model.setPhysicalName(2, tip_sf, "TIP SURFACE")

model.mesh.generate(3)
model.mesh.setOrder(ELEMENTS_ORDER)
volume_entities = []
for entity in extruded_geometry:
if entity[0] == 3:
volume_entities.append(entity[1])
model.setPhysicalName(3, vol, "Mesh volume")
return model

def create_mesh(comm: MPI.Comm, model: gmsh.model, name: str, filename: str, mode: str):
msh, ct, ft = io.gmshio.model_to_mesh(model, comm, rank=0)
msh.name = name
ct.name = f"{msh.name}_cells"
ft.name = f"{msh.name}_facets"
with io.XDMFFile(msh.comm, filename, mode) as file:
msh.topology.create_connectivity(2, 3)
file.write_mesh(msh)
file.write_meshtags(
ct, msh.geometry, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry"
)
file.write_meshtags(
ft, msh.geometry, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry"
)
return (msh, ct, ft)

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)

# Create model
model_name = "Tower"
model = gmsh.model()
model = gmsh_tower(model=model, name=model_name, h=height, d=diameter)
model.setCurrent(model_name)
mesh_file = f"out_gmsh/mesh_rank_{MPI.COMM_WORLD.rank}.xdmf"
domain, cell_markers, facet_markers = create_mesh(
comm=MPI.COMM_SELF,
model=model,
name=model_name,
filename=mesh_file,
mode="w",
)

# -------------- #
# Function Space #
# -------------- #
dim = domain.geometry.dim
V = fem.functionspace(domain, ("Lagrange", ELEMENTS_ORDER, (dim,)))

mt = mesh.meshtags(domain, 2, load_facets, 1)
x = ufl.SpatialCoordinate(domain)

E = fem.Constant(domain, float(material.E))
nu = fem.Constant(domain, float(material.nu))
rho = fem.Constant(domain, float(material.rho))

# --------------------#
# Boundary conditions #
# --------------------#

# Fixed BC
fixed_facets = facet_markers.find(FIXED_TAG)
fixed_surface_dofs = fem.locate_dofs_topological(V, 2, fixed_facets)

u_bc = fem.Function(V)
fixed_bc = fem.dirichletbc(u_bc, fixed_surface_dofs)

# -------------------------#
# linear elastic equations #
# -------------------------#

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

def epsilon(u):

def sigma(u):
return lmbda * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)

du = ufl.TrialFunction(V)
u_ = ufl.TestFunction(V)
a = ufl.inner(sigma(du), epsilon(u_)) * ufl.dx
l = ufl.inner(f, u_) * ds(1)
u = fem.Function(V, name="Displacement")

bcs = [fixed_bc]

problem = LinearProblem(
a, l, u=u, bcs=bcs, petsc_options={"ksp_type": "cg", "pc_type": "gamg", "ksp_rtol": 1e-8}
)

uh = problem.solve()

R = diameter / 2
L = height
A = np.pi * R**2
Fx, Fy, Fz = np.array(load) * A

print("====================================================================")
print("                           DEFORMATION                              ")
print("====================================================================")
Ix = np.pi / 4 * R**4
print(f"         Analytic: {(Fx * L**3) / (3 * material.E * Ix)}")
print(f"         Computed: {u.x.array.max()}")

print("====================================================================")
print("                           REACTIONS                                ")
print("====================================================================")
Rx = -Fx
Ry = -Fy
Rz = -Fz
print(f"    Analytic: Rx = {Rx:.3f}, Ry = {Ry:.3f}, Rz = {Rz:.3f}")

# ====================================================================
Rx = fem.assemble_scalar(fem.form(-sigma(u)[0, 0] * ds(1)))
Ry = fem.assemble_scalar(fem.form(-sigma(u)[0, 1] * ds(1)))
Rz = fem.assemble_scalar(fem.form(-sigma(u)[0, 2] * ds(1)))
print(f"      Direct: Rx = {Rx:.3f}, Ry = {Ry:.3f}, Rz = {Rz:.3f}")

# ====================================================================
residual = ufl.action(a, u) - l
v_reac = fem.Function(V)
virtual_work_form = fem.form(ufl.action(residual, v_reac))

def one(x):
values = np.zeros((1, x.shape[1]))
values[0] = 1.0
return values

def _y(x):
values = np.zeros((1, x.shape[1]))
values[0] = x[1]
return values

def _z(x):
values = np.zeros((1, x.shape[1]))
values[0] = x[2]
return values

u_bc.sub(0).interpolate(one)
fem.set_bc(v_reac.vector, bcs)
Rx = fem.assemble_scalar(virtual_work_form)

u_bc.vector.set(0.0)
v_reac.vector.set(0.0)
u_bc.sub(1).interpolate(one)
fem.set_bc(v_reac.vector, bcs)
Ry = fem.assemble_scalar(virtual_work_form)

u_bc.vector.set(0.0)
v_reac.vector.set(0.0)
u_bc.sub(2).interpolate(one)
fem.set_bc(v_reac.vector, bcs)
Rz = fem.assemble_scalar(virtual_work_form)

print(f"Virtual Work: Rx = {Rx:.3f}, Ry = {Ry:.3f}, Rz = {Rz:.3f}")

print("====================================================================")
print("                             TORQUES                                ")
print("====================================================================")
print(f"    Analytic: Mx = {L*Fy:.3f}, My = {L*Fz:.3f}, Mz = {L*Fx:.3f},")
# ====================================================================

Mx_form = fem.form((x[1] * sigma(u)[1, 2] - x[2] * sigma(u)[1, 1]) * ds(1))
My_form = fem.form((x[2] * sigma(u)[2, 0] - x[0] * sigma(u)[2, 2]) * ds(1))
Mz_form = fem.form((x[0] * sigma(u)[0, 1] - x[1] * sigma(u)[0, 0]) * ds(1))

Mx = fem.assemble_scalar(Mx_form)
My = fem.assemble_scalar(My_form)
Mz = fem.assemble_scalar(Mz_form)
print(f"      Direct: Mx = {Mx:.3f}, My = {My:.3f}, Mz = {Mz:.3f}")

# ====================================================================
u_bc.vector.set(0.0)
v_reac.vector.set(0.0)
u_bc.sub(1).interpolate(_z)
fem.set_bc(v_reac.vector, bcs)
Mx = fem.assemble_scalar(virtual_work_form)

u_bc.vector.set(0.0)
v_reac.vector.set(0.0)
u_bc.sub(0).interpolate(_z)
fem.set_bc(v_reac.vector, bcs)
My = fem.assemble_scalar(virtual_work_form)

u_bc.vector.set(0.0)
v_reac.vector.set(0.0)
u_bc.sub(0).interpolate(_y)
fem.set_bc(v_reac.vector, bcs)
Mz = fem.assemble_scalar(virtual_work_form)

print(f"Virtual Work: Mx = {Mx:.3f}, My = {My:.3f}, Mz = {Mz:.3f}")
``````

Solution:

``````====================================================================
DEFORMATION
====================================================================
Analytic: 6.183389012461246e-05
Computed: 6.168997972160277e-05
====================================================================
REACTIONS
====================================================================
Analytic: Rx = -331.831, Ry = -0.000, Rz = -0.000
Direct: Rx = 0.022, Ry = 0.000, Rz = -335.679
Virtual Work: Rx = -331.815, Ry = -0.000, Rz = -0.000
====================================================================
TORQUES
====================================================================
Analytic: Mx = 0.000, My = 0.000, Mz = 49774.609,
Direct: Mx = 2.511, My = 50291.685, Mz = 0.002
Virtual Work: Mx = -0.000, My = -0.000, Mz = -0.003
``````

And clearly I’m missundestanding something with the orientations of the direct method. And the torques in the Virtual Work.