The correct derivation of the buckling equations of the shear-deformable FGM plates for the extended Kantorovich method

This article presents the derivation of the elastic buckling equations and boundary conditions of shear-deformable plates in the frame of the extended Kantorovich method (EKM). Surveying the literature shows that those stability equations are often obtained using a wrong derivation by confusing them with the linear equilibrium condition. This work aims at providing the correct derivation that is built on the stability of the equilibrium condition. Buckling equations are derived for three different plate theories, namely, the first-order shear deformation plate theory (FSDT), the refined-FSDT, and the refined plate theory (RPT). This article is the first to implement the EKM based on a refined theory. Also, it is the first time to implement the refined-FSDT in buckling analysis. For the generic FGM plates, buckling equations derived based on the FSDT and refined-FSDT are both found to be simple and contain only the lateral displacements/rotations variations. On the other hand, those of the RPT, have coupled lateral and in-plane displacement variations, even if the physical neutral plate is taken as the reference plane. The considered plate is rectangular and under general in-plane loads. The properties of the plate are continuously varying through its thickness which is assumed to change smoothly with a separable function in the two in-plane directions. The von Kármán nonlinearity is considered. The stability equations are derived according to the Trefftz criterion, using the variational calculus. The solution methods of the obtained equations are out of the scope of this article, however, a brief on the solution strategy is presented.


Evolution of the EKM
The extended Kantorovich method (EKM) had evolved from the so called classical Kantorovich method that was introduced by Kantorovich [12] as a solution method to the partial differential equations (PDE) by reducing the d-dimensional PDE that governs the problem into d one dimensional (1D) ordinary differential equations (ODE). Obtaining those ODE's starts by approximating the d-dimensional solution function Wðx 1 ; x 2 ; . . .x d Þ as the multiplication of d uni-dimensional functions w ðkÞ , one in each dimension x k , as shown in Eq. 1.
There are two approaches to obtain the ODE's from the original PDE using the approximated solution function. Both start with the assumption that only one uni-directional function in the approximate solution is unknown, while all the others are given. The first approach is to substitute the approximate solution to the total potential energy integral, then use the variational calculus to derive the governing equation, which is now an ODE. In the second approach, as presented by Grimm and Gerdeen [8], the approximate solution is substituted to the Galerkin form, i.e., the weak form, of the governing PDE and then the fundamental lemma of the variational calculus is implemented to obtain the ODE. Solving the obtained ODE terminates the steps of the classical Kantorovich method. The limitation of the classical Kantorovich method is that its accuracy depends on the initially guessed functions.
To overcome this limitation, Kerr [13] proposed the extended Kantorovich method (EKM), in which the obtained ODE's are solved iteratively until a satisfying convergence is achieved. In the EKM, the first guessed solution, i.e., the trail function, does not have to satisfy any of the boundary conditions. This is the main advantage of the EKM over the other variational methods like Rayleigh-Ritz and Galerkin weighted residual methods. EKM is often reported to provide accurate results with rapid convergence with minimum dependency on the trail function.
Yuan [38] further extended the EKM by introducing the use of multi-term trail functions, as shown in Eq. 2, which describes the solution more accurately.
where n is the number of terms. The multi-term EKM provides more accuracy and the ability to solve additional problems that were not solvable with the single-term EKM, e.g., the buckling of rectangular plate under pure in-plane shear loading, see [9,37]. By using n terms in the approximate solution, the ddimensional PDE is reduced to d sets, each contains n ODE's that have to be solved simultaneously.

Applications of the EKM to the plate buckling problem
Among the first applications of the EKM included the plate bending problem [15], and the plate's eigenvalue problems, i.e., vibration [14] and buckling analyses [1,8,14]. Those researches had been followed by numerous ones in which the EKM was implemented in investigating the plate problems considering various complexities and configurations. A recent review of those articles is presented by Singhatanadgid [33]. It is observed that the vast majority of the researches that applied EKM on the buckling analysis of plates have considered the thin plates, and thus were based on the classical plate theory (CPT). Just a few are based on a shear deformation theory, namely, [17,[26][27][28]. Also, almost all the studies investigated the rectangular plates, except [36,41] in which annular plates are investigated, and [9,10] in which skew plates are considered. The buckling of plates having a variable thickness is investigated using EKM in [4,27]. Buckling of plates under linearly distributed in-plane forces is investigated using EKM in [8,16,18,28], while the case of general in-plane loading is considered in [30]. It is also observed that all of the publications in this topic investigated the buckling under mechanical loading, except [10] in which the thermal loading is considered. The effect of the elastic foundation on the buckling of plates is considered in [8][9][10]. The stability equations in [8, 16-18, 35, 36, 41] are obtained using the Galerkin weighted residual method mentioned in Section 1.1, in which the approximate solution is substituted to the weak form of the twodimensional stability equations. The other researches, i.e. [4,10,[26][27][28][29][30][31], implemented the energy approach of substituting the approximate solution to the integral of the total potential energy then deriving the stability equations using the variational calculus.

Stability and Trefftz Buckling Criterion
Consider two rigid balls at rest, one on the top of a hill and the other on the bottom of a valley, as shown in Fig. 1. Both balls are said to be in equilibrium, i.e., no variation in the potential energy; i.e., where C is the potential energy; and d is the variational operator. Now consider an infinitesimal push is applied to each balls to its left side. The one resting on the bottom of the valley, ball B, will climb up; while the one resting on the top of the hill, ball A, will fall down the hill.
Once the applied push is off, ball B will roll back to its first position, while ball A will continue its drop until it finds another resting position. To differentiate between these two different behaviors, ball B is said to be in a stable equilibrium, while ball A is said to be in an unstable equilibrium. From this classical example, stability can be defined as the persisting of the equilibrium configuration under small disturbance, while the instability is the change of the equilibrium state under small disturbance.
If much bigger side push is applied on the so said stable ball B, it will pass the top and fall down the hill. So, the stable ball B has a limit of the disturbance after which it loses its stability and tries to find another resting position different from the one it started on.
Before reaching the stability limit, the mechanical structure is said to be in equilibrium and stable, i.e., dC ¼ 0 and DC [ 0, while after the stability limit the structure is in equilibrium but unstable, i.e., dC ¼ 0 and DC\0. The boundary point between these two configurations is the stability limit, at which dC ¼ 0 and DC ¼ 0.
The main criteria on which the stability equations can be built are the neutral equilibrium method and the adjacent equilibrium method. The latter criterion is also known as the principle of minimum potential difference or the Trefftz criterion. Both criteria are based on the definition of stability limit as at which the equilibrium configuration of the structure can change under the same total potential energy, i.e., DC ¼ 0. Trefftz criterion finds the stability limit by examining two infinitesimally close different equilibrium states. The total potential energy at the two neighbour equilibrium states can be written as shown in Eq. 4.
where the subscript is put to differentiate between the two adjacent configurations. If DC is expanded in Taylor series, then Eq. 4 can be rewritten as follows.
Neglecting the higher order terms and then taking the first variation of the two sides gives the following equation.
Since both of the configurations are in equilibrium, i.e., dC ¼ 0 and dC 1 ¼ 0, then Eq. 6 reduces to Equation 8 is called the Trefftz criterion, and points to the stationary state of the second variation of the potential energy at the equilibrium state prior to buckling. Herein, Trefftz criterion is written as where d 2 e C is just a compact symbol that means d 2 Cj dC¼0 .

Necessity of this derivation
The important observation that necessitates the writing of this article is that except for the single article by Hassan and Kurgan [10], all the other publications those presented an energy approach based derivation of the buckling equations for EKM had wrongly confused the stability equations with the equilibrium equations. In all those publications, the stability equations had been derived from the state of stationary potential energy, i.e., the equilibrium condition that is shown in Eq. 3. They had derived the linear equilibrium equations and simply solve them as if they were eigenvalue stability equations. However, according to Jones [11], this is fundamentally wrong.
According to the Trefftz criterion, stability equations are derived from the stationary state of the second Fig. 1 Equilibrium vs stability variation of the potential energy at the equilibrium state prior to buckling, as shown in Eq. 8. Detailed discussion on the derivation the stability equation for plates is presented by Jones [11]. An important note is that if wrong notations and misinterpretation were to be tolerated, then the wrongly derived stability equations in the literature would resemble the correct stability equation. So, the solutions of those wrongly derived stability equations are not necessarily wrong.
This article aims at presenting the correct derivation of the static stability (buckling) equations of plates using the EKM. The derivation is provided in detail so fellow researchers can easily redo the derivation for other configurations of their concern. The problem configuration presented here is selected to be reasonably general, so many other interesting problems can be obtained by simplifying the results. The plate is selected to be of rectangular shape. The von Ka´rma´n geometrical nonlinearity is considered. The material of the plate is considered as functionally graded material (FGM), with arbitrary continuously changing properties.
Unlike the derivation presented by Hassan and Kurgan [10] that is based on CPT, the ones presented here are based on the shear-deformation plate theories. Three different theories are adopted. The first is the firstorder shear deformation theory (FSDT). The FSDT is considered because, as stated by Reddy [23] and Qi [22], there is no significant benefit of using higher-order shear deformation theories in the buckling analysis of plates. The second adopted theory is the refined-FSDT first proposed by Qi [22], and recently modified by Nguyen [19]. The third one is the refined plate theory (RPT), first proposed by Shimpi [25], and numerous versions are presented ever since.
To the best of the authors' knowledge, this is the first article to implement EKM based on any refined plate theory. Also, it is the first time to implement the refined-FSDT in any buckling analysis. The first derivation that is based on FSDT is presented in detail, while the latter ones are briefly presented due to the similarity in steps with the first.
One last note to mention is that this article does not deal with the solution of the obtained ODE's. Solution of those equations is a separate topic. Surveying the literature shows that analytical, numerical and semianalytical methods can all be used for that purpose.

The derivation based on the FSDT
The stability equations are derived here according to the Trefftz criterion, shown in Eq. 8. The reasonable steps are to first obtain the total potential energy C of the plate. Then to find the its second variation d 2 C, and impose the equilibrium conditions to the second variation, which gives d 2 e C. Next is to impose the concept of the Kantorovich method, which makes d 2 e C a single-integral functional. The next step is to obtain the variation of to zero to obtain the stability equations and boundary conditions.

Finding C
C of a rectangular plate having length of 2a, width of 2b and thickness h, is given by where U is the strain energy, and V is the potential energy of the applied loads. U and V are given as follows.
Àh=2 r xx xx þ r yy yy þ r zz zz where r ii and s ij are the normal and shear stresses, while e ii and c ij are the normal and shear strains, and N ij are the in-plane applied loads at the reference plane, in the stated directions. In order to find those stress and strain components, the constitutive equations of FSDT are adopted.
where u, v, and w are the displacements of any point of the plate in the directions x, y, and z, respectively. h and u are the rotations of a transverse normal about yÀaxis and xÀaxis, respectively. z 0 is the offset between the mid-plane and some reference plane through the thickness of the plate. u, v, and w are the displacements of any point at the physical neutral plane. Considering the von Ka´rma´n nonlinearities, strains are related to the displacements by the following equations.
Note that the subscripts after commas denote the partial derivatives in the stated directions. Stresses are related to the strains as follows.
where E z ð Þ is Young's elasticity modulus as a function of z. l is the Poisson's ratio. k s is the shear correction factor that accommodates the discrepancy in the energy caused by the FSDT assumption of constant transverse shear stresses [23]. Substituting the constitutive Eqs. 12 to 13 and then to Eq. 14 reveals the stress and strain components in terms of u, v, w, h, u, and their derivatives. Substituting those stresses and strains to Eq. 9 then integrating along z gives C as a double integral functional of u and u x;y ð Þ .
where A, B, and D are given as follows.
D and A are the bending and extension stiffness. For the general case of variable thickness h across the plate, D and A can be expressed as separable functions as shown below.
where D 0 and A 0 are constants; D x and A x are functions of x; D y and A y are functions of y. If the thickness h is constant, , the second variation of the total potential energy, is obtained as follows.
where du, dv, dw, dh, du, du ;x , du ;y , dv ;x , dv ;y , dw ;x , and dw ;y are the variational displacements and their derivatives. Note that all the terms of V vanish when taking the second variation of the total potential energy. So, only U contributes to d 2 C.

Finding d 2 e C
The equilibrium of the flat plate implies that no lateral displacement or rotation occurs prior to the buckling point. Therefore, eliminating w, h, u, and their derivatives from d 2 C in Eq. 18 gives d 2 e C, which is the second variation of the potential energy evaluated at the equilibrium state. where In Eq. 19, the terms multiplied by the integral B are the only ones that couple the variations of the lateral deflection and rotations (dw, dh, du), with those of the in-plane displacements (du; dv). Equating the integral B to zero eliminates that coupling. From Eq. 16, B is controlled by the material properties E ðzÞ and l, and the the offset z 0 of the reference plane from the geometrical mid-plane through the thickness direction. The distance z 0 at which B ¼ 0 defines the neutral plane, at which no stresses or strains occur when the plate is under pure bending [39,40]. Note that when E z ð Þ is constant or symmetric about the mid-plane then z 0 has to vanish and the physical neutral plane coincides with the mid-plane.
Back to Eq. 19, the in-plane displacement terms (u, v) multiplied by ðdw ;x Þ 2 , ðdw ;x Þ 2 and dw ;y dw ;x present nothing but the integrals of the equilibrium state stress components r xx , r yy and s xy , respectively. Those integrals are called the force in-plane resultants N xx , N yy , and N xy , and given by integrating the stress components of Eq. 14 along the thickness. Using the notation shown in Eq. 16, the force in-plane resultants are as follows.
The resultant in-plane forces in the general case are variable through the plate's in-plane dimensions. They can be generally expressed in separated form as follows.
For the case that the plate is under uniform or linearly varying in-plane loads acting in the principle directions x and y, then the resultant in-plane forces will be equal to those applied forces [16,28]. Otherwise, if the applied in-plane loads have nonlinear distributions, then the resultant in-plane forces have to be obtained by solving the pre-buckling equilibrium problem [30]. The resultant forces do not have to be produced by direct mechanical load, as they may be created in plates under thermal loads, given that there is a restriction on expansion/retraction in one direction at least. For more about how to obtain the resultant thermal loads in plates see [6,10]. d 2 e C can be rewritten as follows. where In Eq. 22, the terms in d 2 e C are rearranged in two groups. The first, C L , contains the terms of variation of lateral displacement and rotation and no terms of variation of the in-plane displacements. While the second group, C I , contains only the terms of variation of the in-plane displacement, and has no lateral displacement or rotation terms. In addition, d 2 e C L contains the force resultant components, while the d 2 e C I doesn't. Therefore, a further variation of d 2 e C will result in an uncoupled set of equations. So, it is safe to only consider C L in the remaining steps of the derivation while being sure that this will not affect the resulting stability equations. The considered part in the remaining of the derivation is named d 2 e C L .

Imposing the concept of the Kantorovich method
In this step, the unknown bivariate functions dw ðx;yÞ , dh ðx;yÞ , and du ðx;yÞ are expressed as the summation of many terms, each is a multiplication of univariate functions of different variables, as illustrated in Eq. 2. dw ðx;yÞ , dh ðx;yÞ , and du ðx;yÞ are approximated as follows.
where f , g, and m each is a single-row matrix contains the functions f i ðxÞ, g i ðxÞ and m i ðxÞ, respectively. F, G, and M each is a single-column matrix contains the functions F i ðyÞ, G i ðyÞ and M i ðyÞ, respectively. Note that the function dw ðx;yÞ is approximated using n w terms, while n h terms are used for dh ðx;yÞ , and n u for du ðx;yÞ . n w , n h , and n u are all independent of each other; i.e., each of the functions can be approximated using different number of terms. Substituting the approximations in Eq. 23 to C L in Eq. 22 gives where Assuming the functions of one direction are all known and differentiable, say F, G, and M, then the 2D functional d 2 e C L in Eq. 24 is reduced to a 1D functional of the functions f , g, and m. The stability equations are now obtainable by substituting d 2 e C L to the Trefftz criterion in Eq. 8.
where H's are as follows.
Implementing the integration by parts on each integral in Eq. 25 that contains df 0 , dg 0 , or dm 0 , gives: Since df , dg, and dm are all arbitrary, and to satisfy Eq. 27, all the following stability equations must hold.
In addition, the following boundary conditions must be satisfied at x 2 fÀa; þag.
When both the lateral displacement and rotations are restrained, the edge is said to be clamped ð CÞ. The simply supported edge ð SÞ has the same restraints as the clamped one except that it can freely rotate about its axis. When there is no restriction on both the lateral displacement and rotations, the edge is called free ð FÞ. Those boundary conditions at the edges x 2 fÀa; þag are translated as follows. C : f ¼ 0 and g ¼ 0 and m ¼ 0 S : The buckling loads are obtained by solving the governing Eqs. (28 and 29) as an eigenvalue problem for the critical resultant forces N xx , N yy , and/or N xy .
This completes the derivation of the stability equations and boundary conditions for a given solution in the x direction. The second set of governing equations and boundary conditions for a given solution in the y direction can be obtained following the same steps.

Stability equations and boundary conditions for F ðyÞ , G ðyÞ , and M ðyÞ
To obtain the stability equations and boundary conditions for the other direction, the steps are to be followed starting from Eq. 25 and assuming that f , g and m are the known functions. Substituting this assumption to the Trefftz criterion in Eq. 8 gives: where T's are as follows.
Implementing the integration by parts on each integral in Eq. 31 that contains dF 0 , dG 0 , or dM 0 , gives: Since dF, dG, and dM are all arbitrary, and to satisfy Eq. 33, all the following stability equations must hold.
In addition, the following boundary conditions must be satisfied at y 2 fÀb; þbg.
From Eq. 35, the boundary conditions at the edges y 2 fÀb; þbg can be written as follows C : Note that the obtained FSDT stability equations (28 and 34), and boundary conditions (29 and 35), resemble those found in the available literature, namely [26][27][28] if wrong notations and misinterpretation were to be tolerated. The difference is that the ones presented here are correctly derived based on a stability criteria, while those found in [26][27][28] are wrongly obtained from the equilibrium state by vanishing the first variation of the potential energy. The equations obtainable from the equilibrium state are boundary value problems, while those obtainable based on a stability criterion are eigenvalue problems. The equilibrium equations contain displacement and load terms, while the buckling equations contain the variations of those terms Note that the equations derived in [26][27][28] are for isotropic plates, while those in this article are for a more general case of FGM plates. However, they are comparable to each other by reducing those derived in this article to the special case of isotropic plates. This reduction is obtainable by simply making the elasticity modulus a constant, i.e., E ðzÞ ¼ E, in Eq. 16. Authors have had to compare the obtained equations to the only available isotropic case. As mentioned in Sect. 1.2, the only available works that implemented EKM in buckling analysis of shear deformable plates are Refs. [17,[26][27][28]. Lopatin [17] treats the case of composite sandwich plate, which is different from being continuously changing properties. This leaves out Refs. [26][27][28] to be the only comparable works to those in this article.
Many other aspects can be easily included to the problem, for example, the effects of the elastic foundation, and the skew angle of the plate's geometry. Those aspects are not added to the derivation in this brief because they just add terms to the equations without contribution to the main aim of clarifying the derivation steps.

The derivation based on the refined-FSDT
The idea behind the refined-FSDT that first proposed by Qi [22], is to replace the shear correction factor k s with a suitable shear correction functions imposed to the transverse shear stress and strain in order to dictate their distributions through the thickness. Recently, Nguyen [19] merged that same concept with the one of decomposing the lateral displacement w into bending and shear components (w b ,w s ) to get rid of the dedicated rotation variables, (h, u). The constitutive equations of the refined-FSDT are as follows.
As been mentioned above, the transverse shear strains are corrected by the function g ðzÞ which gets transferred to the transverse shear stresses through the kinematic relations, and no additional correction factor k s is required. Thus transverse shear stresses and strains are given as follows.
Note that the refined-FSDT herein refers to the general form of the theory, so, g ðzÞ used here is general and does not have to be similar to the one provided by Nguyen [19].
Due to the similarity of the steps, and to keep the brevity, the detailed steps of the derivation are omitted, and only a brief description is provided. Starting from the constitutive Eq. 37 considering the new formulae of the transverse shear stresses and strains in Eq. 38, and by implementing the same steps in Sect. 2 of finding the total potential energy C, then its second variation d 2 C, then imposing the prebuckling equilibrium state of zero lateral deflection (w b ¼ w s ¼ 0) to obtain d 2 e C and then approximating (dw b ,dw s ) as then implementing the integration by parts, given that F and R are known, produce the following buckling equations.
The associated boundary conditions that must be satisfied at x 2 fÀa; þag are given below. C : where J's and k's are as follows.
where Q 0 is constant, Q x is a function of x and Q y is a function of y that their multiplication produces the following integral, given that the thickness h is changing through the x and y directions.
The second set of governing equations and boundary conditions can be derived similarly. For the sake of brevity, and due to the similarity with the ones above, those resulting formulations are omitted.

The derivation based on the RPT
The sole aim of this section is to illustrate that for the FGM plates, using the refined plate theories (RPT) provides stability equations having coupled lateral and in-plane displacements variations, even if the physical neutral plane is considered as the reference plane. In RPT, the lateral displacement w is decomposed into into w b and w s . RPT treats the shear correction from down the level of the constitutive equations. The shear component of the lateral displacement imposed with a selected function n ðzÞ contributes to the in-plane displacements. n ðzÞ gets transferred to the transverse shear strains and stresses through the kinematic relations. There are numerous versions of the traditional refined plate theory, e.g., [2,7,25,34]; they normally differ in the adopted shear distribution function, n ðzÞ . Note that RPT herein refers to the form of the theory n ðzÞ is general function. The general form of the constitutive equations of RPT is given as follows.
In the RPT there is no need for the correction factor k s , so the transverse shear stresses and strains are given as follows.
Starting from these constitutive equations and transverse shear strains and by implementing the same steps in Sect. 2 of finding the total potential energy C, then its second variation d 2 C, then imposing the prebuckling equilibrium state of zero lateral deflection (w b ¼ w s ¼ 0), provide d 2 e C as follows.
where where the integrals D, A, and B are the same ones given in Eq. 16, while the integrals Q 1 :::Q 4 are given as follows.
Àh=2 E ðzÞ n ðzÞ n ðzÞ 1 À l 2 dz Àh=2 ðz À z 0 ÞE ðzÞ n ðzÞ 1 À l 2 dz Equation 45 shows that the variation of the bending component of the lateral displacement (dw b ) is coupled with the variations of the in-plane displacements (du, dv) by the terms multiplied by the integral B. In addition, the variation of the shear component of the lateral displacement (dw s ) is coupled with the variations of the in-plane displacements (du, dv) by the terms multiplied by the integral Q 4 . Equating B to zero makes the reference plane to be at the physical neutral plane where no coupling between the bending and the in-plane displacements. However, in order to obtain a set of equations with uncoupled lateral and in-plane displacements, Q 4 has also to be zero. Equating Q 4 to zero necessitates at least one of the following two statements. The first, n ðzÞ is zero. The second, E ðzÞ is constant and the integration of n ðzÞ along the thickness is zero. if n ðzÞ is zero, the RPT reduces to the classical plate theory (CPT). Otherwise, if the E ðzÞ is constant, then the material is homogeneous. Therefore, for buckling analysis of FGM plates, using the RPT will always involve the solution for all the 4 variables (dw b ,dw s ,du,dv), even if the physical neutral plane is selected as the reference plane. In the case of homogeneous plates, the RPT provides much simpler buckling equations consisting only of the lateral displacement components (dw b ,dw s ).

Notes on the resulting formulations
In the buckling equations based on FSDT, the lateral and in-plane displacements/rotations are found uncoupled at the physical neutral plane. This is also found true for the case of the refined-FSDT. In the case of the RPT, the lateral and in-plane displacements are found coupled, even at the physical neutral plane. So, the buckling equations of FGM plates based on the RPT contain all the four variables (dw b , dw s , du, dv) coupled together, while only the lateral displacements/ rotations are contained in the stability equations of the FSDT and the refined-FSDT. The buckling equations based on the refined-FSDT, Eq. 40, consist of only two variables (dw b ,dw s ), while those based on the FSDT, Eq. (28 and 34), consist of three variables (dw,dh,du). Although the number of variables is less in the refined-FSDT stability equations, they are not necessarily simpler than those of the FSDT. That because the former contain forth-order terms, while the latter contain terms with order no higher than two. Actually, the two sets of governing equations are similar in terms of simplicity, since the equivalent first-order system of each has six equations. However, the refined-FSDT is more accurate compared to the FSDT [22]. That results as the refined-FSDT correctly describes the transverse shear strain using an imposed shear correction function g ðzÞ , while the FSDT always assumes the transverse shear strain as a parabolic distribution and compensates for the error using shear correction factor k s .
Another point worth noticing is that the two variables of the stability equations of the refined-FSDT are coupled only at the resultant forces terms. The traditional refined theories come last in this comparison of simplicity.
The main advantage of refined-FSDT over FSDT is the elimination of the need for a shear correction factor k s . Instead, it requires a shear distribution function g ðzÞ , which has to be selected wisely to satisfy the traction conditions at the top and bottom surfaces of the plate, and dictate the desired distribution of the transverse shear strains through the thickness. In the other hand, the shear correction factor k s required by FSDT is a simple numeric value. For homogeneous isotropic plates, the often accepted shear correction factors are (5=6) and (p 2 =12). In such a case, FSDT formulations provide identical results to those of higher-order theories, even for thicker plates. However, these correction factors are not necessarily accurate in the case of the FGM plates. Anyway, using of those correction factors produces marginal error, unless for the case of thicker FGM plates with side-to-thickness ratio a=h less than 10. However, many models are proposed to obtain modified correction factors for FGM, e.g. [3,20,21,24]. In short, the EKM buckling formulations based on the FSDT and the refined-FSDT both are simple and practical to use.

Strategy of solution
As mentioned in Sect. 1.1, EKM reduces the multidimensional eigenvalue PDE into separated onedimensional eigenvalue ODE's that are to be solved iteratively starting by arbitrary assumed solution. Figure 2 shows the EKM algorithm for solving the buckling problem of FGM plate based on FSDT that was derived in Sect. 2. The algorithm starts by assuming an arbitrary solution in the y direction and then solves for the buckling factor k x and solution functions in the x direction. Then, considering the lately found solution in the x direction, the one in the y is solved for. These steps are repeated over and over till no considerable change is obtainable from proceeding the iterations.

Conclusion
The buckling equations and the boundary conditions of rectangular plates are derived in the frame of the EKM. The formulations are based on the FSDT, the refined-FSDT, and the RPT. The first presented derivation, that is based on the FSDT, is claimed to be the correct one versus the wrong one found in the literature. The second and third derivations, those are based on the refined-FSDT and RPT, are the first implementations of EKM to a refined plate theory. It is also the first time to implement the refined-FSDT in buckling analysis. It was concluded that the EKM buckling formulations based on the FSDT and refined-FSDT are equally simple in terms of equivalent firstorder number of equations. The buckling formulations of FGM plates based on the RPT are found coupling the lateral and in-plane displacements even at the neutral plane.
The von Ka´rma´n geometrical nonlinearity is considered. General resultant in-plane forces, thickness variation, and material properties variation through the thickness are considered. The physical neutral plane is taken as the reference plane. The stability equations are derived according to the Trefftz criterion, i.e., the principle of minimum potential energy difference, using the variational calculus. Implementing Trefftz criterion is a valid way to derive the stability equations, opposite to the equilibrium condition that has been wrongly implemented in many other previous publications to derive the stability equations of plates using EKM. This presented article exposes the latter fundamentally wrong approach. An important note is that if wrong notations and misinterpretation were to be tolerated, then the wrongly derived stability equations in the literature would resemble the correct stability equation.
EKM converts the buckling equations and boundary conditions from being 2D eigenvalue PDE into two sets of 1D eigenvalue ordinary differential set of equations (ODE), which are much easier to solve. Additional aspects of the plate problem, like the skew angle, and the effect of the elastic foundation, can easily be added to the derivation. The derivation is provided in detail so fellow researchers can easily redo the derivation for other plate configurations of their concern. By following the same exact steps illustrated here, the stability equations and boundary conditions of the plate based on other plate theories can also be obtained.
In future works, one can compare in detail the presented formulations in terms of accuracy, speed, and convergence. In addition, the effect of the shear correction function n ðzÞ on the accuracy of the refined-FSDT should also be investigated. Future works would also present the derivation of the stability equations based on other stability criteria, e.g. the neutral equilibrium method. A comparison might be conducted between the different derivations in terms of simplicity.