Taking moment from a point off a 3D beam

General
1 

import pyvista
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

# Geometry and material properties
L = 120.0
W = 4.0
H = 8.0

rho = 0.48
delta = W / L
gamma = 0.4 * delta**2
beta = 1.25

E = 2000000  # Modulus of elasticity
nu = 0.3  # Poisson's ratio

lambda_ = E * nu / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))

g = gamma

# Create the mesh
domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([L, W, H])],
                         [20, 6, 6], cell_type=mesh.CellType.hexahedron)
V = fem.functionspace(domain, ("Lagrange", 1, (domain.topology.dim, )))

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

fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)

u_D = fem.Constant(domain, np.array([0, 0, 0], dtype=default_scalar_type))
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets), V)

# Define free end facets for calculating moment
def free_end(x):
    return np.isclose(x[0], L)

free_end_facets = mesh.locate_entities_boundary(domain, fdim, free_end)

# Create a MeshTags object to mark the free end
facet_markers = mesh.meshtags(domain, fdim, free_end_facets, np.full(len(free_end_facets), 1, dtype=np.int32))
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_markers)

# Define the elasticity problem
def epsilon(u):
    return ufl.sym(ufl.grad(u))  # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)

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

fx = 0.1 
fy = -rho * g * 0.5
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(domain, default_scalar_type((0, 0.0, fy)))
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx

problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

# Calculate bending moment Mz
x = ufl.SpatialCoordinate(domain)
Mz = fem.assemble_scalar(fem.form(-x[1] * sigma(uh)[0, 0] * ds(1)))
print(f"Bending moment Mz = {Mz:.6f}")

This is a 3D derivation of the code found here:

Computing consistent reaction forces — Computational Mechanics Numerical Tours with FEniCSx (bleyerj.github.io)

What I am trying to accomplish is at some point on or near the end of the cantilever hopefully towards the bottom center to take the Mz, which in this case is an already known quantity.

In terms of:


Mz = fem.assemble_scalar(fem.form(-x[1] * sigma(uh)[0, 0] * ds(1)))

I get the feeling that ds(1)) may need to change so far now that things moved from 2D to 3D. OK… -x[1] is the co-ordinates of the Y axis of the domain, really what I want is the moment at a specific point. So far there is no errors to report however the Mz seems to be coming out to -0.00 which is an obvious misgnomer.

How can the moment about the Z be taken from a point somewhere near the bottom center of the cantilever?

2

P.S. Just by the way I am looking for things to come out somewhere near Bending moment Mz = -1728.174038:

Finding the moment of a cantilever beam from its lengthwise cross section - General - FEniCS Project

So far I was able to get some scalar type values from facets but they didn’t yet amount to the whole moment at the free end…

3

I was able to make some improvements to the script so far. There is some scalar value of output for the Mz of -473±, however that seems to me much too small so far to be within an acceptable tolerance to what solid mechanics may predict.


import pyvista
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

# Geometry and material properties
L = 120.0
W = 4.0
H = 8.0

rho = 0.48
delta = W / L
gamma = 0.4 * delta**2
beta = 1.25

E = 2000000  # Modulus of elasticity
nu = 0.3  # Poisson's ratio

lambda_ = E * nu / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))

g = gamma

# Create the mesh
domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([L, W, H])],
                         [20, 60, 60], cell_type=mesh.CellType.hexahedron)
V = fem.functionspace(domain, ("Lagrange", 1, (domain.topology.dim, )))
# Define boundary conditions
def clamped_boundary(x):
    return np.isclose(x[0], 0)

fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)

u_D = fem.Constant(domain, np.array([0, 0, 0], dtype=default_scalar_type))
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets), V)

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

# Find facets on the boundary where x = L (free end)
fdim = domain.topology.dim - 1
# Ensure that the facets at the free end (x = L) are correctly marked
free_end_facets = mesh.locate_entities_boundary(domain, fdim, free_end)
facet_markers = mesh.meshtags(domain, fdim, free_end_facets, np.full(len(free_end_facets), 2, dtype=np.int32))

ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_markers)

# Define the elasticity problem
def epsilon(u):
    return ufl.sym(ufl.grad(u))  # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)

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

fx = 0.1 
fy = -rho * g * 0.25
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(domain, default_scalar_type((0, fy, 0.0)))
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx

problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

# Calculate bending moment Mz
x = ufl.SpatialCoordinate(domain)
Mz = fem.assemble_scalar(fem.form(-x[1] * sigma(uh)[0, 0] * ds()))
print(f"Bending moment Mz = {Mz:.6f}")
# Bending moment Mz = -473.749420

The density of the mesh so far is approaching the limitations of what I can accomplish locally on my own system. I noticed that if I increase density much further I get a message output to the terminal that states, “killed”, (more information about the message would be great!).

I don’t quite understand it at this point why -473, that seems to be a large departure from the 2D sample which had looked close enough for all intensive purposes to what solid mechanics may predict…

I noticed so far that the in terms of the 2D samples that the virtual work method does seem to come closer to what solid mechanics may predict. At this point I don’t really have a great handle on the virtual work method and how to employ the use of it outside on the context of the 2D sample that is shown…

Anyone can post this thread, why the Mz that is output here -473± seems to be roughly half of what is predicted in the 2D sample I provided a link to here on this thread is stating?