Hi all
I want to use FEnics to calculate sailing ship rigging.
Starting with basic problem I struggle to have a simple cantilever beam example working.
Especially having the concentrated loads case to match the beam theory.
I have been through these for inspiration:
- Cantilever example with beam element for "CG" order 2 space fails
- How can point loads be applied to the end of the cantilever beam?
- Single 2D Truss mesh
- Linear elasticity with beam element
- No solution for simple cantilever beam with beam element
- Elastic 3D beam structures — Numerical tours of continuum mechanics using FEniCS master documentation
- Linear elastic 2D truss — Computational Mechanics Numerical Tours with FEniCSx
- Packages — FEM on Colab
Although I learnt from them a general understanding of the code, I am still having trouble getting my example to work for concentrated loads. I would greatly appreciate it if someone could help me figure out what I am doing wrong.
Here is the code I have so far:
# Attempt to import google.colab to check if 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
# Install meshio for mesh handling
!pip install meshio
# Importing necessary libraries from FEniCS and others
from dolfin import *
from ufl_legacy import Jacobian, diag
import numpy as np
import meshio
"""# Mesh"""
# Defining the length of the cantilever beam and the number of elements
length = 10.0
N = 80 # 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))
# Define section dimensions
thick = Constant(1)
width = thick#/3
# Define material properties
E = Constant(1e5) # Young's modulus
nu = Constant(0.3) # Poisson's ratio
G = E / 2 / (1 + nu) # Shear modulus
rho = Constant(1e-3) # Density
g = Constant(1) # 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
"""# Functions"""
# 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]))
# 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
dx_shear = dx(scheme="default",metadata={"quadrature_scheme":"default", "quadrature_degree": 4})
# 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
)
"""# Boundaries"""
# Define boundary condition for clamping (fixed end at x=0))
def clamp(x, on_boundary):
return near(x[0], 0.)
# Apply Dirichlet boundary condition
bc = DirichletBC(W, Constant((0, 0, 0, 0, 0, 0)), clamp)
"""# Concentrated load only"""
# Define concentrated load at the end of the beam
P = 1e3 # Point load in N
T = PointSource(W.sub(1), Point(length, 0, 0), -P)
# No distributed loads
l = Constant(0) * w_[1] * dx
# Assemble the stiffness matrix and load vector
K = assemble(k_form)
b = assemble(l)
T.apply(b) # Apply the point load to the load vector
bc.apply(K) # Apply boundary conditions to the stiffness matrix
bc.apply(b) # Apply boundary conditions to the load vector
# Solve for displacements and rotations
u_pt = Function(W, name="Displacement, point load")
solve(K, u_pt.vector(), b)
# Print maximal deflection at the free end of the beam
print("Maximal deflection, point load only:", -u_pt(length, 0, 0)[1])
# Calculating the deflection using the beam theory formula
beam_theory_deflection = (P * length**3) / (3 * EI1)
# Print the result
print("Beam theory deflection:", float(beam_theory_deflection))
# Maximal deflection, point load only: 6.000000000035433
# Beam theory deflection: 40.0
"""# Distributed load only"""
# distributed weight
l = Constant(-rho*S*g)*w_[1]*dx
# Assemble the stiffness matrix and load vector
K = assemble(k_form)
b = assemble(l)
bc.apply(K) # Apply boundary conditions to the stiffness matrix
bc.apply(b) # Apply boundary conditions to the load vector
# Solve for displacements and rotations
u_dist = Function(W, name="Displacement, distributed load")
solve(K, u_dist.vector(), b)
# Print maximal deflection at the free end of the beam
print("Maximal deflection, distributed load only:", u_dist(length, 0, 0)[1])
# Calculating the deflection using the beam theory formula for distributed load
beam_theory_deflection_distributed = (-rho*S*g * length**4) / (8 * EI1)
# Print the result
print("Beam theory deflection for distributed load:", float(beam_theory_deflection_distributed))
"""# Combined loads"""
# Define concentrated load at the end of the beam
P = 1e3 # Point load in N
T = PointSource(W.sub(1), Point(length, 0, 0), -P)
# Define distributed weight
l = Constant(-rho*S*g)*w_[1]*dx
# Assemble the stiffness matrix and load vector
K = assemble(k_form)
b = assemble(l)
T.apply(b) # Apply the point load to the load vector
bc.apply(K) # Apply boundary conditions to the stiffness matrix
bc.apply(b) # Apply boundary conditions to the load vector
# Solve for displacements and rotations
u = Function(W, name="Displacement, combined load")
solve(K, u.vector(), b)
# Print maximal deflection at the free end of the beam
print("Maximal deflection, combined loads:", -u(length, 0, 0)[1])
