How to obtain the moment distribution of a beam fixed at both ends

General
3 

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

L = 120.0  # Extending the domain to the point of interest
H = 8.0
Nx = 40  # Increase the number of elements to maintain mesh quality
Ny = 40

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

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)

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)

fx = 0.0
fy = -0.1
f = fem.Constant(domain, (fx, fy))

msh_cell_type = domain.basix_cell()
Ve = basix.ufl.element(
    'Lagrange', msh_cell_type, degree=2, shape=(domain.geometry.dim,))
V = fem.functionspace(domain, Ve)

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

def left(x):
    return np.isclose(x[0], 0.0)

fdim = domain.topology.dim - 1
marked_values = []
marked_facets = []
facets = mesh.locate_entities_boundary(domain, fdim, left)
marked_facets.append(facets)
marked_values.append(np.full_like(facets, 1))
marked_facets = np.hstack(marked_facets)
marked_values = np.hstack(marked_values)
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(
    domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets]
)

ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tag)
left_dofs = fem.locate_dofs_geometrical(V, left)
u_bc = fem.Function(V)
bcs = [fem.dirichletbc(u_bc, left_dofs)]

u = fem.Function(V, name="Displacement")

problem = dolfinx.fem.petsc.LinearProblem(a, l, bcs, u=u)
problem.solve()

# Moment of inertia for the cross-section
I = (4 * 8**3) / 12

# Initialize bending moment
Mz_total = 0.0

# Loop through sections of the beam to compute the cumulative bending moment
for x_pos in np.linspace(0, L, Nx+1):
    def at_section(x):
        return np.isclose(x[0], x_pos)

    section_facets = mesh.locate_entities_boundary(domain, fdim, at_section)
    
    if section_facets.size > 0:
        marked_section = np.hstack(section_facets)
        section_tag = mesh.meshtags(
            domain, fdim, marked_section, np.full_like(marked_section, 2)
        )
        section_ds = ufl.Measure("ds", domain=domain, subdomain_data=section_tag, subdomain_id=2)
        Mz_section = fem.assemble_scalar(fem.form(-ufl.SpatialCoordinate(domain)[1] * sigma(u)[0, 0] * section_ds))
        Mz_total += Mz_section

# Multiply by the moment of inertia
Mz_total *= I

print(f"Bending moment Mz at x={L} = {Mz_total:.6f}")

# Calculate the theoretical maximum bending moment for comparison
theoretical_max_moment = fy * L * L / 2
print('Theoretical max moment:', theoretical_max_moment)

So far I have dolfinx showing more moment for the distributed load.

I did try to make use of:


#x = ufl.SpatialCoordinate(domain)
#n = ufl.FacetNormal(domain)
#My = fem.assemble_scalar(fem.form(ufl.cross(x, sigma(u)*n)[1]*ds(1)))

however I got a long traceback. Since that didn’t work I took a guess at some method that I thought may have a chance at matching the 720.0 figure of the American Wood Council. There still is a pretty good sized discrepancy.

Have I done anything at all correct here? Could it be shown how to come close to 720.0 as the AWC says, and maybe some details why there is a difference right now of about 200 inch/lbs of moment, and maybe how to remedy that and do this correctly?