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
Nz = 121
Nx = 13
Ny = 13
# Create the mesh
domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([L, W, H])],
[121, 13, 13], 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(3) + 2 * mu * epsilon(u)
fx = 0.1
self_w = -rho * g * 0.25
w = -0.25
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(domain, default_scalar_type((0.0, -0.00011, 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()
'''
https://bleyerj.github.io/comet-fenicsx/tips/computing_reactions/computing_reactions.html#second-method-using-the-work-of-internal-forces
LaTeX:::
\[
M_z = \int_{x=0} (\vec{OM} \times \vec{T}) \times \vec{e}_z dS = -\int_{x=0} (y\sigma_x_x) dS = f_{j_x} - \frac{L^2}{2} = H
\]
'''
# 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}")
theoretical_moment = (w*120*120) / 2
print('theoretical_moment:::', -theoretical_moment)
'''
Bending moment Mz = -1779.556528
theoretical_moment::: -1800.0
'''
So far my intent is to set up a distributed load case where, small case cursive w is equal to 0.25 lbs per inch. The reason why I have modified the Y-axis force of the the force vector f to -0.00011, is that is the only way to get the Mz to come out anywhere close to the theoretical moment…
It seems so far that:
def free_end(x):
return np.isclose(x[0], L)
and:
def free_end(x):
return np.isclose(x[0], 0)
are producing the same results for Mz, however I would expect that they should be different if what I am intending to do is working then Mz should decrease the closer that the free_end facets are to the u_D boundary facets…
f = fem.Constant(domain, default_scalar_type((0.0, -0.00011, 0.0))).
Why is the force having to be scaled down so much to achieve a matching Mz and theoretical moment?
Could I get some help to set things up and get them working correctly for a distributed load case described and computed by, print(‘theoretical_moment:::’, -theoretical_moment)?
Maybe this can help to address it, at least one of the questions I am asking:::
import numpy as np
from mpi4py import MPI
import ufl
from dolfinx import fem, mesh, io
import dolfinx.fem.petsc
from dolfinx.fem.petsc import LinearProblem
import matplotlib.pyplot as plt
from matplotlib.tri import Triangulation
# 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))
# --- Identify unique X-values in the mesh ---
# Get the mesh points (geometry of the domain)
points = domain.geometry.x
# Extract unique X-values
unique_x_values = np.unique(points[:, 0]) # Get unique X values (sorted automatically)
#print("Unique X-values in the mesh:", unique_x_values)
# 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 clamped_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, clamped_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, clamped_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)
# --- Stress Tensor Calculation ---
# Define a function space for stress (DG0 element for discontinuous fields)
W = fem.functionspace(domain, ("DG", 0, (2, 2))) # Tensor space
# Compute stress tensor field from displacement
stress = fem.Function(W, name="Stress")
stress_expr = fem.Expression(sigma(u), W.element.interpolation_points())
stress.interpolate(stress_expr)
# Extract the sigma_xx component
stress_values = stress.x.array # Extract all values from the Function array
# Since stress is a 2x2 tensor, we want to extract the sigma_xx component
# which is the first component of each tensor (index 0)
# The tensor has a structure like [[sigma_xx, sigma_xy], [sigma_yx, sigma_yy]]
# sigma_xx is stored in stress_values at positions corresponding to the first component
stress_xx_values = stress_values[::4] # Take every fourth value starting from 0 (for sigma_xx)
# --- Print stress tensor values at vertices ---
# Get the mesh points
points = domain.geometry.x
# --- Print stress tensor values at vertices where x = 120.0 ---
target_x = 120.0
parray =[]
# Iterate over mesh points and print stress only for points where x = 120.0
for i, point in enumerate(points):
if np.isclose(point[0], target_x):
print(f"Stress at point {point}: {stress_xx_values[i]}")
parray.append([[point], [stress_xx_values[i]]])
moment = 0
for i, parray_item in enumerate(parray):
point = parray_item[0][0]
sigmav = parray_item[1][0]
sumv = sigmav * point[1] * 4.0
moment+=sumv
print('moment:', moment)
NX: 301
Ny: 13
Horizontal reaction Rx = -18.905485
(analytic = -0.0)
--------------------------------------------------
Vertical reaction Ry = -8.102351
(analytic = -120.0)
--------------------------------------------------
Bending moment Mz = -1728.174038
(analytic = 1800.0)
--------------------------------------------------
Stress at point [120. 3.38461538 0. ]: 49.673508529903074
Stress at point [120. 4. 0.]: 65.37257034690427
Stress at point [120. 2.76923077 0. ]: 35.54221486740928
Stress at point [120. 2.15384615 0. ]: 46.50190175001872
Stress at point [120. 1.53846154 0. ]: 15.527938604478791
Stress at point [120. 0.92307692 0. ]: 13.583349708020325
Stress at point [120. 0.30769231 0. ]: 21.751232787833867
Stress at point [120. -0.30769231 0. ]: 31.233044324898543
Stress at point [120. -0.92307692 0. ]: 38.58220388722755
Stress at point [120. -1.53846154 0. ]: 6.376374255143272
Stress at point [120. -2.15384615 0. ]: 15.95482729120869
Stress at point [120. -2.76923077 0. ]: 24.719743311336074
Stress at point [120. -3.38461538 0. ]: -7.883703141590987
Stress at point [120. -4. 0.]: 3.8237744532083573
moment: 2099.415246637295
Here the same beam is mocked-up in 2D. One thing that I noticed is that in terms of using the stress tensor to integrate the moment of force that just small changes to Nx or Ny seems to throw off the computations pretty wildy. Anyone can shed some light on this, why that changing Nx,Ny at all seem to throw off the computations so far?
I had to really jimmy jack the Nx,Ny to achieve the output moment: 2099.415246637295, which does seem a little bit high albeit in a type of range I personally find acceptable in terms of matching the analytic... The reason for using moment+=sumv, is because that is the only way I was able to find, to walk the length of the beam and output the moment at a slice at some location x...Anyone can shed some light on this why 2099.41 is coming out just a little bit too high there?
Here in the 2D sample stress_xx_values values or used, so it seems that when integrating the moment from the stress tensor in 3D space that there is a 3x3 stress tensor to be using which is different from the one shown...
things have to be L there as the clamped end because, as much as I had tried to use 0.0 for the clamped end, I found that in doing so that the computations where no longer good… Anyone can shed some light on this, why the clamped end doesn’t seem to want to walk over to 0.0 where I wanted it originally?
So in terms of doing the same thing in 3D space, it seems there is a different type of stress tensor with different dimensions, however at this point I am still a little bit left out as to how exactly the integration of the stress tensor can or should take place in 3D realm… This is the question I am asking here:
So far I am finding with 3D as with 2D the value for Mz albeit not correct at all yet is also seems to swing fairly wildly depending on Nx,Nz of the mesh, but I don't really understand yet how to set one up so that Mz can come out in a more stable type of way or if there is some logic behind it which Nx, Nz to choose to have the best results for integrating stress tensor for moment to come out...
If there are any answers to these questions and how to relate this 2D sample to doing the same integration for moment in 3D, I would really appreciate to hear about…