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

General
1 

from dolfinx import mesh, fem, default_scalar_type
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 = 160
Ny = 5

# 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.1
fy = 50.0
f = fem.Constant(domain, (fx, -fy / (L * W)))

# Theoretical maximum moment
theoretical_moment_max = (((-fy / L) * (L**2)) / (2))

# 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 left(x):
    return np.isclose(x[0], 0.0)

fdim = domain.topology.dim - 1
marked_values = []
marked_facets = []
# Concatenate and sort the arrays based on facet indices
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]
)


# 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, left)
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
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)))

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 = {-L * H * fy})")
print("-" * 50)

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

print(f"Theoretical M max: {theoretical_moment_max:.6f} inch\lbs")



--------------------------------------------------
Bending moment Mz = -449.826670
        (analytic = 2880000.0)
--------------------------------------------------


Theoretical M max: -3000.000000 inch\lbs

So far the Mz I am getting is nowhere near -3000.00 lbs/inch however I really am not quire sure yet…

So far I have these formulas to go by where w is pounds each inch of material and l is the length in inches of the cantilever beam.

Screenshot 2024-08-01 180937
Screenshot 2024-08-01 180937702×405 22.6 KB

Why this discrepancy between Mz and theoretical moment?

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…

3

I still don’t quite understand why I had to do this:


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

however it seems that doing so is the only way I can get a value to come out close to the theoretical…

If one runs the post-processed stress example for a beam of these dimensions, it seems that things state:


Bending moment Mz = -301.784223
        (analytic = -14400.0)

I don’t really understand it why -14400.0 or what it means so far, as far as I know the theoretical maximum @ L should be -1800±…

I guess I am realizing at this point that I was able to get an Mz on or near that however there still seems to be some difficulties in code to work through…

As much as I had attempted to add another set of facet markers for the free end I was not able to get that working… that is how things had got so in code there…

Ideally I would like to be able to traverse the L from x=0 to tox=L at the segments and see the Mz increase as things approach L…

Anyone can post this thread how to clean up those problems and get things working a bit better than this?