I am trying to replecate the eigenfrequencies from this paper: Modal Analysis of Single Rectangular Cantilever Plate by Mathematically, FEA and Experimental
Ashish R. Sonawane, Poonam S. Talmale
I cant get to achive results anywhere near that after trying a lot of things.
My resulting frequencies seam quite low leading me to belive that my boundry condition isn’t applied correctly. I create the mesh using gmsh.
Here my python code:
from dolfinx.fem import Function, form, functionspace, dirichletbc, locate_dofs_topological
from dolfinx.fem.petsc import assemble_matrix
from slepc4py import SLEPc
from ufl import dx, grad, inner, TrialFunction, TestFunction
from dolfinx.io import XDMFFile, gmshio
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
from ufl import sym, Identity, div
import os
PART = "STEEL_BEAM" # Name of the part being analyzed (used for file naming)
GEO_FILE = "~/Dev/BA/simulation/" + PART + ".geo"
MSH_FILE = "~/Dev/BA/simulation/" + PART + ".msh"
# Material properties
E = 70e9 # Young's modulus [Pa]
nu = 0.35 # Poisson's ratio [-]
rho = 2700 # Density [kg/m³]
# create mesh
try:
os.system("gmsh " + GEO_FILE + " -3 -o " + MSH_FILE)
except Exception:
print("GMSH Geometry (.geo) file could not be converted to a Mesh (.msh) file")
quit()
print("Mesh creation finished.\n")
# mesh import
mesh, cell_tags, facet_tags = gmshio.read_from_msh(PART + ".msh", MPI.COMM_WORLD, 0, gdim=3)
# test and trial functions
deg = 2
V = functionspace(mesh, ("CG", deg, (3,))) # Vector function space for displacement
u = TrialFunction(V)
v = TestFunction(V)
# Lame parameters
mu = E / (2 * (1 + nu))
lambda_ = E * nu / ((1 + nu) * (1 - 2 * nu))
# Strain and stress
def epsilon(u):
return sym(grad(u))
def sigma(u):
return lambda_ * div(u) * Identity(3) + 2 * mu * epsilon(u)
# boundary conditions
fixed_end = facet_tags.find(1) # Assuming region 1 corresponds to the fixed end
u_zero = np.zeros(3, dtype=PETSc.ScalarType) # Zero displacement
bc = dirichletbc(u_zero, locate_dofs_topological(V, 2, fixed_end), V)
# matrix assemblys
k_form = form(inner(sigma(u), epsilon(v)) * dx)
m_form = form(rho * inner(u, v) * dx)
K = assemble_matrix(k_form, [bc])
M = assemble_matrix(m_form, [bc])
K.assemble()
M.assemble()
# setup eigensolver
solver = SLEPc.EPS().create()
solver.setDimensions(4) # number of eigenvalues to compute
solver.setProblemType(SLEPc.EPS.ProblemType.GHEP) # generalized Hermitian eigenvalue problem
st = SLEPc.ST().create()
st.setType(SLEPc.ST.Type.SINVERT)
st.setShift(0.1)
st.setFromOptions()
solver.setST(st)
solver.setOperators(K, M)
# solve the problem
solver.solve()
xr, xi = K.createVecs()
tol, maxit = solver.getTolerances()
nconv = solver.getConverged()
print(f"\nNumber of iterations: {solver.getIterationNumber()}")
print(f"Solution method: {solver.getType()}\n")
print(f"Stopping condition: tol={tol:.4g}, maxit={maxit}\n")
eig_vector = []
eig_freq = []
if nconv > 0:
print("Eigenvalues:")
for i in range(nconv):
k = solver.getEigenpair(i, xr, xi)
fn = np.sqrt(k.real) / (2 * np.pi)
eig_freq.append(fn)
print(f"Mode {i}: {fn} Hz")
vect = xr.getArray()
eig_vector.append(vect.copy())
My geo file for gmsh is:
// ============================================
// Manual Method: Define Points, Lines, Surfaces
// ============================================
mesh_size = 0.05; // Mesh size
l = 1; // Length [m] x-direction
b = 0.05; // Width [m] z-direction
h = 0.005; // Height [m] y-direction
Point(1) = {0, 0, 0, mesh_size};
Point(2) = {0, h, 0, mesh_size};
Point(3) = {0, h, b, mesh_size};
Point(4) = {0, 0, b, mesh_size};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Transfinite Curve {2, 4} = 6 Using Progression 1;
Transfinite Curve {1, 3} = 4 Using Progression 1;
Curve Loop(1) = {1, 2, 3, 4};
Plane Surface(1) = {1};
Recombine Surface {1}; // Quad mesh base surface (1)
Transfinite Surface {1};
Extrude {1, 0, 0} {
Surface{1};
Layers{50};
Recombine;
}
Coherence;
Physical Surface("FixedEnd") = {1};
Physical Volume("SteelBeam") = {1};
// Set mesh options
Mesh.Algorithm = 8; // Frontal-Delaunay for quads
Mesh.MshFileVersion = 2;
// Generate mesh
Mesh 3;
My resulting eigenfrequencies are:
Eigenvalues:
Eigenvalues:
Mode 0: 0.15915494309189526 Hz
Mode 1: 0.15915494309189532 Hz
Mode 2: 0.15915494309189532 Hz
Mode 3: 0.15915494309189537 Hz
Mode 4: 0.1591549430918954 Hz
Mode 5: 0.15915494309189818 Hz
Mode 6: 4.137811871898824 Hz
Mode 7: 25.930992886274172 Hz
Mode 8: 41.14547178779196 Hz
Mode 9: 72.63890011999023 Hz
Mode 10: 142.46413270523485 Hz
For the boundry condition I also tried:
def clamped_boundary(x):
return np.isclose(x[0], 0.0) # e.g., x=0 plane
fixed-end = locate_entities_boundary(mesh, 2, clamped_boundary)
The expected natural freuquencies are (experimental from linked paper):
1: 0.233
2: 1.142
3: 1.737
4: 4.799
5: 8.32
6: 9.36
7: 11.67
8: 15.013
9: 23.617
10: 23.944
With specefied values:
L=1000mm
b=50mm
h=5mm
E=70e3 N/mm^2
nu=0.35
rho=2700kg/m^3