Why do I obtain ~185 near-zero eigenvalues in a clamped 3D elasticity eigenproblem (DOLFINx + SLEPc)?

General
, ,
1

Hello everyone,

I am solving a 3D linear elasticity eigenvalue problem in DOLFINx for the dynamic case of a simple cylindrical rod. The geometry is imported from a STEP file (exported from COMSOL), meshed in Gmsh, and then used in DOLFINx. One end of the rod is fully clamped and the other end is free.

I solve the generalized eigenvalue problem

Ku + \omega^2M u =0

using SLEPc with the solver Krylov–Schur, defining the problem as a GHEP problem. The eigenfrequencies and mode shapes look physically correct. However, the first physical eigenfrequency appears at eigenpair number 186, meaning that I obtain around 185 near-zero eigenvalues before it.

Since the rod is clamped at one end, I would expect no rigid-body modes at all. Even for a fully free body, I would expect only 6 rigid-body modes, not 185.

The solver reports successful convergence, and I have tried changing the tolerance, trying different solvers and refining the mesh. Changing the tolerance to 10^{-14} gave about 265 near-zero eigenvalues, and decreasing the tolerance to 10^{-2} still gave 185 near-zero eigenvalues.

Below is a minimum working example which prints out the eigenfrequencies starting from the first non zero frequency (filtering is done after the eigenvalue solver), as well as plotting the modes from the first non-zero frequency onward. When running the code, the created mesh will appear and the gmsh-window has to be closed before being able to continue running the code. The STEP file “Test_periodic_rod_geo_r025” and the help function file containing the function for building a Gmsh mesh from the STEP file are available in this git-hub repo:
The repository is on GitHub under:
NadineTab02 / FEniCSx_Project

import numpy as np
import matplotlib.pyplot as plt
from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import mesh, fem, plot
from ufl import  dx, inner
from basix.ufl import element
from dolfinx import default_scalar_type
from dolfinx.fem import form, functionspace, Function
from dolfinx.fem.petsc import assemble_matrix
from slepc4py import SLEPc
from dolfinx.io import gmsh as gmshio
import ufl 
import pyvista as pv
from function import import_CAD_to_FEniCSx

#--------------Create mesh from the step file--------------------

lc = 10 # characteristic mesh element size
mesh_info, submesh_info = import_CAD_to_FEniCSx('Test_periodic_rod_geo_r025', 'mesh_simple_case_01', lc)
mesh, cell_tags, facet_tags = mesh_info
dim = mesh.geometry.dim
mesh.geometry.x[:] *= 0.001   # Need to downscale the geometry to original size 1 m

# ---------------Define space function and element----------------------------

deg = 2 
Ve = element("Lagrange", mesh.basix_cell(), deg, shape=(dim,))  
V = functionspace(mesh, Ve)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)



#-------------Material parameters--------------
E = 2.1e14  
rho = 7850
nu =0.3
mu = E / (2*(1 + nu))
lambda_ = E*nu / ((1 + nu)*(1 - 2*nu))

# ------------Boundary tagging------------------------------

facet_tags = mesh_info[2]  #  imported facet_tags
left_tag  = 2001  # start of rod (surface), clamped
right_tag = 1001  # end of rod (surface), free    #correct tags now, correspond to left starting at z=0, and right is at z=L
fdim = mesh.topology.dim - 1
left_facets  = facet_tags.find(left_tag)
right_facets = facet_tags.find(right_tag)

ds = ufl.Measure("ds", domain=mesh, subdomain_data=facet_tags)

#-------------Define boundary condition------------------ 

u_D = np.array([0.0, 0.0, 0.0], dtype=default_scalar_type)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, left_facets), V) #clamped at the left #always one degree lower in fdim, since we apply bc on a surface (whole structure is volume)



#-------------------Define weak form for linear elasticity and assemble matrices--------------------------

def epsilon(u):
    return 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T) 

def sigma(u):
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(dim) + 2 * mu * epsilon(u)


a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx   # stiffness matrix
m = rho * inner(u, v) * dx  #mass matrix


K = assemble_matrix(form(a), bcs=[bc])
K.assemble()
M = assemble_matrix(form(m), bcs=[bc])
M.assemble()

#---------------Solve the eigenvalue problem----------------

eps = SLEPc.EPS().create(mesh.comm)
eps.setOperators(K,M)
eps.setProblemType(SLEPc.EPS.ProblemType.GHEP)

tol = 1e-6
eps.setTolerances(tol=tol, max_it=1000)
eps.setType(SLEPc.EPS.Type.KRYLOVSCHUR)

st = eps.getST()

# Set shift-and-invert transformation
st.setType(SLEPc.ST.Type.SINVERT)

eps.setDimensions(nev=200)  
eps.setFromOptions()
eps.solve()
eps.view()
eps.errorView()
nconv = eps.getConverged()

#--------------Filter out rigid body modes and print out the eigenfrequencies---------------------------------

modes = []
frequencies = []
f_tol = 1

for i in range(eps.getConverged()):
    val = eps.getEigenvalue(i)
    if val.real <=0:
        continue
    omega = np.sqrt(val.real)
    f = omega /(2*np.pi)
    if f < f_tol:
        continue
    print(f"Mode {i}: omega^2 = {val}, omega = {omega:.6f} rad/s, f = {f:.3f}  Hz")
    modes.append(i)
    frequencies.append(f)


#---------------------Plot the modes----------------------------------------------

chunk_size = 4

scale_factor = 5.0   # multiply displacements for visibility
show_scalar_bar = True
view_angle = "xz"  

for start in range(0, len(modes), chunk_size):
        chunk = modes[start:start + chunk_size]
        ncols = 2
        nrows = int(np.ceil(len(chunk) / ncols))
        p = pv.Plotter(shape=(nrows, ncols), window_size=(1600, 1200))

        for j, mode_index in enumerate(chunk):
            row = j // ncols
            col = j % ncols
            p.subplot(row, col)
        
            # Get the eigenvector from SLEPc
            eh = Function(V)  
            eps.getEigenvector(mode_index, eh.x.petsc_vec)
            
            dim = V.dofmap.index_map_bs  
            u = eh.x.array.reshape(-1, dim)

            # Build PyVista grid
            topology, cell_types, geometry = plot.vtk_mesh(mesh)
            grid = pv.UnstructuredGrid(topology, cell_types, geometry)
        
            grid.point_data["u"] = u 
            grid.point_data["|u|"] = np.linalg.norm(u, axis=1)
            
            arrows = grid.glyph(orient="u", scale="|u|", factor=0.2, tolerance=0.01)  

            # Plot
            p.add_mesh(arrows,scalars="|u|",cmap="viridis",show_scalar_bar=True)
            p.add_text(f"Mode {mode_index}", font_size=10)
            p.view_isometric()
            p.show_axes()
        
        p.show()

My question is: what could cause such a large number of near-zero eigenvalues in a clamped 3D elasticity eigenproblem where rigid-body modes should not exist? Could this indicate an issue with how Dirichlet boundary conditions are applied to K and/or M, or possibly a problem with the assembled matrices or the solver?

Thank you in advance!

2

As you haven’t given us the file

we can’t reproduce your example.

I would guess the issue is similar to:

where the solution is to remove the rows and columns that have the clamped conditions by using

K_sub, (idx_set_row, idx_set_row) = extract_sub_matrix(K, V, [bc])
M_sub, _ = extract_sub_matrix(M, V, [bc])

as described in the aforementioned post.
Please let me know if this helps.

All the best,
Jørgen

1 Like
3

Thank you for your quick reply! I was not able to insert the link to the github-repo in my original post, but the function file together with the STEP file are here: GitHub - NadineTab02/FEniCSx_Project

I will try your suggested method and see if it solves the problem.

/Nadine

4

No need to check the github-repo, your solution worked perfectly fine, thank you so much!

/Nadine

1 Like