Rayleigh-Ritz Method for Composite Plates

How the plate and sandwich panel tools solve a laminated plate: one Legendre basis for every edge mix, buckling and vibration as eigenvalue problems, deflection as a linear system, and a convergence check on every result.

The Rayleigh-Ritz method approximates the deflection of a composite plate by a finite sum of admissible functions, each of which already meets the geometric edge conditions, and picks the coefficients that make the plate's total potential energy stationary^{1, 2}. Buckling and free vibration then become symmetric eigenvalue problems, and deflection becomes a symmetric linear system.

Four tools run on it. Plate Buckling, Plate Deflection and Plate Vibration are documented on the plate reference, and the Sandwich Panel on the sandwich panel reference. Those pages cover the inputs and what each tool reports. This page covers the method underneath: the energy, the basis and how each edge condition enters it, the eigen solves, and how the number of terms is chosen and checked.

Why plates need an approximate method

Very few plates have a closed form. The classic one is a plate simply supported on all four edges with D₁₆ = D₂₆ = 0, where a double sine series satisfies the plate equation exactly^{1, 3}. Clamp one edge or free another and the sine series no longer meets the edge conditions. Use a laminate with ±θ plies at different distances from the mid-plane, a quasi-isotropic [0/±45/90]_s included, and the bending-twisting terms D₁₆ and D₂₆ no longer vanish, so the series fails again. How large those terms are against D₁₁ and D₂₂ is tabulated as the bending-twisting parameters on the engineering constants page.

The Ritz method handles both. The edges go into the basis functions, and the full bending stiffness D from the laminate's ABD matrix goes into the energy. D₁₆ and D₂₆ are then part of the solution, not an error to be estimated afterwards³.

The energy the method minimises

The total potential energy of the plate is the bending strain energy U, plus the potential V_N of the in-plane resultants, minus the work W of the transverse loads, all in N·mm^{1, 3}:

The plate occupies 0 ≤ x ≤ a, 0 ≤ y ≤ b, with a and b in mm. The deflection w (mm) is measured along the load, so a positive pressure gives a positive w. The bending strain energy uses the full D (N·mm):

The D₁₆ and D₂₆ terms pair a curvature with the twist w_,xy. They sit in the energy from the start, so the solve carries bending-twisting coupling at full strength.

The in-plane resultants N_x, N_y and N_xy (N/mm) act through the slope of the deflected surface¹:

They are tension positive. Compression makes V_N negative, and that negative energy is what destabilises the plate. The pressure q (MPa) does work through w:

A patch load adds one more term of the same kind: its force P (N) divided by the patch area on the plate, integrated over that area. Each problem then takes its own part of Π. Buckling uses U + V_N with no transverse load. Deflection uses U − W with no in-plane load. Free vibration adds the kinetic energy of the moving plate to U + V_N, where V_N comes from an optional in-plane prestress.

The Legendre basis

With the non-dimensional coordinates ξ = x/a and η = y/b, the deflection is a double sum of m terms along x and n along y:

Each one-dimensional function is a boundary factor times a Legendre polynomial:

P_i is the Legendre polynomial of degree i on [−1, 1]. The exponent p belongs to the edge at x = 0 and q to the edge at x = a, and each follows that edge's condition (see Edge conditions and boundary exponents). Y_j(η) is built the same way from the edges at y = 0 and y = b. The factor ξ^p makes every X_i vanish at ξ = 0 when p ≥ 1, and its slope as well when p = 2, whatever P_i does. That is all the method asks of an admissible function.

The polynomials come from the Bonnet recurrence, with t = 2ξ − 1, P₀ = 1 and P₁ = t:

Each step takes three terms at any degree, and the first two derivatives follow from the same recurrence. That keeps the basis cheap to evaluate at high order.

The first four basis functions along one edge pair of a thin plate, each scaled to a peak |X| of 1. Left, simply supported at both ends: every curve is zero at ξ = 0 and ξ = 1 and free to slope there. Right, clamped at ξ = 0 and free at ξ = 1: every curve leaves ξ = 0 with zero value and zero slope, and the free end is left unconstrained.

Why this basis. Bhat⁴ replaced the beam characteristic functions of the classical plate Ritz analysis with orthogonal polynomials, and showed that they converge well on plates with free, simply supported and clamped edges. Bhat built the set by Gram-Schmidt orthogonalisation of a starting beam polynomial. This basis takes the idea but starts from the boundary factor alone, and builds the higher members by the Legendre recurrence. Gram-Schmidt loses orthogonality in floating point at high order. The recurrence does not, and it keeps the stiffness matrix far better conditioned than the plain monomials ξ^p+i(1 − ξ)^q of the same degree would.

The choice does not change the answer. The Ritz solution depends only on the space the functions span, here the boundary factor times every polynomial of degree below m, so any basis of that space gives the same result in exact arithmetic. What the Legendre form buys is reaching that result in floating point. Beam characteristic functions are not used: each edge pair needs its own transcendental function, and their hyperbolic terms cancel catastrophically in floating point at high mode numbers.

Edge conditions and boundary exponents

The basis enforces the geometric conditions, zero deflection and zero slope or rotation, exactly. The natural conditions on moment and shear are never imposed: the minimum-energy solution approaches them as terms are added^{1, 2}. So each edge condition comes down to one exponent per field.

A thin plate clamped at x = 0, free at x = a and simply supported on both y edges. Its boundary factors are ξ² along x and η(1 − η) along y.

Thin plate: one field

The thin (Kirchhoff) plate has one field, w, and one exponent per edge. This table belongs to the thin-plate solve:

n is the edge normal, M_n the edge moment (N·mm/mm) and V_n the Kirchhoff edge shear (N/mm).

Mindlin plate: three fields

The Mindlin plate has three fields: w and the rotations ψ_x and ψ_y of the plate's normals (see Transverse shear: the Mindlin plate). On an edge, one rotation turns about the edge (the normal rotation) and the other along it (the tangential rotation). On an x = const edge ψ_x is the normal rotation and ψ_y the tangential one, and on a y = const edge the roles swap. Each field takes its own exponent¹:

The clamped edge takes exponent 1 on w here and 2 on the thin plate. The slope condition has not gone away. It has moved from w onto the normal rotation, which is now a field of its own. The simply supported edge is the hard one, holding the tangential rotation as well as w, which matches the thin plate's simply supported edge in the thin limit. The soft simply supported edge, with that rotation free, is not offered.

On both plates the natural conditions hold only in the limit, so zero moment on a free or simply supported edge is approached, not met, at any finite number of terms. Where two free edges meet, the corner condition is natural too, and it is the slowest to converge. Moments along a free edge converge more slowly than w.

Edge sets with no solution

With no clamped edge and at most one simply supported edge, the basis admits a rigid-body motion w = α + βx + γy. It has zero curvature and so zero strain energy, the stiffness matrix K is singular, and there is no solution to report². Four free edges admit all three motions. One simply supported edge with three free ones admits the rotation about that edge. The tools reject both before they assemble anything.

Every other set works, including a cantilever plate clamped on one edge (CFFF). Zero curvature everywhere forces w to be linear, and a linear function that vanishes with its normal slope along one edge, or vanishes along two different edges, is zero. The Mindlin plate follows the same rule: zero strain energy there needs zero curvature and zero shear strain, which leaves the same rigid-body motions.

Stiffness from one-dimensional integrals

Each basis function is a product X_i(ξ)Y_j(η), so every term of U splits into one integral over ξ times one over η¹:

with (r) the r-th derivative. Substituting the series into U gives U = ½c^TKc. The D₁₁ term of K and the D₁₆ coupling term show the pattern:

The D₁₂, D₂₂, D₂₆ and D₆₆ terms follow the same way, each a product of two one-dimensional tables, and K is symmetric. An X_i is a polynomial of degree at most m + 3, so Gauss-Legendre quadrature with m + 4 points, exact to degree 2m + 7, integrates every product exactly. The tables carry no quadrature error.

Buckling as an eigenvalue problem

The applied pattern N̂ = (N̂_x, N̂_y, N̂_xy) in N/mm is scaled by a load factor λ. Its potential V_N = ½λc^TK_Gc defines the geometric stiffness K_G, built from the same tables:

Stationarity of U + V_N gives the buckling eigenproblem^{1, 3}:

K is positive definite once the edge set has a solution, but K_G is in general indefinite. The problem is therefore reduced to a standard symmetric one through the Cholesky factor of K⁵:

Read the spectrum like this:

• The lowest positive λ is the buckling load: the plate buckles at λ times the applied pattern.
• A negative λ means the reversed pattern buckles the plate, at |λ| times the pattern.
• μ = 0 is a shape the load does not excite, such as a shape with no slope along x under pure N_x. It has no finite λ and is dropped.
• No positive λ means the pattern does not buckle the plate. Pure tension gives only negative λ.

Convergence from above. The bases for m and m + 2 terms are nested, so by the Rayleigh-Ritz minimax property each positive λ falls as terms are added¹. A truncated series makes the plate too stiff, and too few terms overestimate the buckling load. That is the unconservative direction. Higher modes converge more slowly than the first.

Hand check. For a simply supported plate with D₁₆ = D₂₆ = 0 under N_x = −N₀, the sine series is exact^{1, 3}:

Here m counts half-waves along x, not series terms. Take a [0/90]_s T300/5208 laminate of four 0.125 mm plies, with D₁₁ = 1670.6, D₁₂ = 30.2, D₂₂ = 331 and D₆₆ = 74.7 N·mm, on a 500 × 500 mm plate. One half-wave gives the minimum, and N₀ = (π²/500²)(1670.6 + 2(30.2 + 2 × 74.7) + 331) = 0.0932 N/mm. The Ritz series converges to that value from above.

Deflection as a linear system

With no in-plane load, stationarity of U − W gives a linear system for the coefficients¹:

The pressure term of the load vector f is q·ab times the integral of X_i over ξ times the integral of Y_j over η, and a patch load adds its own term. K is the same positive definite matrix as in buckling, so the system is solved with its Cholesky factor.

Energy bound. At the Ritz solution Π = −½f^Tc, and the Ritz minimum lies above the exact one. The load-weighted deflection f^Tc therefore rises toward its exact value as terms are added¹. The truncated plate is too stiff, and too few terms underestimate the load-weighted deflection.

The tools take a patch load rather than a point load, because under a concentrated force neither the thin plate's moments nor the Mindlin plate's deflection would converge. The plate reference explains why.

Free vibration and repeated frequencies

For free vibration the kinetic energy of the moving plate joins the strain energy. On the thin plate, with I₀ the mass per unit area (t/mm²) and I₂ the rotary inertia per unit area (t), it is the translation plus the rotation of the normals⁶:

The dot is the time derivative. Substituting the series gives a mass matrix M from the same one-dimensional tables, with I₀ on the products of the functions and I₂ on the products of their first derivatives. On the Mindlin plate the rotary inertia acts on ψ_x and ψ_y instead of on the slope of w. For a harmonic motion, stationarity gives the vibration eigenproblem^{2, 6}:

ω is in rad/s and f in Hz. K_G(N₀) is the geometric stiffness of an optional in-plane prestress N₀ (N/mm, tension positive), with no load factor. The modes are normalised in the mass, so a mode is a shape with no physical amplitude.

The reduction factors the stiffness, not the mass⁵:

M is singular on a Mindlin plate with I₂ = 0, because its rotation blocks vanish, while K + K_G is positive definite whenever the plate has not buckled. The reduction never inverts M, so that case solves as it stands.

A prestress at or beyond the buckling load takes away that positive definiteness, and the Cholesky factorisation meets a non-positive pivot. The tools report a plate buckled by its prestress, with no frequencies. It is a result, not an error.

Frequencies converge from above, like buckling loads. The nested bases make each ω fall as terms are added, and the higher modes converge more slowly than the fundamental.

Repeated frequencies

A square plate with the same condition on all four edges, or on opposite pairs that map onto each other when x and y swap, has pairs of equal frequencies by symmetry. Leissa's closed form for the simply supported isotropic plate without rotary inertia (I₂ = 0) shows it²:

Here m and n count half-waves along x and y, and ρ is the mass per unit area. On a square, modes (1,1), (1,2), (2,1), (2,2), (1,3) and (3,1) come out in the ratio 1 : 2.5 : 2.5 : 4 : 5 : 5, so modes 2 and 3 are one pair and modes 5 and 6 another. For an isotropic 200 × 200 mm plate, 1 mm thick, with E = 70 GPa, ν = 0.3 and a density of 2.7 g/cm³, D = Eh³/12(1 − ν²) = 6410.3 N·mm, ω₁₁ = 760.4 rad/s and f₁ = 121.0 Hz. With m = n terms the Ritz discretisation has the same symmetry, so its eigenvalues repeat to round-off.

Deflated Lanczos

A plain Lanczos run cannot resolve a repeated root. A single-vector Lanczos run builds its search space from one start vector, and that space meets a two-dimensional eigenspace in only one direction. In exact arithmetic it finds one copy of the pair and never the other⁵. Round-off may bring the second copy in late, but not reliably.

The solver therefore runs Lanczos with full reorthogonalisation and deflates each run against the eigenpairs already found: every new vector is made orthogonal to them as well as to the earlier vectors. Once the first run has the six lowest frequencies, a confirming restart begins from a fresh vector, deflated against all six. A value below the current sixth is a missed copy of a repeated root. It is locked in, the old sixth drops out, and another restart confirms. A value at or above the sixth ends the search, so a plate with no repeated root in its first six modes costs two runs. Buckling asks for one eigenpair and skips the restart, because a missed copy cannot change the lowest value.

Transverse shear: the Mindlin plate

Thin-plate theory assumes the normals to the mid-plane stay straight and normal. The Mindlin plate, also called first-order shear deformation theory, keeps them straight but lets them rotate on their own^{1, 6}. That adds two fields, the rotations ψ_x and ψ_y (rad), to the in-plane displacements u and v (mm):

The curvatures now come from the rotations, and each transverse shear strain is the difference between a rotation and the slope of w:

The strain energy gains a shear part U_s, carried by the transverse shear stiffness A_s (N/mm): the ply shear stiffnesses summed over the ply thicknesses t_k (mm) and multiplied by the shear correction factor 5/6¹:

Each field uses the Legendre form of the basis with its own exponents, from the Mindlin exponent table. K_G and the transverse loads act on the w block only: the in-plane resultants act through the slope of w, and the loads do work on w. The Mindlin plate has three times the unknowns per term.

In the thin limit, as A_s grows against D/a², the shear energy forces γ → 0, so ψ_x = w_,x and ψ_y = w_,y, and the Mindlin plate recovers the thin plate. A low-order displacement element can lock there, its shear energy stiffening the plate far beyond the thin solution. The polynomial basis spans ψ = ∇w and does not lock.

The tools solve every plate both ways. The difference between the two answers is the shear effect they report, explained under transverse shear and the thin plate. A sandwich panel enters the same solve through the shear stiffness of its core, covered on the sandwich panel reference.

Series terms and the convergence check

No tool asks for a number of terms. Each sets m and n from the plate's aspect ratio, by theory and load:

The longer side gets more terms because buckles and modes have more half-waves along it. The Mindlin plate has three fields per term, so its counts are set on their own.

The Mindlin result is then solved again at (m + 2) × (n + 2), and flagged as not converged when the two solves differ by more than 1 %. The thin-plate comparison is not checked again. What gets compared depends on the problem:

• Buckling: the lowest positive λ. If only one of the two solves finds that the load buckles the plate, the result is not converged.
• Deflection: |w_max|, the moment field and q_max, the largest transverse shear resultant, each with its own flag.
• Vibration: each of the six frequencies on its own, mode i of the check against mode i of the result.

The numbers on screen come from the m × n solve. Because the bases are nested, the buckling load and the frequencies can only fall, and the load-weighted deflection can only rise, between the two solves. The check measures that last increment, not the remaining error, so a slowly converging case can pass while still more than 1 % from its limit. Two modes that swap order between the solves compare different shapes and can raise a flag that is not a convergence failure. How the chips appear on each result is under series terms and the convergence check on the plate reference.

Plate Vibration on a simply supported quasi-isotropic plate. Each of the six modes carries its own convergence chip, the bending + shear frequencies sit beside the bending only ones, and the Series terms row gives m × n for the Mindlin solve.

Numerical solution

This solver runs in double precision (float64), unlike the rest of the calculation layer, which works in arbitrary precision. There is no generalized symmetric eigensolver in that arithmetic, and an arbitrary-precision solve of a few hundred unknowns is too slow to repeat on every input change and again for the convergence check. The Legendre basis keeps the round-off down.

K is factorised by Cholesky. The lowest eigenpairs come from the Lanczos method with full reorthogonalisation⁵, with the reduced operator S or R applied as a back substitution, a matrix product and a forward substitution, never formed as a matrix. The start vector is fixed and seeded, so the same input gives the same modes. A run stops when the residual of each wanted Ritz value, the eigenvalue estimate Lanczos produces, is below 10⁻¹⁰ of the largest |Ritz value|.

Assumptions and limits

• Bifurcation of a perfect plate. The buckling load carries no imperfection knock-down, and post-buckling behaviour is not modelled.
• Uniform pre-buckling resultants. The in-plane resultants are taken as given and uniform over the plate. How the load is introduced and how the edges are restrained in plane are not modelled.
• Hard simply supported edges only. A simply supported edge holds the tangential rotation as well as w.
• First-order shear deformation. The shear strain is constant through the thickness, with a single correction factor of 5/6¹. On a very thick or strongly mixed stack that factor is an approximation.
• No in-plane and transverse coupling in deflection. The deflection solve ignores the stiffening of in-plane tension and the softening of compression. The vibration solve does couple its prestress in, through K + K_G(N₀).
• Vibration inertia. I₂ is taken about the mid-plane and the coupling inertia I₁ is dropped. I₁ vanishes when the density is symmetric about the mid-plane, so the omission matters only for an unsymmetric laminate. No damping and no forced response.
• Natural conditions in the limit. Zero moment and shear on free edges, and the free-corner condition, hold only as the number of terms grows.

Author: Rick Schrijver

References

1. Reddy, J.N. (2004). Mechanics of Laminated Composite Plates and Shells: Theory and Analysis, 2nd ed. CRC Press. ISBN: 978-0849315923. (Classical laminated plate theory and its energy functional, Rayleigh-Ritz solutions for bending, buckling and vibration of rectangular plates, first-order shear deformation with the 5/6 shear correction) doi:10.1201/b12409
2. Leissa, A.W. (1969). Vibration of Plates. NASA SP-160. Washington, DC: NASA. (Rayleigh-Ritz method for rectangular plates with mixed edges, natural boundary conditions; Eq. 4.20, p. 44 for the simply supported isotropic frequencies) apps.dtic.mil
3. Whitney, J.M. (1987). Structural Analysis of Laminated Anisotropic Plates. Technomic. ISBN: 978-0877625186. (Energy of anisotropic plates with D₁₆ and D₂₆, buckling of laminated plates)
4. Bhat, R.B. (1985). Natural frequencies of rectangular plates using characteristic orthogonal polynomials in Rayleigh-Ritz method. Journal of Sound and Vibration, 102(4), 493-499. (Orthogonal polynomial bases for plates with mixed edges) doi:10.1016/S0022-460X(85)80109-7
5. Golub, G.H. & Van Loan, C.F. (2013). Matrix Computations, 4th ed. Johns Hopkins University Press. ISBN: 978-1421407944. (Cholesky factorisation, the symmetric-definite generalized eigenproblem, and the Lanczos method with full reorthogonalisation and its behaviour on repeated eigenvalues)
6. Mindlin, R.D. (1951). Influence of rotatory inertia and shear on flexural motions of isotropic, elastic plates. Journal of Applied Mechanics, 18, 31-38. (The plate theory with transverse shear deformation and rotary inertia)

Frequently Asked Questions