Still trying to learn FEniCS basics, I am trying to apply uniform distributed load on a portion of a 1D cantilever beam in 3D space.
This is how I locate the station :
def load_portion(x, on_boundary):
return on_boundary and start_x <= x[0] <= end_x
… mark the facets :
facets = MeshFunction("size_t", mesh, 1, 0) # MeshFunction for facets
AutoSubDomain(load_portion).mark(facets, 1) # mark 1 for portion_load
… define the differential element of integration :
ds_facets = Measure("ds",
domain=mesh,
subdomain_data=facets,
metadata=metadata)
… and finally apply load :
distributed_load = Constant(-rho*S*g)
l_form = (distributed_load # intensity
*w_[1] # y direction
*ds_facets(1)) # on mark 1 facets
… but it doesn’t match beam theory so far :
maximal deflection = -4.02e-05
beam theory = -7.04e-05
My main source is Static analysis with external distributed mass
Could you please help me spot my mistake ?
My full code here bollow :
# running in Google Colab
try:
import google.colab
except ImportError:
# If not in Colab, import necessary FEniCS packages
import ufl_legacy
import dolfin
else:
# If in Colab, ensure FEniCS is installed
try:
import ufl_legacy
import dolfin
except ImportError:
# Download and install FEniCS if not already installed
!wget "https://fem-on-colab.github.io/releases/fenics-install-real.sh" -O "/tmp/fenics-install.sh" && bash "/tmp/fenics-install.sh"
import ufl_legacy
import dolfin
!pip install meshio # mesh handling
# !apt-get install -qq xvfb libgl1-mesa-glx
# !pip install pyvista -qq
# Importing necessary libraries from FEniCS and others
from dolfin import *
from ufl_legacy import Jacobian, diag
import numpy as np
import meshio # mesh handling
# import pyvista
"""# Mesh"""
# Defining the length of the cantilever beam and the number of elements
length = 10.0
N = 10 # number of elements
# Creating an array of points for the mesh (N+1 points in 3D space)
points = np.zeros((N+1, 3))
points[:, 0] = np.linspace(0, length, N+1)
# Defining the cells (elements) of the mesh
cells = np.array([[i, i+1] for i in range(N)])
# Writing the mesh to an XDMF file using meshio
meshio.write("mesh.xdmf", meshio.Mesh(points, {"line": cells}))
# Reading the mesh from the XDMF file into a FEniCS Mesh object
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
infile.read(mesh)
"""# Sections
Same as [CTV24_FEniCS-elastic-3D-beam-structure-tutorial_2024-05-16.ipynb](https://colab.research.google.com/drive/1PFxZcjQ5UaI9WL-cUiI8h9JzDJt5qgqf#scrollTo=hXViYLk3nLwo)
"""
# Function to calculate the tangent vector of the mesh elements
def tangent(mesh):
t = Jacobian(mesh)
return as_vector([t[0,0], t[1, 0], t[2, 0]])/sqrt(inner(t,t))
# Calculate the tangent vector of the mesh
t = tangent(mesh)
ez = as_vector([0, 0, 1]) # unit vector in the z-direction.
a1 = cross(t, ez) # first cross-sectional vector
a1 /= sqrt(dot(a1, a1)) # normalize
a2 = cross(t, a1)
a2 /= sqrt(dot(a2, a2))
# Values
thick_val = 1
width_val = thick_val
E_val = 1e5 # Young's modulus
rho_val = 1e-3 # Density
g_val = 1 # Gravity
EI1_val = E_val * width_val * thick_val**3 / 12
S_val = thick_val * width_val
# Define section dimensions
thick = Constant(thick_val)
width = Constant(width_val)
# Define material properties
E = Constant(E_val)
nu = Constant(0.3) # Poisson's ratio
G = E / 2 / (1 + nu) # Shear modulus
rho = Constant(rho_val) # Density
g = Constant(g_val) # Gravity
# Define section properties
S = thick * width # Cross-sectional area
ES = E * S # Axial rigidity
EI1 = E * width * thick**3 / 12 # Bending rigidity about axis 1
EI2 = E * width**3 * thick / 12 # Bending rigidity about axis 2
GJ = G * 0.26 * thick * width**3 # Torsional rigidity
kappa = Constant(5. / 6.) # Shear correction factor
GS1 = kappa * G * S # Shear rigidity in direction 1
GS2 = kappa * G * S # Shear rigidity in direction 2
"""# Function Space"""
# Defining the function space for the displacement and rotation fields
Ue = VectorElement("CG", mesh.ufl_cell(), 2, dim=3)
Te = VectorElement("CG", mesh.ufl_cell(), 1, dim=3)
W = FunctionSpace(mesh, MixedElement([Ue, Te]))
"""# Boundaries"""
# Locate clamp
def clamp(x, on_boundary):
return near(x[0], 0.)
# Apply clamp
bc = DirichletBC(W, Constant((0, 0, 0, 0, 0, 0)), clamp)
# Locate right_end
def right_end(x, on_boundary):
return near(x[0], length)
# Locate partial load station
end_frac = 0.8 # as a fraction of length
start_frac = 0.5 # as a fraction of portion_load_end
end_x = end_frac * length
start_x = start_frac * end_x
def load_portion(x, on_boundary):
return on_boundary and start_x <= x[0] <= end_x
# Create a MeshFunction objects to store the marks
boundaries = MeshFunction("size_t", mesh, 0, 0) # MeshFunction for vertices
AutoSubDomain(right_end).mark(boundaries, 1) # mark 1 for right_end
facets = MeshFunction("size_t", mesh, 1, 0) # MeshFunction for facets
AutoSubDomain(load_portion).mark(facets, 1) # mark 1 for portion_load
"""# Functions"""
# Defining the test and trial functions
u_ = TestFunction(W)
du = TrialFunction(W)
(w_, theta_) = split(u_)
(dw, dtheta) = split(du)
# Define generalized strains as a function of displacement and rotation
def tgrad(u):
return dot(grad(u), t)
def generalized_strains(u):
(w, theta) = split(u)
return as_vector([
dot(tgrad(w), t), # Axial strain
dot(tgrad(w), a1) - dot(theta, a2), # Shear strain in direction 1
dot(tgrad(w), a2) + dot(theta, a1), # Shear strain in direction 2
dot(tgrad(theta), t), # Torsional strain
dot(tgrad(theta), a1), # Bending strain about axis 1
dot(tgrad(theta), a2)]) # Bending strain about axis 2
# Define generalized stresses as a function of generalized strains
def generalized_stresses(u):
return dot(diag(as_vector([ES, GS1, GS2, GJ, EI1, EI2])),
generalized_strains(u))
# Defining the variational form
Sig = generalized_stresses(du) # stresses
Eps = generalized_strains(u_) # strain
# will be used to integrate the shear terms
metadata = {"quadrature_scheme": "default", "quadrature_degree": 4}
dx_shear = dx(scheme="default", metadata=metadata)
ds_vertices = Measure("ds",
domain=mesh,
subdomain_data=boundaries,
metadata=metadata)
ds_facets = Measure("ds",
domain=mesh,
subdomain_data=facets,
metadata=metadata)
# Define the bilinear form (stiffness matrix) of the beam
k_form = (
sum([Sig[i] * Eps[i] * dx for i in [0, 3, 4, 5]]) # normal and bending terms
+ (Sig[1] * Eps[1] + Sig[2] * Eps[2]) # shear terms
* dx_shear) # Integrate
"""# Concentrated load only"""
# apply load
f = 1e3 # concentrated load
l_form = (inner(-f, # f intensity
u_[1]) # along y
*ds_vertices(1)) # on mark 1 vertices
# Assemble the equation
K, b = assemble_system(k_form, l_form, bc)
# Solve the equation
u = Function(W)
solve(K, u.vector(), b)
# maximal deflection at the free end of the beam
max_deflection = u(length, 0, 0)[1]
# Calculate value with beam theory formula
beam_theory_deflection = -(f * length**3) / (3 * EI1_val)
# Print result
print("Point load only")
print("maximal deflection = ", max_deflection) # = -40.21200000025746
print("beam theory = ", beam_theory_deflection) # = -40.0
"""# Distributed load only"""
# distributed weight
distributed_load = Constant(-rho*S*g)
l_form = distributed_load*w_[1]*dx
# Assemble the stiffness matrix and load vector
K, b = assemble_system(k_form, l_form, bc)
# Solve for displacements and rotations
u = Function(W, name="Displacement, distributed load")
solve(K, u.vector(), b)
# Extract maximal deflection at the free end of the beam
max_deflection = u(length, 0, 0)[1]
# Value for beam theory
beam_theory_deflection = (-rho_val*S_val*g_val * length**4) / (8 * EI1_val)
# Print result
print("Distributed load only")
print("maximal deflection = ", max_deflection) # = -0.00015106000000100
print("beam theory = ", beam_theory_deflection) # = -0.00015
"""# Combined loads"""
# Concentrated load
l_form_concentrated = (inner(-f, u_[1]) *ds_vertices(1))
# Distributed load
distributed_load = Constant(-rho*S*g)
l_form_distributed = distributed_load*w_[1]*dx
# Combined load
l_form = l_form_concentrated + l_form_distributed
# Assemble the stiffness matrix and load vector
K, b = assemble_system(k_form, l_form, bc)
# Solve for displacements and rotations
u = Function(W, name="Displacement, distributed load")
solve(K, u.vector(), b)
# Extract maximal deflection at the free end of the beam
max_deflection = u(length, 0, 0)[1]
# Value using the beam theory formula
beam_theory_deflection = (
-(f * length**3) / (3 * EI1_val) +
(-rho_val * S_val * g_val * length**4) / (8 * EI1_val))
# Print result
print("Combined Loads")
print("maximal deflection = ", max_deflection) # = -40.21215106025746
print("beam theory = ", beam_theory_deflection) # = -40.00015
"""# Partially Distributed Load"""
# Distributed load on a portion of the beam
distributed_load = Constant(-rho*S*g)
l_form = (distributed_load # intensity
*w_[1] # y direction
*ds_facets(1)) # on mark 1 facets
# Assemble the stiffness matrix and load vector
K, b = assemble_system(k_form, l_form, bc)
# Solve for displacements and rotations
u = Function(W, name="Displacement, distributed load")
solve(K, u.vector(), b)
# Calculating the deflection using the beam theory formula
# https://www.engineersedge.com/beam_bending/beam_bending24.htm
distributed_load_val = -rho_val * S_val * g_val
beam_theory = ( distributed_load_val
*( 4 * length * ( end_x**3 - start_x**3 )
- ( end_x**4 - start_x**4))
/ ( 24 * EI1_val))
# Print the result
print("Partially Distributed Load")
print("maximal deflection = ", u(length, 0, 0)[1]) # = -4.021200000025747e-05
print("beam theory = ", beam_theory) # = -7.04e-05