Structural Modal Analysis

General
, , , ,
2

one should remove the strongly enforced dofs from the system, i.e. something along the lines of:

I would do something along the lines of:

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, gmsh as gmshio
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
from ufl import sym, Identity, div
import os

import dolfinx
def extract_sub_matrix(A:PETSc.Mat, V: dolfinx.fem.FunctionSpace, bcs:list[dolfinx.fem.DirichletBC])->PETSc.Mat:
    num_dofs_local = V.dofmap.index_map.size_local*V.dofmap.index_map_bs
    num_ghosts = V.dofmap.index_map.num_ghosts*V.dofmap.index_map_bs

    not_dirichlet = np.full(num_dofs_local+num_ghosts, 1, dtype=np.int32)

    for bc in bcs:
        not_dirichlet[bc.dof_indices()[0]] = 0

    idx_set = np.flatnonzero(not_dirichlet).astype(np.int32)
    idx_set_row = idx_set[idx_set<num_dofs_local]
    idx_set_col = idx_set[idx_set<num_dofs_local]
    isx_row = PETSc.IS(A.getComm()).createGeneral(idx_set_row)
    isx_col = PETSc.IS(A.getComm()).createGeneral(idx_set_col)
    lgm_row = A.getLGMap()[0].applyIS(isx_row)
    lgm_col = A.getLGMap()[1].applyIS(isx_col)

    A_sub = A.createSubMatrix(lgm_row, lgm_col)
    A_sub.assemble()
    return A_sub, (idx_set_row, idx_set_col)


PART = "STEEL_BEAM"  # Name of the part being analyzed (used for file naming)

GEO_FILE = PART + ".geo"
MSH_FILE = 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
try:
    from dolfinx.io import gmshio
    mesh, cell_tags, facet_tags = gmshio.read_from_msh(PART + ".msh", MPI.COMM_WORLD, 0, gdim=3)
except ImportError:
    from dolfinx.io import  gmsh as gmshio
    mesh_data = gmshio.read_from_msh(PART + ".msh", MPI.COMM_WORLD, 0, gdim=3)
    mesh = mesh_data.mesh
    cell_tags = mesh_data.cell_tags
    facet_tags = mesh_data.facet_tags
# 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()
K_sub, (idx_set_row, idx_set_row) = extract_sub_matrix(K, V, [bc])
M_sub, _ = extract_sub_matrix(M, V, [bc])

# 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_sub, M_sub)

# solve the problem
solver.solve()

xr, xi = K_sub.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 = []

d = Function(V)
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())
        d.x.array[idx_set_row] = vect

        with dolfinx.io.VTXWriter(d.function_space.mesh.comm, f"EigenVector{i}.bp", [d]) as bp:
            bp.write(0.0)

which would yield:

Mode 0: 4.137863516946098 Hz
Mode 1: 25.930998859980065 Hz
Mode 2: 41.14547192302046 Hz
Mode 3: 72.63890229130486 Hz
Mode 4: 142.46413512919722 Hz
Mode 5: 152.11218623471086 Hz
Mode 6: 235.7921718294653 Hz
Mode 7: 254.833591000687 Hz
Mode 8: 352.78470736837915 Hz

I’m not an expert in modal analysis, maybe @jackhale or @michalhabera can comment on anything obviously wrong?

1 Like