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!

2

This formulation of the biharmonic problem will only be valid for C^1 elements, which are not in FEniCSx currently. I can recommend switching to a DG formulation or a mixed formulation.

One DG formulation is shown in this demo and I have some code here that implements it for clamped or simply supported boundary conditions. In addition, you might want to look into Dirichlet BCs affecting eigenvalues, as discussed here.

Since you are working on a rectangle, the mixed formulation (code here) should work fine and might be easier to implement, but on non-convex domains it will give you wrong results.

Hope this helps.

3

Hi Ottarhellan,

Thank you so much for your response. I will definitely look into these references.