Help solving a clasical cantiveler beam in comparition with theory

General
1

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
nthreads = 1
os.environ["OMP_NUM_THREADS"] = str(nthreads)
os.environ["OPENBLAS_NUM_THREADS"] = str(nthreads)
os.environ["MKL_NUM_THREADS"] = str(nthreads)

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

# loads
load = (10, 0, 0)
# --------------------#


# ------- #
# 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.add(name)
    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()

    fixed_sf = model.addPhysicalGroup(2, [circle], tag=FIXED_TAG)
    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])
    tip_sf = model.addPhysicalGroup(2, [other_boundary_entities[-1]], tag=TIP_SURFACE_TAG)
    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])
    vol = model.addPhysicalGroup(3, volume_entities, tag=VOLUME_TAG)
    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,)))

load_facets = facet_markers.find(TIP_SURFACE_TAG)
mt = mesh.meshtags(domain, 2, load_facets, 1)
metadata = {"quadrature_degree": 4}
ds = ufl.Measure("ds", subdomain_data=mt, subdomain_id=1, metadata=metadata)
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 #
# --------------------#
f = ufl.as_vector(load)

# 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):
    return ufl.sym(ufl.grad(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
Fx, Fy, Fz = load


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.

2

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)

is a surface load and not a resultant force, so your Fx should be instead load[0]*surface where surface is the cross-section area.

1 Like
3

@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
nthreads = 1
os.environ["OMP_NUM_THREADS"] = str(nthreads)
os.environ["OPENBLAS_NUM_THREADS"] = str(nthreads)
os.environ["MKL_NUM_THREADS"] = str(nthreads)

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

# loads
load = (10, 0, 0)
# --------------------#


# ------- #
# 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.add(name)
    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()

    fixed_sf = model.addPhysicalGroup(2, [circle], tag=FIXED_TAG)
    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])
    tip_sf = model.addPhysicalGroup(2, [other_boundary_entities[-1]], tag=TIP_SURFACE_TAG)
    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])
    vol = model.addPhysicalGroup(3, volume_entities, tag=VOLUME_TAG)
    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,)))

load_facets = facet_markers.find(TIP_SURFACE_TAG)
mt = mesh.meshtags(domain, 2, load_facets, 1)
metadata = {"quadrature_degree": 4}
ds = ufl.Measure("ds", subdomain_data=mt, subdomain_id=1, metadata=metadata)
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 #
# --------------------#
f = ufl.as_vector(load)

# 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):
    return ufl.sym(ufl.grad(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.

4

Back here to post the solution, at least for the direct integration method, The problem with the orientation of the reaction forces is that my fixed surface is normal to the Z axis so the correct formulation is

R_x &= \int_S -\sigma_{zx} \, dS \\
R_y &= \int_S -\sigma_{zy} \, dS \\
R_z &= \int_S -\sigma_{zz} \, dS

So in this case the solution is:

Rx = fem.assemble_scalar(fem.form(-sigma(u)[2,0] * ds(1)))
Ry = fem.assemble_scalar(fem.form(-sigma(u)[2,1] * ds(1)))
Rz = fem.assemble_scalar(fem.form(-sigma(u)[2,2] * ds(1)))

But for a general solution (remembering this is MWE) we need to use the traction vector which is $$\mathbf{T} &= \sigma \cdot \mathbf{n}$$ where n is the surface normal vector. So we can do:

x, y, z = ufl.SpatialCoordinate(domain)
n = ufl.FacetNormal(domain)
Tx, Ty, Tz = ufl.dot(sigma(u), n)
Rx = fem.assemble_scalar(fem.form(-Tx * ds(1)))
Ry = fem.assemble_scalar(fem.form(-Ty * ds(1)))
Rz = fem.assemble_scalar(fem.form(-Tz * ds(1)))

Also I have the wrong torque orientation in my analytical form, so I use the general torque form too:

x, y, z = 0, 0, L
Mx = y * Fz - z * Fy
My = z * Fx - x * Fz
Mz = x * Fy - y * Fx

And now everything is in the right place. I still don’t know how to implement the torque in the virtual work, but for now, this is enough.