Hey everyone,
I’m currently trying to run a modal analysis for a simple 3D elastic body using FEniCSx and SLEPc, but I’m having trouble understanding why my results do not make physical sense.
I’m trying to solve the Eigenvalue problem:
K*u = omega**2*M*u
I want to solve it using the SLEPc PEP solver and I also tested it with the EPS solver for comparison. In the future i would like to also add a damping matrix C, it is currently set to zero in the PEP solver. To my understanding, the results should be similar between the two Solvers, but they are different and not feasable.
When I increase the geometry, the eigenfrequencies increase while they should decrease.
The geometry should normally be t = 2e-3, h = 35e-3, 50e-3 and the first natural frequency should be around 1400 Hz based on experiments, but my solver returns a Eigenvalue of about 1.0.
I’m pretty sure my Dirichlet boundary conditions are correct, the metal piece should be fixed at the bottom.
I checked the norms and the symmetry of my matrices, they seem about fine, but feel free to double check.
Does anyone know, where my problem could be? Should the EPS and PEP solvers give identical results when the damping matrix is set to zero, or am I already mistaken there?
I would be happy to share my code snippet with you, if that helps. Thanks for any hints or sanity checks!
- Robin
from mpi4py import MPI
import numpy as np
import ufl
from dolfinx import mesh, fem
import dolfinx
import dolfinx.fem.petsc
from slepc4py import SLEPc
#Define Object
t = 2e+3
h = 35e+3
w = 50e+3
domain = mesh.create_box(
MPI.COMM_WORLD,
[[0.0, 0.0, 0.0], [w, h, t]],
[12,8,2],
cell_type=mesh.CellType.hexahedron,
)
tdim = domain.topology.dim
#functionspace
V = fem.functionspace(domain, ("Lagrange", 2, (domain.geometry.dim,)))
#Material Data
rho = 2.7e3 # kg/m^3
young = 71e9 # Pa
nu = 0.30 # Poisson
mu = young / (2.0 * (1.0 + nu))
lmbda = young * nu / ((1.0+nu)*(1.0-2.0*nu))
#Dirichlet
def bottom(x):
return np.isclose(x[1], 0.0)
clamped_facets = dolfinx.mesh.locate_entities_boundary(domain, tdim - 1, bottom)
clamped_nodes = dolfinx.mesh.locate_entities(domain, 0, bottom)
num_facets = domain.topology.index_map(tdim - 1).size_local
markers = np.zeros(num_facets, dtype=np.int32)
clamped = 1
markers[clamped_facets] = clamped
facet_marker = dolfinx.mesh.meshtags(domain, tdim - 1, np.arange(num_facets, dtype=np.int32), markers)
clamped_dofs = dolfinx.fem.locate_dofs_topological(V, facet_marker.dim, facet_marker.find(clamped))
u_clamped = dolfinx.fem.Constant(domain, (0.0, 0.0, 0.0))
bcs = [
dolfinx.fem.dirichletbc(u_clamped, clamped_dofs, V)
]
def epsilon(u):
return ufl.sym(ufl.grad(u))
def sigma(u):
return lmbda*ufl.tr(epsilon(u))*ufl.Identity(tdim) + 2.0*mu*epsilon(u)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
k = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
m = rho * ufl.inner(u,v) * ufl.dx
k_compiled = dolfinx.fem.form(k)
K = dolfinx.fem.petsc.assemble_matrix(k_compiled, bcs=bcs)
K.assemble()
m_compiled = dolfinx.fem.form(m)
M = dolfinx.fem.petsc.assemble_matrix(m_compiled, bcs=bcs)
M.assemble()
#Damping
alpha = 0.0
beta = 0.0
C = K.copy()
C.scale(beta)
C.axpy(alpha, M, True)
C.assemble()
#Solve
#Pep
pep = SLEPc.PEP().create(domain.comm)
pep.setOperators([M,C,K])
pep.setProblemType(SLEPc.PEP.ProblemType.GENERAL)
pep.setType(SLEPc.PEP.Type.QARNOLDI)
pep.setTolerances(tol=1e-8, max_it=200)
pep.setFromOptions()
pep.solve()
#eps
eps = SLEPc.EPS().create(domain.comm)
eps.setOperators(M, K)
eps.setProblemType(SLEPc.EPS.ProblemType.GHEP)
pep.setTolerances(tol=1e-8, max_it=200)
pep.setFromOptions()
eps.solve()
#Analyse Output
nconv = pep.getConverged()
print("converged eigenpairs:", nconv)
for i in range(min(nconv, 10)):
Eigenvalue = pep.getEigenpair(i)
omega_re = np.real(Eigenvalue)
omega_im = np.imag(Eigenvalue)
omega = abs(omega_im)
freq = np.sqrt(omega) / (2*np.pi)
zeta = -omega_re / np.sqrt(omega_re**2 + omega_im**2)
if omega_re < -1e-10:
print(f"Mode {i}: freq = {freq:.3f} Hz, damping ratio = {zeta:.4f}")
elif omega_re > 1e-10:
print(f"Mode {i}: freq = {freq:.3f} Hz, damping ratio = {zeta:.4f}")
else:
print(f"Mode {i} : freq = {freq:.3f} Hz, damping ratio = {zeta:.4f}")
nconv = eps.getConverged()
print("converged eigenpairs:", nconv)
for i in range(min(nconv, 10)):
Eigenvalue = eps.getEigenpair(i)
omega = Eigenvalue
freq = np.sqrt(omega) / (2*np.pi)
print(f"Mode {i}: freq = {freq:.3f} Hz")