dolfinx 0.8.0…
This is the gist example:
import numpy as np
from mpi4py import MPI
import ufl
from dolfinx import mesh, fem, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
import matplotlib.pyplot as plt
L = 120.0 # Beam length in inches
W = 4.0 # Cross-sectional width (and height here; for a square beam)
E = 2000000.0 # Young's modulus in psi
nu = 0.3 # Poisson's ratio
mu = E/(2*(1+nu))
lambda_ = (E * nu) / ((1 + nu) * (1 - 2 * nu))
rho = 30.0/12**3 # Density in lb/in^3 (converted from 30 lb/ft^3)
w_load = -50.0 # lbs/in along the beam length
w_vol = w_load / L / (W*W)
delta = W / L
gamma = 0.4 * delta**2
g = gamma
f_body_total = (0, 0, w_vol + (-rho * g))
domain = mesh.create_box(MPI.COMM_WORLD,
[np.array([0, 0, 0]), np.array([L, W, W])],
[20, 6, 6],
cell_type=mesh.CellType.hexahedron)
V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim, )))
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 = np.array([0, 0, 0], dtype=default_scalar_type)
bc = fem.dirichletbc(u_D,
fem.locate_dofs_topological(V, fdim, boundary_facets),
V)
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_trial = ufl.TrialFunction(V)
v_test = ufl.TestFunction(V)
a = ufl.inner(sigma(u_trial), epsilon(v_test)) * ufl.dx
L_form = ufl.dot(fem.Constant(domain, default_scalar_type(f_body_total)), v_test) * ufl.dx \
+ ufl.dot(T, v_test) * ds
problem = LinearProblem(a, L_form, bcs=[bc],
petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
print("Displacement field (array):", uh.x.array)
print("Minimum displacement:", min(uh.x.array))
This is the second part of the gist that integrates:
x_sections = [0.0, 6.0, 120.0]
print("Cross-sectional positions (inches):", x_sections)
tol = 0.5
bending_moments = []
u = fem.Function(V, name="Displacement")
for x_sec in x_sections:
def facets_at_xsec(x):
return np.isclose(x[0], x_sec, atol=tol)
facets = mesh.locate_entities(domain, 2, facets_at_xsec)
print(facets)
if len(facets) == 0:
print(f"No cells found near x = {x_sec}")
bending_moments.append(0.0)
continue
markers = np.full(len(facets), 1, dtype=np.int32)
cell_tags = mesh.meshtags(domain, domain.topology.dim, facets, markers)
dx_section = ufl.Measure("dx", domain=domain, subdomain_data=cell_tags)
z_c = W / 2
X = ufl.SpatialCoordinate(domain)
z = X[2]
sigma_xx = sigma(uh)[1, 1]
My_integrand = (z - z_c) * sigma_xx
My_form = My_integrand
My = fem.assemble_scalar(fem.form(My_form * ds(dx_section)))
bending_moments.append(My)
print(f"Bending moment M_y at x = {x_sec} is approximately: {My} in-lb")
In regards to the second part that integrates, sigma(uh) has a ufl.shape of (3,)… I have these components to go by as a reference so far:
• σₓₓ = σ(u)[0, 0]
• σₓᵧ = σ(u)[0, 1]
• σₓ𝓏 = σ(u)[0, 2]
• σᵧₓ = σ(u)[1, 0]
• σᵧᵧ = σ(u)[1, 1]
• σᵧ𝓏 = σ(u)[1, 2]
• σ𝓏ₓ = σ(u)[2, 0]
• σ𝓏ᵧ = σ(u)[2, 1]
• σ𝓏𝓏 = σ(u)[2, 2]
Output:
Displacement field (array): [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 2.99149928e-03
-2.21514691e-08 -1.34378060e-01]
Minimum displacement: -0.13437807301956423
Cross-sectional positions (inches): [0.0, 6.0, 120.0]
[ 2 8 14 20 38 44 50 56 86 92 98 104 110 155 161 167 173 179
185 248 254 260 266 272 278 284 372 378 384 390 496 502 508 623 629 752]
Bending moment M_y at x = 0.0 is approximately: 134022.83141157418 in-lb
[ 3 10 16 23 40 47 52 59 88 95 101 106 113 157 164 170 176 181
188 250 257 263 269 275 280 287 375 381 387 393 499 505 511 626 632 755]
Bending moment M_y at x = 6.0 is approximately: 134022.83141157418 in-lb
[1960 1961 1962 1963 2072 2073 2074 2075 2172 2173 2174 2175 2176 2258
2259 2260 2261 2262 2263 2327 2328 2329 2330 2331 2332 2333 2376 2377
2378 2379 2406 2407 2408 2424 2425 2433]
Bending moment M_y at x = 120.0 is approximately: 134022.83141157418 in-lb
So I have the same moment coming out @ each individual section. The ending cell numbers are 752,755,2433 of the tree sections in question…I also tried to use sigma(uh)[1,1] however I have the same difficulty with the same quantity coming out at each section.
Why is the quantity being computed always the same, shouldn’t it be changing at each cross section? Could it be explained how to use the stress components correctly in this scenario to obtain a proper moment at each cross section?
As far as I know the moment at the free end should be about -3000.0 inch/lbs approximately… That’s by the formula w=-50.0 * x=120.0**2 /2… which is the American Wood Council’s version of things where w=lbs/inch x=length in inches…