Kirchhoff plate (biharmonic) eigenproblem in FEniCSx gives near-zero frequencies

General
1

Hi everyone,

I am trying to compute natural frequencies of a simply supported rectangular plate using a Kirchhoff plate model in FEniCSx with SLEPc.

For my parameters:

a = 1.0      # plate length (m)
b = 0.5      # plate width (m)
t = 0.002    # plate thickness (m)
E = 2.1e11   # Young's modulus (Pa)
nu = 0.3     # Poisson's ratio (-)
rho = 7850   # density (kg/m^3)
plate_nx= 100 # Elements in x
plate_ny= 50 # Elements in y

My analytical result is for the first frequency ≈ 24.59 Hz. However, SLEPc returns first frequency ≈ 0.002 Hz.

domain = mesh.create_rectangle(
        comm, [[0.0, 0.0], [a, b]], [plate_nx, plate_ny], cell_type=mesh.CellType.triangle)
    V = fem.functionspace(domain, ("Lagrange", 2))

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

    D = young * t**3 / (12.0 * (1.0 - nu**2))

    # Correct Kirchhoff plate weak form (biharmonic)
    a_form = D * ufl.inner(ufl.grad(ufl.grad(u)), ufl.grad(ufl.grad(v))) * ufl.dx
    m_form = rho * t * ufl.inner(u, v) * ufl.dx

    # Simply supported boundary conditions
    def simply_supported_boundary(x):
   
        tol = 1e-8
        return (np.isclose(x[0], 0.0, atol=tol)
                | np.isclose(x[0], a, atol=tol)
                | np.isclose(x[1], 0.0, atol=tol)
                | np.isclose(x[1], b, atol=tol))

    dofs = fem.locate_dofs_geometrical(V, simply_supported_boundary)
    bc = fem.dirichletbc(PETSc.ScalarType(0.0), dofs, V)
  # Assemble matrices with standard Dirichlet handling
    K = assemble_matrix(fem.form(a_form), bcs=[bc])
    K.assemble()
    # K.shift(1e-8)  # Small value for numerical stability
    M = assemble_matrix(fem.form(m_form), bcs=[bc])
    M.assemble()
# Setup eigensolver
    solver = SLEPc.EPS().create(domain.comm)
    solver.setProblemType(SLEPc.EPS.ProblemType.GHEP)
    
    solver.setType(SLEPc.EPS.Type.KRYLOVSCHUR)
    solver.setDimensions(num_modes, max(4 * num_modes, 40))
    solver.setOperators(K, M)
    solver.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_REAL)
    solver.setTarget(shift)
    st = solver.getST()
    st.setType(SLEPc.ST.Type.SINVERT)
    st.setShift(shift)
    solver.setTolerances(1e-8, 500)
    solver.setFromOptions()
    solver.solve()

    nconv = solver.getConverged()
    reason = solver.getConvergedReason()
  1. Is this formulation valid in FEniCSx?
  2. Are the near-zero frequencies due to using (C^0) Lagrange elements for a fourth-order problem?

Any guidance or references would be appreciated.

Thanks!