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

General
1

I have set up a mock-up for the uniformly distributed load:


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

# Define domain dimensions and material properties
L = 120 #inches
W = 4 #inches
H = 8 # iches 
rho = 0.0624 #lb per ft ^3
delta = W / L
gamma = 0.4 * delta**2

# Young's modulus and Poisson's ratio
# Lames parameters
E = 2000000 # PSI
nu = 0.0
mu = E / (2 * (1 + nu))
lambda_ = E * nu / ((1 + nu) * (1 - 2 * nu))

# Gravity
g = gamma

w = 6.66 # lbs per foot

I = (W * H**3) / 12

# Create a denser mesh and define function space
domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([L, W, H])],
                         [40, 12, 12], cell_type=mesh.CellType.hexahedron)
V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim, )))

# Define boundary conditions
def clamped_boundary_left(x):
    return np.isclose(x[0], 0)

def clamped_boundary_right(x):
    return np.isclose(x[0], L)

fdim = domain.topology.dim - 1
boundary_facets_left = mesh.locate_entities_boundary(domain, fdim, clamped_boundary_left)
boundary_facets_right = mesh.locate_entities_boundary(domain, fdim, clamped_boundary_right)

u_D = np.array([0, 0, 0], dtype=default_scalar_type)
bc_left = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets_left), V)
bc_right = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets_right), V)

# Define the load and weak form
T = fem.Constant(domain, default_scalar_type((0, 0, 0)))
ds = ufl.Measure("ds", domain=domain)

def epsilon(u):
    return ufl.sym(ufl.grad(u))

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

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(domain, default_scalar_type((0, 0, -w / (W * H))))  # Distributed load in the z-direction
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds

# Solve the linear problem
problem = LinearProblem(a, L, bcs=[bc_left, bc_right], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

# Calculate the moment as the curl of the displacement field
#V_curl = fem.FunctionSpace(domain, ("Lagrange", 1))
V_curl = fem.functionspace(domain, ("Lagrange", 1, (3,)))  # Vector function space for 3D vectors
# Define a function to store the curl
curl_uh = fem.Function(V_curl)

# Project the curl onto the function space
w = ufl.TestFunction(V_curl)
v_curl = ufl.TrialFunction(V_curl)
a_curl = ufl.inner(v_curl, w) * ufl.dx
L_curl = ufl.inner(ufl.curl(uh), w) * ufl.dx

# Solve the projection problem
curl_problem = LinearProblem(a_curl, L_curl)
curl_uh = curl_problem.solve()

# Compute the moment
moment = curl_uh

# Extract the moment values, I guess they should be in lbs per inch...
moment_values = moment.x.array.reshape((-1, 3))
max_moment = np.max(moment_values[:, 2])  # Assuming we are interested in the z-component of the moment

# Extract the displacement values
displacement_values = uh.x.array.reshape((-1, 3))
max_displacement = np.max(np.linalg.norm(displacement_values, axis=1))

# Print the results
print(f'Maximum moment (z-component): {max_moment}')
print(f'Maximum displacement: {max_displacement}')

print('Solution and moment distribution were obtained.')


# Analytical results for comparison
awc_max_displacement = (5 * 6.66 * 120.0**4) / (384 * 2000000 * I)
awc_max_moment = (6.66 * 120.0**2) / 8

# Check against the American Wood Council formula worksheet for a 
# uniformly distributed load...
print(f'AWC max displacement: {awc_max_displacement}')
print(f'AWC max moment: {awc_max_moment}')

The moment values coming out of dolfinx are not the same as what the AWC formulas say. The displacement value is closer however it is still far enough away that I would have to say some adjustment is needed. Ideally I would like to have the values coming out of dolfinx to be as precisely close to what the AWC says as possible. I understand that might mean adjusting some parameters. So far I must have some type of incorrect setup for obtaining moment I would guess.

‘’’

Maximum moment (z-component): 6.547203128877776e-15
Maximum displacement: 0.010254790287438496
Solution and moment distribution were obtained.
AWC max displacement: 0.05268164062499999
AWC max moment: 11988.0

‘’’

Can I get some assistence to get these output values to line up?

((wx)/2)(l-x) is the induced moment I thin, in inch lbs, where x is a distance along the length from 0. If it is not possible to get the moment out can you show me how to induce it in addition to w? Does moment have to be set or can it just be gotten somehow by the curl method?

2

There is no way the bending moment can be computed as the curl of a displacement field…
You are using a 3D elastic code, the resulting moment over a surface S around the origin is given by:
\boldsymbol{M} = \int_S \boldsymbol{x} \times \boldsymbol{\sigma n} dS
which contains 3 values : a torsion component and two bending components. To have the bending moment in direction 1, you have to compute something like

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

if the surface S is indexed by surf_idx.

but this will not be strictly in equilibrium with the loading. To have consistent reaction force, I suggest you follow this procedure:

https://bleyerj.github.io/comet-fenicsx/tips/computing_reactions/computing_reactions.html#second-method-using-the-work-of-internal-forces

which can be easily extended to 3D.

1 Like
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?

4

Why do you multiply by the moment of inertia ?
Please note that fy is a volume distributed force so your theoretical formula is not correct. You also work in 2D so you need to multiply integrals by the beam width in the third direction

1 Like
5

OK… In terms of multiplying by the moment of inertia it was a wrong guess at this point.

So you are telling me to multiply by the third direction, so if w=4, h=8 and z=120.0 (length, L), I believe so far you mean the z.

So by multiplying by the third direction I interpret that so far to mean to multiply by the distance between segments in the Z direction that are taken.

I appreciate you, trying to make up for some lost time… Very helpful link, I hope to have a much better understanding of it. Sometimes I need a bit of extra time to look things over.