Kirchhoff  Love Plate Theory: First-Order Analysis, Second-Order Analysis, Plate Buckling Analysis, and Vibration Analysis Using the Finite Difference Method

: This paper presents an approach to the Kirchhoff  Love plate theory (KLPT) using the finite difference method (FDM). The KLPT covers the case of small deflections, and shear deformations are not considered. The FDM is an approximate method for solving problems described with differential equations. The FDM does not involve solving differential equations; equations are formulated with values at selected points of the structure. Generally in the case of KLPT, the finite difference approximations are derived based on the fourth-order polynomial hypothesis (FOPH) and second-order polynomial hypothesis (SOPH) for the deflection surface. The FOPH is made for the fourth and third derivative of the deflection surface while the SOPH is made for its second and first derivative; this leads to a 13-point stencil for the governing equation. In addition, the boundary conditions and not the governing equations are applied at the plate edges. In this paper, the FOPH was made for all of the derivatives of the deflection surface; this led to a 25-point stencil for the governing equation. Furthermore, additional nodes were introduced at plate edges and at positions of discontinuity (continuous supports/hinges, incorporated beams, stiffeners, brutal change of stiffness, etc.), the number of additional nodes corresponding to the number of boundary conditions at the node of interest. The introduction of additional nodes allowed us to apply the governing equations at the plate edges and to satisfy the boundary and continuity conditions. First-order analysis, second-order analysis, buckling analysis, and vibration analysis of plates were conducted with this model. Moreover, plates of varying thickness and plates with stiffeners were analyzed. Finally, a direct time integration method (DTIM) was presented. The FDM-based DTIM enabled the analysis of forced vibration of structures, with damping taken into account. In first-order, second-order, buckling, and vibration analyses of rectangular plates, the results obtained in this paper were in good agreement with those of well-established methods, and the accuracy was increased through a grid refinement.


Introduction
This paper describes the application of Fogang's model [1] based on the finite difference method (FDM), used for the EulerBernoulli beam, to the KirchhoffLove plate. The Kirchhoff-Love plate theory (KLPT) is a twodimensional mathematical model that is used to determine the stresses and deformations in thin plates subjected to forces and moments. This theory is an extension of EulerBernoulli beam theory and was developed in 1888 by Love using assumptions proposed by Kirchhoff [2]. KLPT is governed by the GermainLagrange plate equation; this equation was first derived by Lagrange in December 1811 in correcting the work of Germain [3] who provided the basis of the theory. The analytical approach to KLPT consists of solving the governing equations that are expressed via means of partial differential equations, and satisfying the boundary and continuity conditions. However, solving the twodimensional partial differential equations may be difficult, especially in the presence of axial forces, an elastic Winkler foundation, plate of varying thickness, etc. Numerical methods permit therefore to overcome solving the differential equations. A considerable volume of literature has been published on numerical methods for KirchhoffLove plate analysis. For rectangular plates, Navier [4] in 1820 introduced a simple method for the analysis when a plate is simply supported along all edges; the applied load and the deflection were expressed in terms of Fourier components and double trigonometric series, respectively. Another approach was proposed by Lévy [5] in 1899 for rectangular plates simply supported along two opposite edges; the applied load and the deflection were expressed in terms of Fourier components and simple trigonometric series, respectively. More recently, Kindelan et al. [6] presented a method to obtain optimal finite difference formulas which maximize their frequency range of validity. Both conventional and staggered equispaced stencils for first and second derivatives were considered. Onyia et al. [7] presented the elastic buckling analysis of rectangular thin plates using the single finite Fourier sine integral transform method. Pisacic et al. [8] developed a procedure of calculating deflection of rectangular plate using a finite difference method, programmed in Wolfram Mathematica; the system of equations was built using the mapping function and solved with solve function. Li et al. [9] applied the generalized finite difference method to simulate the bending behavior of functionally graded plates; the governing equations and constrained boundary conditions were derived based on the first-order shear deformation theory and Hamilton's principle. Ferreira et al. [10] proposed a meshless strategy using the Generalized Finite Difference Method upon substitution of the original fourth-order differential equation by a system composed of two second-order partial differential equations; mixed boundary conditions, variable nodal density and curved contours were explored.
In the classical plate analysis using the FDM, the finite difference approximations are derived based on fourth order polynomial hypothesis (FOPH) and second-order polynomial hypothesis (SOPH) in each direction of the deflection surface; the FOPH is made for the fourth and third derivative of the deflection surface while the SOPH is made for the second and first derivative, leading to a 13-point stencil for the governing equation.. In addition, points outside the plate are generally not considered; the boundary conditions are applied at the plate edges and not the governing equations.
Consequently, the non-application of the governing equations at the plate edges together with the different polynomial hypotheses for the deflection surface have led to inaccurate results, making the FDM less interesting in comparison to other numerical methods such as the finite element method. In this paper, a model based on FDM was presented. This model consisted of formulating the differential equations with finite differences and introducing additional points at plate edges and at positions of discontinuity (continuous supports/hinges, incorporated beams, stiffeners, brutal change of stiffness's, etc.). The introduction of additional points allowed us to apply the governing equations at the plate edges and to satisfy the boundary and continuity conditions. Furthermore, the finite difference approximations were derived In these equations, w(x,y,z) is the displacement in z-direction, D is the flexural rigidity of the plate, E is the elastic modulus, d is the plate thickness, and  is the Poisson's ratio. The resulting equilibrium equation is given by (3) where q(x,y) denotes the applied transverse load per unit area. Substituting Equations (1a-c) into Equation (3) (4)) has fourth order derivatives; consequently, the deflection surface is approximated around a point of interest i as a fourth degree polynomial in each direction.
Thus, the deflection surface in x-direction i.e. can be described with values of deflections at grid points as follows: The shape functions f j (x) (j = i-2; i-1; i; i+1; i+2) can be expressed using the Lagrange interpolating polynomials: Applying the Lagrange interpolation polynomials (Equations (5-6)) the first, second, third, and fourth derivatives related to a five-point stencil can be expressed, respectively, as follows (7e-h) However, FDAs of the first and second derivatives are sometimes expressed using second order polynomial hypothesis in each direction of the deflection surface as follows (7i-j) Given the spacings x = h and y =h, the derivative  2 w/xy is then obtained using Equation (7j) as follows (7k)

Governing equation, moments, and Kirchhoff shear forces using a 13-point stencil
Equations (8a-b) also applied for the bending terms  4 w/x 4 and  4 w/y 4 . The twisting term  4 w/x 2 y 2 defined using Equation (7i) is expressed by means of the following 9-point stencil

Governing equation, moments, and Kirchhoff shear forces at nodes e, ee, s, and ss
The governing equation (Equation (4)) is described using the 13-point stencil of Equation (9g). The bending and twisting moments per unit length, and the Kirchhoff shear forces per unit length, with Poisson's ratio  = 0, are calculated using Equations (9i-m).

Governing equation, moments, and Kirchhoff shear forces at nodes eee and sss
The 25-point stencil and the 13-point stencil can be considered as described in Section 2.1.2.1. In case of a 25-point stencil, the governing equation is formulated using Equation (8e); the bending/twisting moments per unit length, and the Kirchhoff shear forces per unit length, with Poisson's ratio  = 0, are calculated using Equations (9a-e). In case of a 13point stencil, the governing equation is formulated using Equation (9g); the bending/twisting moments per unit length, and the Kirchhoff shear forces per unit length, with Poisson's ratio  = 0, are calculated using Equations (9i-m).

Finite difference approximations in the vicinity of free plate angles
Four additional nodes are introduced at the free plate angle, as represented in Figure 4. Alternatively, node nn could be considered instead of node ww, according to       (4), at any edge node is expressed using the 13-point stencil of Equation (9g).
The bending moments with respect to outer normal and tangential directions are expressed using Equations (1a-b), respectively, as follows The slope w/n of the deflection surface with respect to the outer normal direction, the curvatures - 2 w/n 2 and - 2 w/t 2 with respect to outer normal and tangential directions, respectively, are expressed using following formulas developed in Timoshenko et al. [11] (14a-c) The partial derivatives w/x and w/y in Equation (14a) are expressed using Equations (7d) or (7j),  2 w/x 2 and  2 w/y 2 using Equation (7i), and  2 w/xy using Equation (7k).
The Kirchhoff shear force, with respect to the outer normal to the skew edge, is determined as follows (15) where V x and V y are Kirchhoff shear forces with respect to Cartesian coordinates, expressed using Equations (9l-m).

Finite difference approximations of loadings
Let us determine here the FDA of the transverse load per unit area in the case of a varying distributed load q(x,y). Let (x i ,y i ) be the Cartesian coordinates of node i. This FDA, denoted by q i , can be taken as the average load around the node of interest and is then expressed as follows:    The governing equation, FDA of Equation (4), is applied at any regular node () of the plate. While applying the governing equation to the first two nodes at each end of the loading, the FDA of the twisting term will be formulated as follows (18) whereby the partial derivatives with respect to y will have non-centered FDAs according to Equation (7f). The derivatives with respect to y will have non-centered FDAs while applying Equation (19c-d) to the first two nodes at each end of the loading; Equations (7f) and (7i) will be then considered.
An adjustment of the continuity equations is made in case of a linear hinge (no continuity of the slope, m xx,il = m xx,ir = 0), a linear support (W il = W ir = 0, no equation (19d)), or a spring.

Monolithically connected beams or stiffeners
A beam in y-direction is assumed monolithically connected to the plate, its torsion stiffness being considered or not. Let linearly distributed forces (p) and moments (m  , positive counterclockwise) be applied on the beam.

Beam without torsion stiffness
The model represented in Figure 7a where EI b is the bending stiffness of the beam. Two boundary conditions must be satisfied at each beam's end; therefore, two additional nodes are introduced at each beam's end, with the deflections as unknowns. The term d 4 w b /dy 4 is expressed using Equation (8b).

Beam with Saint Venant torsion stiffness
The angles of twist of the beam are denoted by  b . Equations (19a-b) and (20a-b) are applied; in addition, the following equations are set where GI t is the torsional stiffness of the beam. Two boundary conditions on bending and one boundary condition on torsion must be satisfied at each beam's end; therefore, two additional nodes are introduced at each beam's end, one of which having two unknowns (deflection and angle of twist) and the other one unknown (deflection). The term d 2  b /dy 2 is expressed using Equation (7i).

Beam with Saint Venant torsion stiffness and warping torsion stiffness
Equations (19a-b), (20a-b), and (21a) are applied; in addition, the following equation is set where EI  is the warping torsion stiffness of the beam. Two boundary conditions on bending and two boundary conditions on torsion must be satisfied at each beam's end; therefore, two additional nodes are introduced at each beam's end, both of which having two unknowns (deflection and angle of twist). The term d 4  b /dy 4 and d 2  b /dy 2 are expressed using Equations (7a) and (7c), respectively.

Concentrated support or spring at node i
The unknown support reaction S i can be converted into a transverse load per unit area q i as follows where x and y are the grid spacing in x-and y-directions, respectively. So, there is one unknown (S i ) more and one more boundary condition (the deflection at node i is zero). As mentioned in Section 2.1.3.1, a local grid refinement would deliver more accurate results at the expense of intensive calculations. The case of a concentrated spring support with stiffness K is analyzed similarly, as follows (24)

Elastic Winkler foundation over an area of the plate
Given an elastic Winkler foundation having the stiffness k w . Equation (4) becomes (25) To account for the elastic Winkler foundation, k W h 4 /D is added to the term associated to the node of interest in the stencil describing the governing equations (Equations (8e), (9g), (10), and (11)). Equation (8e) i.e. becomes

Local grid refinement
A local grid refinement can be considered at positions where concentrated load or bending moment is applied, or at a concentrated support or spring. Refinement nodes () are then introduced around the node of interest, as represented in

Finite difference approximations for an isotropic Kirchhoff-Love plate 2.2.2.1 Finite difference approximations at an interior node
Let the nodes spacings be x = h and y = h in x-and y-direction, respectively. Following parameters are defined (29a-c) Similarly to Section 2.1.2.1, the analysis is conducted considering a 25-point stencil and a 13-point stencil.

Governing equation, moments, and Kirchhoff shear forces using a 13-point stencil
The stencil notation of the governing equation (Equation (27)) of the plate is as follows The bending/twisting moments per unit length are determined using Equations (9i-k), and the transverse forces per unit length using Equations (28a-b), (9l-m), (7j), and (29a-b), as follows (33)

Finite difference approximations in the vicinity of supported plate angles
The node distribution was represented in Figure 3.

Governing equation of the plate, bending moments, and twisting moments at node k
Applying the formulas developed in first-order analysis (Section 2.1. (27) yields the stencil notation of the governing equation (Equation (27)) of the plate as follows

Governing equation, moments, and Kirchhoff shear forces at nodes e, ee, s, and ss
The governing equation (Equation (27)

Governing equation, moments, and Kirchhoff shear forces at nodes eee and sss
The 25-point stencil and the 13-point stencil can be considered as described in Section 2.

Finite difference approximations in the vicinity of free plate angles
The node distribution was represented in Figure 4.

Governing equation of the plate, bending moments, and twisting moments at node k
Applying the formulas developed in first-order analysis (Section 2.1. (27) yields the stencil notation of the governing equation (Equation (27)) of the plate as follows

2.3) and Equations (29a-c) into Equation
The bending/twisting moments per unit length are determined using Equations (9i-k). The transverse force per unit length T x is determined using Equation (33a), and T y using Equation (12)

Governing equation, moments, and Kirchhoff shear forces at nodes e, ee, s, and ss
The governing equation, Equation (27), is given by Equation (32). The bending/twisting moments per unit length and the transverse forces per unit length are determined using Equations (9i-k) and (33), respectively.

Governing equation, moments, and Kirchhoff shear forces at nodes eee and sss
The dispositions of Section 2.2.2.2 applied.

Finite difference approximations at skew edges
The node distribution of Section 2. The transverse force T n , with respect to the outer normal to the skew edge, is determined as follows (37) where T x and T y are transverse forces with respect to Cartesian coordinates, expressed using Equations (33).

Analysis at positions of discontinuity
The dispositions of Section 2.2.2 applied; however the Kirchhoff shear forces are replaced by the transverse forces.

Vibration analysis of an isotropic Kirchhoff-Love plate 2.3.1 Free vibration analysis of a plate with constant thickness
Our focus here is to determine the eigenfrequencies of the plate. Damping is not considered. A second-order analysis is conducted; the first-order analysis can be deduced. The deflection surface is denoted by w  (x,y,t). The governing equation can be expressed as follows: where  is the plate's mass per unit volume and d is the plate thickness. A harmonic vibration being assumed, w * (x,y,t) can be expressed as follows: x y x y w x y t w x y t w x y t x x y y Here,  is the circular frequency of the plate. Substituting Equation (39a) into Equation (38) yields (39b) The FDAs of Equation (39b) are obtained by adding -d 2 h 4 /D to the term associated to the node of interest in the stencil describing the governing equations. Therefore, Equation (30) i.e. becomes The bending/twisting moments per unit length, the shear forces and transversal forces per unit length are calculated similarly to previous sections.

Effect of a concentrated mass
We analyzed the dynamic behavior of a plate carrying a concentrated mass at node i; the mass M p is defined as follows (41) where m p is the dimensionless mass, and x and y = x are the node spacings in x-and y-direction, respectively.
The FDA of the governing equation at node i is obtained by adding -m p d 2 h 4 /D to the term associated to the node of interest in the stencil.   The governing equation (Equation (44)) can be formulated with FDM at any point of the plate at any time t. As described earlier, additional points are introduced to satisfy the boundary and continuity conditions. Thus, the plate deflection w * (x,y,t) can be determined, and the efforts be calculated with appropriate formulas developed earlier.

Analysis of isotropic Kirchhoff-Love plates of variable stiffness
The flexural rigidity D of the plate may vary throughout the plate, and is then denoted D(x,y). Substituting D(x,y) into Equations (1a-c) yields the bending and twisting moments as follows (46) The substitution of Equation (46) into (3) yields the following governing equation The Kirchhoff shear forces are expressed using Equations (46) as follows ,  Results and discussions

Rectangular plate simply supported along all edges and subjected to a uniform load
A rectangular plate simply supported along all edges and subjected to a uniformly distributed load was analyzed. The plate dimensions in x-and y-directions are denoted by a and b, respectively, and the uniform load by q.
The bending moments at the middle of the plate, depending on the ratio b/a, are m xx,m = qa 2 /N x , and m yy,m = qa 2 /N y .
Details of the analysis and results are presented in Appendix A and in the supplementary material "Rectangular plate simply supported and subjected to a uniform load." Table 1 lists the results obtained with Timoshenko [11] using the solution by Lévy [5] (exact results) and those obtained in the present study.

Rectangular SSFF plate subjected to a uniformly distributed load
A rectangular plate simply supported along edges x = 0 and x = a, fixed along edges y = 0 and y = b, and subjected to a uniformly distributed load was analyzed. The bending moments, depending on the ratio a/b, are m yerm = -qb 2 /N yerm , m xm = qb 2 /N xm , and m ym = qb 2 /N ym ; m yerm is the moment at the middle of the fixed edge, m xm and m ym are the bending moments at the middle of the plate in x-and y-directions, respectively. Analysis and results are detailed in Appendix B and in the supplementary material "Rectangular SSFF plate subjected to a uniformly distributed load." Table 2 displays the results obtained with Courbon [13] using the solution of Lévy [5] (exact results) and those obtained in the present study.

Rectangular plate simply supported along all edges and subjected to a non-uniform heating
A rectangular plate simply supported along all edges and subjected to a non-uniform heating was analyzed. The plate dimensions in x-and y-directions were denoted by a and b, respectively. The bending moments at the middle of the plate, depending on the ratio a/b and Poisson's ratio  = 0, are m xx,m = D T N x , and m yy,m = D T N y . Detailed analysis and results are presented in Appendix C and in the Supplementary Material "Rectangular plate simply supported along all edges and subjected to a non-uniform heating." Table 3 lists the results obtained by Fogang [14] using the solution by Lévy [5] and the Fourier sine transform method, and those obtained in the present study.

Rectangular plate simply supported along all edges and subjected to a uniformly distributed load and a compressive force
A rectangular plate simply supported along all edges and subjected to a uniformly distributed load q and compressive forces N x was analyzed. The compressive forces were applied along edges x = 0 and x = a, whereby a and b were the plate dimensions in x-and y-directions, respectively. The axial force N x is such that N x a 2 / 2 /D = -1.0. The bending moments at the middle of the plate, depending on the ratio b/a, are m xx,m = qa 2 /C x , and m yy,m = qa 2 /C y . Details of the analysis and results are presented in Appendix D and in the supplementary material "Rectangular plate simply supported and subjected to a uniform load and a compressive force." Table 4 lists the results obtained with Timoshenko [11] using the solution by Navier [4] (exact results) and those obtained in the present study.

Plate buckling of a rectangular plate simply supported along all edges and subjected to a compressive force
The plate buckling analysis of a rectangular plate simply supported along all edges and subjected to compressive forces N x was conducted. The compressive forces were applied along edges x = 0 and x = a, whereby a and b were the plate dimensions in x-and y-directions, respectively. The critical axial force N x,cr is defined as follows N x,cr = -k 2 D/b 2 ; k is the plate buckling factor. Bryan [15] derived following expressions of the buckling load and buckling factor (49) Details of the results are presented in the supplementary material "Plate buckling of a rectangular plate simply supported and subjected to a compressive force." Table 5 displays the buckling factors, depending on the ratio a/b, obtained using Equation (49) by Bryan [15] and those obtained in the present study. Table 5 Buckling factor k for the plate subjected to compressive forces along edges x = 0 and x = a The results of the present study show good agreement with the results of Bryan [15], and the accuracy is increased through a grid refinement. However, special attention should be taken when solving the eigenvalue problem, since confusion between the modes is easily done; an efficient tool may be needed for this purpose.

Vibration analysis 3.3.1 Free vibration analysis of a rectangular plate simply supported along all edges
The free vibration analysis of a rectangular plate simply supported along all edges was conducted. The plate dimensions in x-and y-directions were denoted by a and b, respectively. The vibration frequency  is defined as follows: The circular frequency for this special case is widely known in the literature (see Wikipedia [16]); its general formulation and the first mode (m = n = 1) are given as follows (51) Details of the results are presented in the supplementary material "Free vibration analysis of a rectangular plate simply supported along all edges." Table 6 displays the free vibration factors  * (Equation (50)), depending on the ratio a/b, obtained using Equation (51b) by Wikipedia [16] and those obtained in the present study. Table 6 Coefficients  * of natural frequencies (first mode) of a rectangular plate The results of the present study show good agreement with the exact results, and the accuracy is increased through a grid refinement. However, special attention should be taken when solving the eigenvalue problem, since confusion between the modes is easily done; an efficient tool may be needed for this purpose.

Conclusions
The FDM-based model developed in this paper enabled, with relative easiness, first-order analysis, second-order analysis, and vibration analysis of KirchhoffLove plates. The results showed that the calculations conducted as described in this paper were accurate. The 13-pt stencil and 25-pt stencil presented deliver good results, whereby the results of 13-pt stencil were better for a given grid. A new 25-pt stencil should be investigated in in future research, whereby a fourth-order polynomial hypothesis for the deflection surface would be considered at plate angles and at skew edges.
The following aspects were not addressed in this study but could be analyzed with the model in future research: However, some study limitations should be acknowledged 

Large deformation theory
Supplementary Materials: The following files were uploaded during submission:  "Rectangular plate simply supported along all edges and subjected to a uniform load,"  "Rectangular SSFF plate subjected to a uniformly distributed load,"  "Rectangular plate simply supported along all edges and subjected to a non-uniform heating,"  "Rectangular plate simply supported and subjected to uniform load and compressive force,"  "Plate buckling of a rectangular plate simply supported and subjected to a compressive force,"  "Free vibration analysis of a rectangular plate simply supported along all edges." Only the 8 x 8 element discretization was uploaded in these supplementary files.

Author Contributions:
Funding: