Finding the moment of a cantilever beam from its lengthwise cross section

General
2 

from dolfinx import mesh, fem
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import numpy as np

# Define the beam properties
L = 120.0  # Length in inches
H = 8.0    # Height in inches
W = 4.0    # Width in inches
I = (W * H**3) / 12  # Moment of inertia in inches^4

# Mesh properties
Nx = 301  # Finer mesh along the length
Ny = 13  # Finer mesh along the height

print('NX:', Nx)
print('Ny:', Ny)

# Create the mesh
domain = mesh.create_rectangle(
    MPI.COMM_WORLD,
    [(0.0, -H / 2), (L, H / 2)],
    (Nx, Ny),
    diagonal=mesh.DiagonalType.crossed,
)

# Material properties
E = fem.Constant(domain, 2000000.0)
nu = fem.Constant(domain, 0.3)
mu = (E / (2 * (1 + nu)))
lmbda = (E * nu) / ((1 + nu) * (1 - 2 * nu))

# Define strain and stress
def eps(v):
    return ufl.sym(ufl.grad(v))

def sigma(v):
    return lmbda * ufl.tr(eps(v)) * ufl.Identity(2) + 2.0 * mu * eps(v)

# Applied load
fx = 0.0
fy = 0.25 * W
fy_analytic = 0.25
f = fem.Constant(domain, (fx, -fy))  # Uniformly distributed load over the length

# Define function space
V = fem.functionspace(domain, ("P", 2, (2,)))

# Define variational problem
du = ufl.TrialFunction(V)
u_ = ufl.TestFunction(V)
a = ufl.inner(sigma(du), eps(u_)) * ufl.dx
l = ufl.inner(f, u_) * ufl.dx

# Boundary conditions
def free_end(x):
    return np.isclose(x[0], L)

fdim = domain.topology.dim - 1

# Locate facets on the left (fixed end)
left_facets = mesh.locate_entities_boundary(domain, fdim, free_end)
marked_values = np.full_like(left_facets, 1)
facet_tag = mesh.meshtags(domain, fdim, left_facets, marked_values)

# Measure for boundary integration
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tag)

# Apply Dirichlet BCs (fixing the left end)
left_dofs = fem.locate_dofs_geometrical(V, free_end)
u_bc = fem.Function(V)
bcs = [fem.dirichletbc(u_bc, left_dofs)]

# Solve the problem
u = fem.Function(V, name="Displacement")
problem = LinearProblem(a, l, bcs, u=u)
problem.solve()

# Calculate reaction forces and moments at the fixed end
x = ufl.SpatialCoordinate(domain)
Rx = fem.assemble_scalar(fem.form(-sigma(u)[0, 0] * ds(1)))
Ry = fem.assemble_scalar(fem.form(-sigma(u)[1, 1] * ds(1)))
Mz = fem.assemble_scalar(fem.form(-x[1] * sigma(u)[0, 0] * ds(1)))
Mz*=W

# Print results
print(f"Horizontal reaction Rx = {Rx:.6f}")
print(f"             (analytic = {-L * H * fx})")
print("-" * 50)

print(f"Vertical reaction Ry = {Ry:.6f}")
print(f"           (analytic = {-fy * L})")
print("-" * 50)

print(f"Bending moment Mz = {Mz:.6f}")
print(f"        (analytic = {(fy_analytic * L**2) / 2})")
print("-" * 50)



Horizontal reaction Rx = -18.905485
             (analytic = -0.0)
--------------------------------------------------
Vertical reaction Ry = -8.102351
           (analytic = -120.0)
--------------------------------------------------
Bending moment Mz = -1728.174038
        (analytic = 1800.0)
--------------------------------------------------

OK… The anallytic and the bending moment Mz seem to be psuedo similar. I followed the instructions I was given in the other post:

How to obtain the moment distribution of a beam fixed at both ends - General - FEniCS Project

This is the best I could do so far. It seems that “fy” will already be distributed evenly. I would guess that may be a bit much and would need to be toned down in terms of loading.

Have I done this OK? Is there anything I have missed? So far I would propose this a potential solution to the problem of finding the beam’s moment in this cantilever case…My apologies, it took some time to sort through everything and arrive at this potential answer so far…