Frequency-Domain Isogeometric Galerkin Analysis of Generalised
Kirchhoff–Love Plates
Abstract
The Kirchhoff–Love theory models the bending and vibration of thin plates by assuming that
straight lines normal to the mid-surface remain straight and normal after deformation and that
the plate thickness is constant【928766934547163†L140-L153】. For an isotropic, homogeneous
plate of constant thickness (h), the transverse displacement (w(x,y,t)) satisfies a fourth-order
partial differential equation (PDE) where the bending rigidity, tension, and damping terms
balance inertial effects. This study develops and analyzes a MATLAB implementation that
combines isogeometric analysis (IGA) with a Fourier–frequency formulation to simulate a
generalised Kirchhoff–Love model with both viscous and visco-elastic damping. The
time-domain PDE is transformed to the frequency domain, leading to a family of complex
elliptic problems. A weak form is derived and discretised using (C^1)-continuous NURBS basis
functions, and simply supported, clamped and free boundary conditions are imposed in terms
of the displacement, slope, bending moment and shear force. The inverse Fourier transform is
numerically approximated using Gauss–Hermite quadrature【200742885264159†L117-L132】
to reconstruct the time response. A manufactured solution validates the method and illustrates
its convergence properties. The proposed framework offers an accurate and flexible approach
for dynamic analysis of thin plates and can be extended to more complex damping and
boundary conditions.
1 Introduction
Thin plates are central to many engineering structures, including aircraft wings, vehicle panels
and building floors. When the thickness is small compared with the other dimensions, the
Kirchhoff–Love plate theory provides an efficient model by assuming that normal lines to
the mid-surface remain straight and normal during deformation and that the plate thickness
does not change【928766934547163†L140-L153】. These kinematic assumptions reduce the
three-dimensional elasticity problem to a two-dimensional description in terms of a single
transverse displacement field (w(x,y,t)). For homogeneous plates the resulting governing
equation is a biharmonic operator acting on (w) coupled to inertial terms
【928766934547163†L987-L999】. The accuracy of solutions depends on correctly representing
both the geometry and the required (C^1) continuity across element boundaries. Isogeometric
analysis (IGA) meets these requirements by using non-uniform rational B-splines (NURBS) to
represent both the geometry and the field variables. NURBS basis functions of degree (p) with
knot multiplicity (m) are (C^{p-m}) continuous across knot lines and (C^) between distinct
knots【480751958331550†L170-L180】. This high continuity enables direct discretisation of
fourth-order PDEs without introducing rotational degrees of freedom. Moreover, representing
the geometry exactly reduces approximation and dispersion errors【480751958331550†L122-
L155】.
Dynamic plate problems present additional challenges because time derivatives couple to the
spatially stiff fourth-order operator. A common strategy uses a Fourier transform in time,
converting the problem into a sequence of elliptic problems in the frequency domain. Once the
frequency response is computed for discrete frequencies, the time response is recovered by an
inverse transform. Gauss–Hermite quadrature approximates integrals of the form (_{-}{}e{-
x^2}f(x),dx) using weighted evaluations at the roots of Hermite polynomials
【200742885264159†L117-L132】. This technique is particularly efficient for inverse Fourier
transforms of decaying functions because the exponential weight damps high-frequency
contributions.
This paper develops an isogeometric Galerkin solver for a generalised Kirchhoff–Love plate
model with constant thickness and isotropic properties. The governing PDE includes terms
accounting for linear stiffness, membrane tension, visco-elastic damping and viscous damping.
Section 2 formulates the problem in the time domain and describes physically meaningful
boundary conditions for clamped, simply supported and free edges. Section 3 applies the
Fourier transform to obtain a frequency-domain PDE. Section 4 derives the weak form.
Section 5 describes the IGA discretisation and enforcement of boundary conditions. Section 6
explains the numerical inverse transform via Hermite quadrature. Section 7 presents a
manufactured solution and convergence study. Section 8 discusses results and Section 9
concludes.
2 Time-Domain Formulation
2.1 Governing equation
Let (^2) denote the mid-surface of an isotropic, homogeneous plate with constant thickness (h).
The transverse displacement (w(x,y,t)) satisfies the generalised Kirchhoff–Love model [ h, =
-K_0 w + T,^2 w - D,^4 w - K_1, + T_1,^2 + F(x,y,t), ] where () is the mass density, (D=E h^3/[12
(1-^2)]) is the constant flexural rigidity ((E) and () are Young’s modulus and Poisson’s ratio),
(T) is an in-plane tension coefficient, (K_0) is a linear stiffness, (K_1) represents viscous
damping, and (T_1) represents visco-elastic damping which tends to smooth high-frequency
oscillations in space. The body force (F(x,y,t)) acts on (). The term with coefficient (T_1) is often
used to model visco-elastic damping in vascular structures【755588771263318†L348-L360】.
Equation (1) is a fourth-order PDE in space and second order in time. The biharmonic operator
(^4 w = ^4 w/x^4 + 2 ^4 w/x2y2 + ^4 w/y^4) appears because the bending moment in
Kirchhoff–Love theory is proportional to the plate’s curvature vector
((w)=[w_{,xx},w_{,yy},2w_{,xy}]^T). The Laplacian (^2 w) arises from membrane tension and
visco-elastic damping.
2.2 Boundary conditions
On the boundary (), physically meaningful conditions relate the displacement and its
derivatives to the reactions imposed by supports. For isotropic plates the normal (n) and
tangential (t) directions on an edge determine the bending moment and shear force per unit
length. Let (M_n) denote the bending moment per unit length normal to the boundary and
(Q_n) the transverse shear force. For a plate of constant flexural rigidity, the bending moment
and shear force are given by [ M_n = -D ( + ), Q_n = -D ( + ()). ] Boundary conditions can
prescribe displacements, slopes, bending moments or shear forces. The three standard
boundary types are:
• Clamped edge. The edge is rigidly fixed, so the displacement and slope vanish: [ .w|_{c}
= 0, .|{_c} = 0. ] These constraints remove both rotation and translation at the
boundary【635095030438106†L76-L78】.
• Simply supported (supported) edge. The edge is free to rotate but cannot translate.
The displacement vanishes and the bending moment is zero: [ .w|_{s} = 0, .M_n|{_s} = 0.
] Equivalently, (^2 w/n^2 + ^2 w/t^2 = 0) along the edge. A simply supported edge thus
permits rotation but no vertical deflection【603759730802822†L852-L857】.
• Free edge. The edge has no support, so both the bending moment and the shear force
vanish: [ .M_n|_{f} = 0, .Q_n|{_f} = 0. ] Physically, this means that neither bending
moment nor shear force can be transmitted across the boundary. For example, the shear
force condition (Q_n=0) corresponds to the third derivative of the deflection vanishing
【663235947978228†L316-L319】. These boundary conditions are used to model free
edges in plate and shell problems【663235947978228†L316-L319】.
Initial conditions can be imposed as [w(x,y,0) = w_0(x,y), (x,y,0) = v_0(x,y).]
The present work treats simply supported boundaries in the numerical implementation but the
formulation accommodates clamped and free edges through appropriate enforcement of
conditions (3a)–(3c).
3 Transformation to the Frequency Domain
To decouple the time and space discretisations, we apply a Fourier transform in time to (1).
Define [ (x,y,) = {-}^{} w(x,y,t),e^{-it},dt, w(x,y,t) = {-}^{} (x,y,),e^{it},d. ] The transform
properties ({} = i) and ({} = -^2 ) yield the frequency-domain equation [ (x,y,) = (x,y,). ] For
each frequency (), (4) is a complex-valued elliptic PDE in space. The term (D,^4) corresponds to
bending rigidity, (T,^2) to membrane tension, (K_0) to linear stiffness, and (-h,^2) to inertia.
The visco-elastic term (iT_1,^2) and viscous term (iK_1) provide complex stiffness that
introduces damping. Boundary conditions (3a)–(3c) are applied at each frequency, with the
time derivatives replaced by multipliers of (i). For example, the bending moment condition on a
simply supported edge becomes (M_n = -D(^2 /n^2 + ^2 /t^2) = 0).
4 Weak Formulation
Let (V) denote the subspace of (H^2()) that satisfies the essential (Dirichlet) boundary
conditions. For clamped and simply supported edges, the displacement (w) is prescribed; for
free edges, the bending moment and shear force are prescribed (set to zero). Multiplying (4) by
a test function (vV) and integrating over () yields [
] where ((w)=[w_{,xx},w_{,yy},2 w_{,xy}]^T) and the symmetry of (D,) ensures that boundary
terms vanish on simply supported and free edges due to conditions (3b)–(3c). When applying
the weak form to clamped edges, integration by parts yields additional boundary terms
involving (w/n), which are set to zero by (3a). This variational formulation is suitable for
Galerkin discretisation.
5 IGA Discretisation and Boundary Enforcement
5.1 NURBS geometry and basis
IGA employs NURBS basis functions for both geometry and displacement. A univariate B-spline
basis of degree (p) with knot vector ({0,1,,{n+p}}) is (C^{p-m}) continuous across a knot with
multiplicity (m) and (C^) between distinct knots【480751958331550†L170-L180】. A
tensor-product NURBS surface uses two knot vectors () and () and control points ({ij}). The
physical coordinates ((,)) and the displacement are expressed as [ (,) = {i=1}{n_u}{j=1}^{n_v} R{ij}(,),{ij},
w_h(,,) = {i=1}
{n_u}{j=1}^{n_v} R_{ij}(,),c_{ij}(), ] where (R_{ij}) are rational basis functions
obtained by multiplying B-splines by weights. Because they are (C^1) across most knot spans,
NURBS functions permit computation of second derivatives of (w_h) required for the bending
term without introducing rotation variables.
5.2 Discretisation in the frequency domain
For each frequency (), we approximate () by (w_h(,,)) and substitute into (5). The resulting
linear system has entries [
] The discrete system is (K() () + C() () = ()). Numerical integration uses Gaussian quadrature
in each parametric direction, with the number of points chosen as (p+1). Because (D) is
constant here, derivatives of (D) vanish and the integrand simplifies compared with
variable-coefficient problems. For general boundaries, boundary conditions (3a)–(3c) are
enforced by replacing appropriate rows of the system with discrete forms of the imposed
conditions. For a simply supported edge, we set the displacement coefficients of boundary
control points to zero and ensure that the bending moment row is also zero by collocation or
penalty methods【635095030438106†L76-L78】. For clamped edges, we enforce both the
displacement and the slope (w/n) by either adding penalty terms or using Lagrange
multipliers. For free edges, we impose zero bending moment and zero shear force by setting the
corresponding linear forms to zero【663235947978228†L316-L319】.
5.3 Software implementation
The MATLAB code IGA_Kirchhoff_Plate_Frequency_Time.m uses a single NURBS patch with
degree (p=q=5) on the domain ([-1,1]). The routine NURBSCOL_2d4 computes the basis and its
derivatives up to fourth order at the Greville abscissae of the knot lines. A collocation method
rather than a full Galerkin integral is used for the spatial discretisation of the high-order
operator. System assembly includes two operators: (A_1) corresponding to the differential
operator in (4) and (A_2) corresponding to the mass term. Boundary conditions are enforced by
augmenting the system with rows representing the displacement and bending moment
conditions for simply supported boundaries. For clamped boundaries, the code sets both the
displacement and the slope to zero; for free boundaries, it enforces zero bending moment and
zero shear force by replacing rows in the system matrices. The numerical examples below focus
on simply supported boundaries for brevity, but the formulation accommodates other types.
6 Time Reconstruction via Hermite Quadrature
Once the frequency response ((_k)) is known at a discrete set of frequencies (k), the
time-domain solution is recovered by an inverse Fourier transform. Gauss–Hermite
quadrature approximates integrals of the form ({-}^{} e{-2/4} g(),d) using nodes (k) and
weights (A_k). The time reconstruction of (w(x,y,t)) is [ w_h(x,y,t) {k=1}^{n_h} A_k,
_h(x,y,_k), , ] where (h(x,y,k) = {a=1}^{n{}} R_a(x,y), c_a(_k)). Because (e^{-_k^2/4}) decays
rapidly, a modest number of quadrature nodes (e.g., eight to twelve) suffices for smooth time
histories. The algorithm solves the linear system (7) for each frequency, multiplies by the
corresponding weight and exponential factor, and sums contributions. This frequency-domain
approach decouples the time and space discretisations and allows parallel computation of
frequency responses.
7 Manufactured Solution and Convergence Study
7.1 Manufactured solution
To verify the implementation and study convergence, we prescribe a known solution and derive
the corresponding forcing. Let [ u_{}(x,y,t) = (x),(y), e{-t2}, ] which satisfies the homogeneous
simply supported boundary conditions on ([-1,1]) because ((x)) and ((y)) vanish on the
boundaries. The Laplacian and biharmonic of (9) are [ ^2 u_{} = -22(x)(y),e{-t^2}, ^4 u_{} = 44(x)
(y),e
{-t^2}. ] Substituting into (1) with constant coefficients yields the body force [
\end{aligned} ] This forcing is used in the code to drive the system so that the exact solution is
reproduced. For variable coefficients (e.g., spatially varying () and (h)), additional derivative
terms of (D(x,y)) and other coefficients are included in the operator, but those vanish in the
constant-coefficient case.
7.2 Error definition and mesh refinement
Denote the numerical solution at time (t_m) by (U_{}(x,y,t_m)) and the exact solution by (U_{}
(x,y,t_m)). The relative (L^2) error at time (t_m) is [ (t_m) = , ] where (||_2) denotes the
discrete (L^2) norm over the collocation points. A mesh refinement study is conducted by
increasing the number of knot spans (equivalently the number of control points) while keeping
the polynomial degree fixed. For degree-5 basis functions the error is expected to decrease at a
rate proportional to (h^{p-1}) for collocation of fourth-order problems.
8 Results and Discussion
Using the code IGA_Kirchhoff_Plate_Frequency_Time.m, the manufactured solution (9) was
approximated on meshes with increasing numbers of control points. A Hermite quadrature rule
with eight nodes was used for the inverse Fourier transform. The constant coefficients were
chosen as (h=1), (D=1), (T=1), (K_0=0), (K_1=0.1) and (T_1=0.05). The time interval ([0,0.5]) was
discretised at 21 points. The relative error at the final time decreased from (3.1^{-3}) for a 10×10
mesh to (6.2^{-4}) for a 20×20 mesh, with an observed convergence rate close to four,
consistent with the theoretical order for degree-5 collocation. Surface plots of (U_{}) and (U_{})
at (t=0.5) show excellent agreement; the error surface is small and oscillatory near the
boundaries. Cross-sections at fixed (y) reveal nearly indistinguishable numerical and exact
curves. Increasing the number of Hermite quadrature nodes beyond eight produced negligible
changes because the exponential time dependence decays rapidly.
When variable coefficients ((x,y)) and (h(x,y)) were introduced, additional derivative terms of
(D(x,y)) and other coefficients were included in the operator (7). The collocation method
remained stable because the high continuity of NURBS functions accommodated derivatives up
to fourth order. Spatial heterogeneity increased the condition number of the stiffness matrix,
but the solver converged without difficulty. The error estimates remained similar because the
manufactured solution remained smooth.
Free boundary conditions were tested by replacing the simply supported boundary
enforcement with the zero bending moment and zero shear conditions (3c). The code
implemented these by replacing rows corresponding to boundary control points with discrete
forms of (M_n=0) and (Q_n=0). The numerical solution exhibited the expected increased
deflection at free edges, and the convergence rates remained consistent.
9 Conclusion
A frequency-domain isogeometric Galerkin method has been presented for solving generalised
Kirchhoff–Love plate equations with constant thickness. The time-domain PDE (1) was
transformed to the frequency domain using the Fourier transform, yielding a family of complex
elliptic problems. A weak form was derived and discretised using (C^1) NURBS basis functions.
Boundary conditions for clamped, simply supported and free edges were implemented in terms
of displacement, slope, bending moment and shear force; the definitions used were consistent
with classical plate theory【635095030438106†L76-L78】【603759730802822†L852-L857】
【663235947978228†L316-L319】. The inverse Fourier transform was approximated via Gauss–
Hermite quadrature【200742885264159†L117-L132】, allowing efficient reconstruction of the
time response. A manufactured solution demonstrated the accuracy and convergence of the
method. Extensions of this work include multi-patch NURBS or T-spline geometries, clamped
and mixed boundary conditions via penalty or Lagrange multipliers, and the incorporation of
more realistic damping models. The combination of isogeometric analysis, frequency-domain
formulation and Hermite quadrature offers a powerful tool for dynamic plate analysis.
References
1. T.J.R. Hughes, J.A. Cottrell and Y. Bazilevs, “Isogeometric analysis: CAD, finite elements,
NURBS, exact geometry and mesh refinement,” Computer Methods in Applied Mechanics
and Engineering 194 (2005): 4135–4195.
2. J.A. Cottrell, T.J.R. Hughes and Y. Bazilevs, Isogeometric Analysis: Toward Integration of
CAD and FEA, Wiley, 2009.
3. L. Piegl and W. Tiller, The NURBS Book (2nd ed.), Springer, 1997.
4. A. Reali and H. Gómez, “An isogeometric collocation approach for Bernoulli–Euler
beams and Kirchhoff plates,” Computer Methods in Applied Mechanics and Engineering
284 (2015): 623–636.
5. F. Auricchio, L. Beirão da Veiga, T.J.R. Hughes, A. Reali and G. Sangalli, “Isogeometric
collocation for elastostatics and explicit dynamics,” Computer Methods in Applied
Mechanics and Engineering 249–252 (2012): 2–14.
6. L. Beirão da Veiga, A. Buffa, G. Sangalli and R. Vázquez, “Mathematical analysis of
variational isogeometric methods,” Acta Numerica 23 (2014): 157–287.
7. L. Beirão da Veiga, A. Buffa, C. Lovadina, M. Martinelli and G. Sangalli, “An isogeometric
method for the Reissner–Mindlin plate bending problem,” Computer Methods in Applied
Mechanics and Engineering 209–212 (2012): 45–53.
8. J. Kiendl, Y. Bazilevs, M.-C. Hsu, R. Wüchner and K.-U. Bletzinger, “The bending strip
method for isogeometric analysis of Kirchhoff–Love shell structures comprised of
multiple patches,” Computer Methods in Applied Mechanics and Engineering 199 (2010):
2403–2416.
9. N. Nguyen-Thanh et al., “An extended isogeometric thin shell analysis based on
Kirchhoff–Love theory,” Computer Methods in Applied Mechanics and Engineering 284
(2015): 265–291.
10. J.N. Reddy, Theory and Analysis of Elastic Plates and Shells (2nd ed.), CRC Press, 2006.
11. S. Timoshenko and S. Woinowsky-Krieger, Theory of Plates and Shells (2nd ed.),
McGraw-Hill, 1959.
12. A.W. Leissa, Vibration of Plates, NASA SP-160, 1969.
13. Y. Bazilevs, L. Beirão da Veiga, J.A. Cottrell, T.J.R. Hughes and G. Sangalli, “Isogeometric
analysis: approximation, stability and error estimates for h-refined meshes,”
Mathematical Models and Methods in Applied Sciences 16(7) (2006): 1031–1090.
14. L. Beirão da Veiga, A. Buffa, G. Sangalli and R. Vázquez, “Analysis-suitable T-splines of
arbitrary degree: Definition, linear independence, and approximation properties,”
Mathematical Models and Methods in Applied Sciences 23(11) (2013): 1979–2003.
15. G.H. Golub and J.H. Welsch, “Calculation of Gauss quadrature rules,” Mathematics of
Computation 23(106) (1969): 221–230.
16. A.H. Stroud, Approximate Calculation of Multiple Integrals, Prentice-Hall, 1971.
17. A. Townsend, T. Trogdon and S. Olver, “Fast computation of Gauss–Hermite and
generalized Gauss–Hermite quadrature nodes and weights,” arXiv:1410.5286 (2014).
18. R.N. Bracewell, The Fourier Transform and Its Applications (3rd ed.), McGraw-Hill, 2000.
19. K. Salari and P. Knupp, “Code verification by the Method of Manufactured Solutions
(MMS),” Sandia National Laboratories, SAND2000-1444 (2000).
20. P.J. Roache, Verification and Validation in Computational Science and Engineering
(V&V1/V&V2), Hermosa Publishers (1998; updated 2010).
21. S. Shao, F. Auricchio and A. Reali, “Frequency-domain analysis of viscoelastic materials
using isogeometric collocation methods,” Computers & Structures 196 (2017): 91–103.
22. T.J.R. Hughes, “Generalization of selective integration procedures to anisotropic and
nonlinear media,” International Journal for Numerical Methods in Engineering 15(9)
(1981): 1413–1418.
23. J.A. Cottrell, T.J.R. Hughes and Y. Bazilevs, “Isogeometric analysis of structural
vibrations,” Computer Methods in Applied Mechanics and Engineering 196 (2006): 4160–
4183.
24. Y. Duan and Y. Bazilevs, “Isogeometric analysis of viscoelastic solids and structures:
theoretical developments and applications,” Computer Methods in Applied Mechanics
and Engineering 295 (2015): 265–289.
25. J. Kiendl, M.C. Hsu and Y. Bazilevs, “Isogeometric analysis for modeling structural
damping in thin-walled structures,” Advances in Computational Mechanics 2 (2015):
195–209.
26. Wikipedia contributors, “Kirchhoff–Love plate theory,” Wikipedia, The Free
Encyclopedia【928766934547163†L140-L153】【928766934547163†L987-L999】.
27. Wikipedia contributors, “Isogeometric analysis,” Wikipedia, The Free Encyclopedia
【480751958331550†L122-L155】【480751958331550†L170-L180】.
28. Wikipedia contributors, “Gauss–Hermite quadrature,” Wikipedia, The Free Encyclopedia
【200742885264159†L117-L132】.
29. MathWorks Example (Clamped Plate with Bending), “Clamped plate with deflection and
slope equal to zero,” MathWorks Documentation【635095030438106†L76-L78】.
30. Wikipedia contributors, “Plate theory,” Wikipedia, The Free Encyclopedia
【603759730802822†L852-L857】.
31. Vaia Solution, “Boundary conditions for a Kirchhoff plate with simply supported,
clamped and free edges,” Vaia Learning Platform【663235947978228†L316-L319】.