Free vibration of a simply supported plate

Natural frequencies and mode shapes from a 3D FEniCSx model, checked against Navier’s plate solution

FEniCSx

modal analysis

SLEPc

validation

Modal analysis of a thin rectangular plate: a 3D solid model in FEniCSx with a SLEPc eigensolver, validated against Navier’s Kirchhoff-plate solution.

Published

October 10, 2026

A thin rectangular plate, simply supported on all four edges and left to vibrate freely. Its lowest natural frequencies and mode shapes come from a 3D solid model in FEniCSx, solved as a generalised eigenvalue problem with SLEPc, and are checked against Navier’s classical solution for a Kirchhoff plate.

Figure 1: Mode shapes of the mid-plane deflection \(w\), normalised to a unit peak: the colour shows the shape, and the surface oscillates as \(w \cos\omega t\), exaggerated for display. Pick a mode, drag to rotate, hover to read values (the animation holds still).

Modes extracted 10 lowest, by shift-and-invert

Largest frequency error 0.36% FE vs Navier, all modes

Lowest MAC 0.999993 mode-shape correlation

The model

Plate	600 × 400 × 5 mm, simply supported on all four edges
Material	linear elastic, \(E\) = 210 GPa, \(\nu\) = 0.3, \(\rho\) = 7,850 kg/m³
Mesh	structured hexahedra, 30 × 20 in-plane and 2 layers through the thickness, with quadratic (27-node) Lagrange elements
Eigensolver	SLEPc Krylov–Schur, shift-and-invert about zero with a MUMPS LU factorisation, 10 eigenpairs

Analytical solution

In Kirchhoff plate theory, the transverse deflection \(w(x, y, t)\) of a plate of thickness \(h\) obeys

\[ D\,\nabla^4 w + \rho h\,\frac{\partial^2 w}{\partial t^2} = 0, \qquad D = \frac{E h^3}{12\,(1 - \nu^2)}. \tag{1}\]

A vibration in a single mode, \(w = W(x, y)\cos\omega t\), separates space from time and leaves an eigenvalue problem for the mode shape \(W\) and its angular frequency \(\omega\):

\[ D\,\nabla^4 W = \rho h\,\omega^2\, W. \tag{2}\]

On a simply supported \(a \times b\) plate, every edge has zero deflection and zero bending moment: \(W = \partial^2 W/\partial x^2 = 0\) on \(x = 0, a\) and \(W = \partial^2 W/\partial y^2 = 0\) on \(y = 0, b\). Navier’s double-sine solution meets all of them. Its mode shapes are

\[ W_{mn}(x, y) = \sin\frac{m\pi x}{a}\,\sin\frac{n\pi y}{b}, \qquad m, n = 1, 2, \dots, \tag{3}\]

with \(m\) half-waves along \(x\) and \(n\) along \(y\), separated by \(m - 1\) and \(n - 1\) straight nodal lines (white in Figure 1). Substituted into Equation 2, they give the natural frequencies

\[ \omega_{mn} = \pi^2\left(\frac{m^2}{a^2} + \frac{n^2}{b^2}\right)\sqrt{\frac{D}{\rho h}}, \qquad f_{mn} = \frac{\omega_{mn}}{2\pi}. \tag{4}\]

The theory neglects transverse shear and rotary inertia, which lower the frequencies of a real plate by \(O\big((h/\lambda)^2\big)\) for a wavelength \(\lambda\).

FEniCSx implementation

The plate is meshed as a 3D elastic solid rather than a plate. Seeking free vibrations \(\underline{u}(\underline{x}, t) = \underline{\phi}(\underline{x})\,e^{i\omega t}\) turns the weak form of elastodynamics into a generalised eigenvalue problem:

\[ \int_\Omega \underline{\underline{\sigma}}(\underline{\phi}) : \underline{\underline{\varepsilon}}(\underline{v})\,\mathrm d\Omega = \omega^2 \int_\Omega \rho\, \underline{\phi} \cdot \underline{v} \,\mathrm d\Omega \quad \forall\, \underline{v} \qquad \Longrightarrow \qquad \underline{\underline{K}} \underline{\Phi} = \omega^2 \underline{\underline{M}} \underline{\Phi}. \tag{5}\]

Condensed excerpts from the scripts (... marks omitted lines):

3D/simulation.py

1displacement_space = fem.functionspace(domain, ("P", element_degree, (geom_dim,)))
...
lame_lambda = young_modulus*poisson_ratio/((1+poisson_ratio)*(1-2*poisson_ratio))
shear_modulus = young_modulus/(2*(1+poisson_ratio))

def strain(displacement):
    return ufl.sym(ufl.grad(displacement))

def stress(strain_tensor):
    return lame_lambda*ufl.tr(strain_tensor)*ufl.Identity(geom_dim) + 2*shear_modulus*strain_tensor

# Free vibration: K u = omega^2 M u
2stiffness_form = ufl.inner(stress(strain(trial)), strain(test))*ufl.dx
3mass_form = density*ufl.inner(trial, test)*ufl.dx

1: Quadratic vector Lagrange space (element_degree = 2): with 2 element layers it resolves bending through the thickness.
2: The stiffness form is the one of static linear elasticity.
3: Consistent mass matrix, weighted by the density \(\rho\).

3D/simulation.py

# Simply supported: the edge faces are held transversally but rotate freely. The
# mid-plane edge lines are held normal to the edge, which suppresses the in-plane
# rigid body motion without restraining the bending rotation.
transverse_space = displacement_space.sub(2)
supported_edge_dofs = fem.locate_dofs_topological(
    transverse_space, bc_dim, tags.find(supported_edge_tag)
1)

def mid_plane_edge_condition(normal_component, edge_coordinates):
    def on_mid_plane_edge(x):
        on_edge = np.zeros(x.shape[1], dtype=bool)
        for edge_coordinate in edge_coordinates:
            on_edge |= np.isclose(x[normal_component], edge_coordinate, atol=geometry_tolerance)
        return on_edge & np.isclose(x[2], 0.0, atol=geometry_tolerance)

    component_space = displacement_space.sub(normal_component)
2    mid_plane_edges = mesh.locate_entities_boundary(domain, edge_dim, on_mid_plane_edge)
    mid_plane_edge_dofs = fem.locate_dofs_topological(component_space, edge_dim, mid_plane_edges)
    return fem.dirichletbc(np.float64(0.0), mid_plane_edge_dofs, component_space)

boundary_conditions = [
    fem.dirichletbc(np.float64(0.0), supported_edge_dofs, transverse_space),
    mid_plane_edge_condition(0, (0.0, plate_length)),
    mid_plane_edge_condition(1, (0.0, plate_width)),
]

1: \(u_z = 0\) on the four edge faces: the transverse support of the plate.
2: Mesh edges (1D entities) on the mid-plane perimeter. Fixing only the in-plane normal displacement there leaves the bending rotation free, as a simple support should.

3D/simulation.py

1stiffness_matrix = assemble_matrix(fem.form(stiffness_form), bcs=boundary_conditions, diag=1)
stiffness_matrix.assemble()
mass_matrix = assemble_matrix(fem.form(mass_form), bcs=boundary_conditions, diag=0)
mass_matrix.assemble()

eigenvalue_solver = solve_generalized_eigenproblem(
    stiffness_matrix,
    mass_matrix,
    eigenvalue_count=vibration_mode_count,
2    problem_type=SLEPc.EPS.ProblemType.GHEP,
)
angular_frequencies_squared, vibration_modes = extract_eigenpairs(
    eigenvalue_solver, displacement_space
)

1: Constrained rows get a unit diagonal in \(\underline{\underline{K}}\) and a zero one in \(\underline{\underline{M}}\): they map to infinite eigenvalues, far from the shift, and never pollute the spectrum.
2: Generalised Hermitian problem: symmetric \(\underline{\underline{K}}\), positive semi-definite \(\underline{\underline{M}}\).

tools/eigensolver.py

eigenvalue_solver.setType(SLEPc.EPS.Type.KRYLOVSCHUR)
...
eigenvalue_solver.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_MAGNITUDE)
eigenvalue_solver.setTarget(shift)
spectral_transform = eigenvalue_solver.getST()
1spectral_transform.setType(SLEPc.ST.Type.SINVERT)
spectral_transform.setShift(shift)
linear_solver = spectral_transform.getKSP()
linear_solver.setType("preonly")
linear_solver.getPC().setType("lu")
2linear_solver.getPC().setFactorSolverType("mumps")

1: Shift-and-invert about \(\sigma = 0\): the lowest frequencies become the largest eigenvalues of \((\underline{\underline{K}} - \sigma\underline{\underline{M}})^{-1}\underline{\underline{M}}\), which Krylov–Schur finds first.
2: \(\underline{\underline{K}} - \sigma\underline{\underline{M}}\) is factorised once with MUMPS and reused at every iteration.

3D/post.py

1node_coordinates = displacement_space.tabulate_dof_coordinates()
...
mid_plane_nodes = np.isclose(node_coordinates[:, 2], 0.0, atol=1e-6*plate_thickness)
mid_plane_x = node_coordinates[mid_plane_nodes, 0]
mid_plane_y = node_coordinates[mid_plane_nodes, 1]

2def modal_assurance_criterion(first_shape, second_shape):
    return np.dot(first_shape, second_shape)**2/(
        np.dot(first_shape, first_shape)*np.dot(second_shape, second_shape)
    )

candidate_half_wave_counts = [
    (half_wave_count_x, half_wave_count_y)
    for half_wave_count_x in range(1, vibration_mode_count + 1)
    for half_wave_count_y in range(1, vibration_mode_count + 1)
]
...
for mode_index, (fe_natural_frequency, vibration_mode) in enumerate(
    zip(fe_natural_frequencies, vibration_modes)
):
3    fe_mid_plane_deflection = vibration_mode[mid_plane_nodes, 2]
    mac_per_candidate = [
        modal_assurance_criterion(
            fe_mid_plane_deflection,
            mode_shape(half_wave_count_x, half_wave_count_y, mid_plane_x, mid_plane_y, plate_length, plate_width),
        )
        for half_wave_count_x, half_wave_count_y in candidate_half_wave_counts
    ]
4    best_candidate = int(np.argmax(mac_per_candidate))
    half_wave_count_x, half_wave_count_y = candidate_half_wave_counts[best_candidate]

1: Coordinates of every node of the quadratic space. With an even number of element layers, the mid-plane \(z = 0\) is a layer of nodes: a 61 × 41 grid.
2: The MAC of Equation 6.
3: \(\underline{\phi}{}_i\): the transverse component \(u_z\) of the eigenvector at the mid-plane nodes.
4: The candidate with the largest MAC labels the FE mode with its \((m, n)\).

Results

The modal assurance criterion (MAC) was used to compare the form of the analytical and FE modes. The MAC compares two mode shapes \(\underline{\phi}\) and \(\underline{\psi}\) with the formula :

\[ \text{MAC}(\underline{\phi}, \underline{\psi}) = \frac{\big(\underline{\phi}\cdot\underline{\psi}\big)^2} {\big(\underline{\phi}\cdot\underline{\phi}\big)\big(\underline{\psi}\cdot\underline{\psi}\big)} \tag{6}\]

It is the squared cosine of the angle between the two vectors, so the sign and scale of either do not matter: it is 1 for proportional shapes and 0 for orthogonal ones, as any two different Navier modes are. For each FE mode \(i\):

Sample both shapes at the same points. \(\underline{\phi}{}_i\) holds the transverse displacement \(u_z\) of FE mode \(i\) at the 2,501 nodes of the mid-plane \(z = 0\), and \(\underline{\psi}{}_{mn}\) holds \(W_{mn}\) (Equation 3) at the same nodes.
Compare with every candidate. Compute \(\text{MAC}(\underline{\phi}{}_i, \underline{\psi}{}_{mn})\) for all \(m, n \le\) 10.
Keep the best match. FE mode \(i\) takes the \((m, n)\) with the largest MAC, and its frequency is compared with \(f_{mn}\) (Equation 4).

The code is in the Mode pairing tab above.

Table 1: Natural frequencies: 3D finite elements against Navier’s solution (Equation 4), with the MAC of each pair (Equation 6).

Mode	(m, n)	Navier [Hz]	FE [Hz]	Δf	MAC
1	(1, 1)	110.98	110.72	−0.23%	1.000000
2	(2, 1)	213.42	212.85	−0.27%	0.999999
3	(1, 2)	341.47	341.36	−0.03%	0.999999
4	(3, 1)	384.15	383.56	−0.16%	0.999998
5	(2, 2)	443.91	442.76	−0.26%	0.999998
6	(3, 2)	614.64	612.42	−0.36%	0.999996
7	(4, 1)	623.18	623.03	−0.02%	0.999998
8	(1, 3)	725.62	726.75	+0.15%	0.999998
9	(2, 3)	828.06	827.42	−0.08%	0.999995
10	(4, 2)	853.67	850.74	−0.34%	0.999993

Figure 2: MAC (Equation 6) of every FE mode against the lowest Navier modes, both in order of frequency.

Every mode is identified. Each FE mode matches exactly one Navier mode: the diagonal MAC is at least 0.999993 and every off-diagonal term rounds to zero at four decimals. The FE spectrum has Navier’s mode order, even for the closest pair, (3, 2) and (4, 1), only 10.6 Hz apart.
Frequencies agree within 0.36%, from −0.36% to +0.15% of Navier’s. Two effects of similar size compete: the transverse shear and rotary inertia that Kirchhoff theory neglects lower the frequencies of the 3D model, while the discretisation raises them, since conforming elements with a consistent mass matrix overestimate every natural frequency, the more so for short wavelengths.