Elastic-Plastic Buckling and Postbuckling Finite Element Analysis of Plates Using Higher-Order Theory

In this paper, some of the displacement-based plate theories are used to investigate the elasticplastic analysis of plates in the framework of the ̄nite element method including the buckling and postbuckling e®ects with the focus on the general third-order plate theory (GTPT). The plate calculation results were compared with the results obtained using 64-nodes solid elements involving Lobatto integration scheme. The problem is solved using the Newton–Raphson method applying modi ̄ed Cris ̄eld constant arc-length procedure. The results show good agreement of results and the GTPT can be acknowledged to ful ̄ll essential criteria for application to the elastic-plastic analysis of thin and thick plates.


Introduction
Many structural elements in various branches of technology can be modeled using a plate model. Buckling is one of the possible collapse modes of slender structural elements and therefore much attention is paid in the literature to the buckling of plates. Important thing to know in order to avoid the structural failure is of course the buckling stress. However, in some cases, it is also useful to follow the postbuckling structural response. Typically, we deal with elastic buckling although for large loading and speci¯c proportions of stockier structural elements buckling may occur in plastic range. Behavior of low-carbon steels can be modeled with a bilinear constitutive law and explicit yield stress. However, the stainless steels and aluminum alloys exhibit nonlinear behavior with the plastic strains appearing even for relatively small stresses.
Many papers focused on comparison of incremental and deformation theories. Shrivastava 1 analyzed the inelastic buckling by including the e®ects of transverse shear by deformation and incremental theories of plasticity. Owen and Figueiras 2 used a thick shell formulation accounting for shear deformation based on a degenerate three-dimensional continuum element. Watanabe and Kondo 3 formulated the elastic-plastic incremental tangent sti®ness matrix without numerical integration over the area of the element deriving increments of the nodal displacements from the yield condition. Ore and Durban 4 presented a linear buckling analysis of annular elasto-plastic plates under shear loads using both plasticity theories. Papadopoulos and Taylor 5 derived the discrete¯eld equations from a nonlinear version of the Hu-Washizu variational principle and applied them for a¯nite element analysis of elastoplastic Reissner-Mindlin plates. Tugcu 6 analyzed the buckling of elastic-plastic plates subject to biaxial loading using a bifurcation-of-equilibrium approach. Durban and Zuckerman 7 carried out the analysis of rectangular plates under uniaxial loading. The problem of in°uence of the plate slenderness on the buckling type (elastic versus elastic-plastic) was addressed by Betten and Shin. 8 Wang et al. 9 derived analytical elastic/plastic stability criteria for the case of the elastic-plastic buckling of thin and thick plates employing the constitutive equations given by Chakrabarty. 10 Wang and Aung 11 investigated the buckling of thick plates using the Ritz method. The Ritz method was also applied by Smith et al. 12 who investigated the plastic buckling of steel plates applying the classical plate theory (CPT). The paper by El-Sawy et al. 13 is an example of application of the¯nite element method to the analysis of the elastic-plastic buckling. Wang and Huang 14 and Zhang and Wang 15 employed the di®erential quadrature method. This approach was used also by Hasrati et al. 16 for postbuckling analysis of compressed plates. Kadkhodayan and Maarefdoust analyzed the elastic/plastic buckling of thin rectangular 17 and skew 18 plates under various loads and boundary conditions using Generalized Di®erential Quadrature discretization technique. Ruocco 19 solved the di®erential equations de¯ning the buckling problem of sti®ened plates analytically using the Kantorovich technique.
Various deformation theories have been developed for plates. However, in reference to the elastic-plastic analysis, mostly the CPT and the¯rst-order shear deformation theory (FSDT) were used, in the framework of either analytical or numerical approach.
The drawbacks of the CPT and FSDT are well known and have been thoroughly discussed in the literature. These problems can be overcome applying higher-order shear deformation plate theories (HOPT) which o®er accurate solution and allow to avoid problems related to¯rst-order theories. The HOPT use higher-order polynomials in the expansion of the displacement components through the thickness of the plate. According to the assumptions of HOPT, the restriction on warping of the cross section is relaxed and allows for variation in the thickness direction of the plate. Unlike the FSDT, the HOPT require no shear correction factors.
The formulations in the framework of the HOPT can be categorized with respect to various aspects such as the form of expansion of displacements, particularly throughout the thickness, precise ful¯lling shear stress free conditions at the top and bottom surfaces and accounting for extensibility of the transverse normal to the mid-surface.
Many of the theories neglect transverse normal strains. [20][21][22][23][24][25][26][27] One of the¯rst formulations which included extensibility of the transverse normal to the mid-surface, signi¯cant in the analysis of thick plates and laminated plates with laminae having widely di®erent elastic constants, was proposed by Reissner. 28 Lo et al. 29 derived a theory of plates accounting for the e®ects of transverse shear deformation, transverse normal strain, and a nonlinear distribution of the in-plane displacements with respect to the thickness coordinate. Reissner also combined the distribution of shear stresses through thickness and extensibility of the normal to the midsurface. 30 Other examples of theories including the extensibility can also be found in Refs. 31-35. Reddy developed a third-order theory with the displacement¯eld capable to describe the quadratic variation of transverse shear strains and stresses and vanishing of transverse shear stresses on the top and bottom surfaces of a plate. 20,36 Leung et al. 37 proposed a third-order plate theory releasing the shear-free condition on the top and bottom surfaces. The theory was also applied by Saidi et al. 38 to axisymmetric bending and buckling of functionally graded solid circular plates. Jemielita 39 reviewed the governing equations of the re¯ned theories of plates including those developed by Reissner, 40 Schmidt, 41 Levinson, 42 Murthy 43 and Kączkowski 44 including also his own approach. 45 Other examples can be found in Refs. 46-49. Reddy and Kim 50 formulated a general third-order theory for the plate deformations (GTPT) and demonstrated that it contains the existing plate theories: the general third-order theory with tangential traction free surfaces, the Reddy thirdorder theory, the¯rst-order plate theory and the CPT. Later they applied the GTPT for the¯nite element analysis of functionally graded plates with a modi¯ed couple stress e®ect and the von K arm an nonlinearity. 51 The¯rst to apply the HOPT for the elastic-plastic analysis of bending of plates were Kant et al. 52 In their study, they used an HOPT model which included distortion of the transverse normal and used it for the elasto-plastic analysis of plate bending, employing the von Mises criterion, associated°ow rule and isotropic hardening. The nonlinear equations were solved using the modi¯ed Newton-Raphson method. In their development, they also employed the re¯ned higher-order theory formulated by Kant 53 and Kant et al. 54 which was also applied by Pandya and Kant to orthotropic plates. 55 In this paper, some of the displacement-based plate theories are applied to investigate the elastic-plastic analysis of plates in the framework of the¯nite element method including the buckling and postbuckling e®ects with the focus on the GTPT. The results are compared to those available in the literature and the authors' own three-dimensional elastic and elastic-plastic analyses using solid¯nite elements.

First-order shear deformation theory
We begin the presentation with the FSDT based on the displacement¯eld uðx; y; zÞ ¼ u m ðx; yÞ þ z x ðx; yÞ; vðx; y; zÞ ¼ v m ðx; yÞ þ z y ðx; yÞ; wðx; y; zÞ ¼ w m ðx; yÞ; ð1Þ where x and y are coordinates related to the plane where the midsurface is positioned, z is the coordinate perpendicular to x and y, u m , v m are in-plane displacements, w m is an out-of-plane displacement, x and y are rotations of a transverse line and u, v and w are displacement of an arbitrary point. Equation (1) shows that the following Kirchho® assumptions are valid: the in-plane displacements vary linearly through the plate thickness and the transverse displacement (in the direction perpendicular to the mid-surface) is constant through the thickness which is equivalent to inextensibility of the normal to the midsurface and to vanishing of the normal stresses zz thus, the constitutive relationships are based on the plane stress problem. Regarding the transverse shear stresses the consequence of assuming the displacements according to Eq. (1) is that the transverse normal lines do not remain orthogonal to the midsurface after deformation, thus allowing for transverse shear strains and stresses. The stresses are, however, uniform throughout the thickness which is contrary to the distribution of the transverse shear stresses obtained from 3D elasticity theory which is quadratic and the stresses vanish at the upper and lower plate surfaces. The in°uence of the uniform shear stress is compensated by introducing a shear correction parameter modifying the shear constitutive relationship.
Since the displacements u, v and w are dependent on the variables u m , v m , w m , x and y application of the theory in the framework of the¯nite element method calls for C 0 continuity.
Using the displacement¯eld given by Eq. (1), the strains including the nonlinear von K arm an strain components can be derived as follows: M. Taczała et al.

General third-order plate theory
The third-order plate deformation theories with vanishing surface tractions such as the general third-order theory with tangential traction free surfaces or the Reddy third-order theory have certain drawbacks related to the fact that their formulations require the C 0 -interpolation for u m ; v m ; x ; y and the Hermite interpolation for w m ; z ; ' z . Moreover, in some cases, these theories result in unsymmetrical¯nite element sti®ness matrices even for a linear case. 50 The problem related to various interpolations was addressed by Pandya and Kant 56 who proposed a method of developing an isoparametric displacement¯nite element formulation including the conditions for vanishing of the transverse shear partly during de¯ning the dis-placement¯eld as well as when formulating the shear rigidity matrix and used it for the laminated composite plates. We also note that the condition of vanishing shear stresses at the top and bottom surfaces expressed in an easy manner via shear strains in the elastic range, in the elasto-plasticity range could become more complex due to greater number of nonzero terms in the constitutive matrix.
Reddy and Kim 50 proposed the formulation free from the described limitations. The displacement¯eld for the GTPT is as follows: uðx; y; zÞ ¼ u m ðx; yÞ þ z x ðx; yÞ þ z 2 ' x ðx; yÞ þ z 3 x ðx; yÞ; vðx; y; zÞ ¼ v m ðx; yÞ þ z y ðx; yÞ þ z 2 ' y ðx; yÞ þ z 3 y ðx; yÞ; wðx; y; zÞ ¼ w m ðx; yÞ þ z z ðx; yÞ þ z 2 ' z ðx; yÞ: Assuming the in-plane displacements u, v in the form of the cubic polynomial and the out-of-plane displacement w in the quadratic polynomial with respect to z we obtain a quadratic variation of the transverse shear in this direction with all the displacements contributing to this distribution. The formulation was employed to derive the equations of motion with the use of the modi¯ed couple stress theory for functionally graded material (FGM) plates. The same formulation was also presented 51 for the analysis of the bending de°ections of FGM plates. In both cases, the von K arm an nonlinear strains were considered. Similar approach was proposed also by the other authors. [57][58][59] In Eq. (3), we have eleven generalized displacements: displacements at the midsurface u m ; v m ; w m , rotations of the transverse normal x ; y as well as higher order displacements which have more complex physical interpretation z ; ' x ; ' y ; ' z ; x ; y . For instance, z is a constant term in the expression for strain " z (and the total strain at the mid-surface): where ' z is a multiplier of the linear term of the strain variation The assumed displacement¯eld given by Eq. (3) allows for the parabolic variation of transverse shear strains. The cubic variation of in-plane displacements causes the transverse normal to deteriorate from the straight form while the quadratic variation of out-of-plane displacement implies extension through the thickness thus leading to varying thickness of the plate and emerging direct stresses in the direction of z-coordinate. The von K arm an nonlinear strain-displacement relations obtained by assuming small strains and moderately large rotations can be approximated as follows: @w @y 2 % @w m @y 2 ; @w @x @w @y % @w m @x @w m @y ; ð6Þ with all the remaining quadratic terms equal to 0, leading to the formulation as follows: We note that applying the shear-free condition becomes straightforward in the case of adopting these assumptions. Developing the formulation toward more complex geometry of the model such as sti®ened plates would call for including all nonlinear terms as they can play an important role when the°exural-torsional or tripping buckling modes occur. 60,61 Considering all relevant aspects we selected the general third-order shear deformation theory described by Eq. (3) to be used to in the elastic-plastic analysis of plates. Kim and Reddy 51 developed a¯nite element model of the GTPT which requires the use of C 1 -continuity. However, the model was used to the modi¯ed couple stress theory accounting for the microstructural e®ect through a single length scale parameter involving the second derivatives of the variables. That problem is not covered in the present analysis and thus adopting the C 0 -approximation for all 11 generalized displacements appears to be an appropriate approach. 50

Three-dimensional analysis
While isotropic as well as laminated composites and FGM plates are typically analyzed using plate elements, for some cases it can be necessary to apply threedimensional solid elements with one or several solid elements through the thickness. The use of such elements or 3D analytical solutions may be particularly advantageous when transverse shear e®ects are predominant and the normal transverse (through thickness) stress should be accounted for. Applying the solid elements, we obtain results which can be used as reference values to verify accuracy of the results obtained with the two-dimensional idealization with using a plate theory. As an example Werner 62 presented a Navier-type three-dimensionally exact solution for bending of linear elastic isotropic rectangular plates. Also Demasi 63 developed a Navier-type method for¯nding the exact three-dimensional solution for isotropic thick and thin rectangular plates. The formulation was next extended to the case of multilayered plates composed of isotropic layers. 64 Zenkour 65,66 applied the 3D analytical elasticity solutions as a reference for functionally graded plates. A comparison of results according to the Mindlin and Reddy plate models with three-dimensional solutions for elastic simply supported isotropic and transverse inextensible rectangular plate can be found in Ref. 67.
The 3D analysis allows for the application of a fully three-dimensional elastic and elastic-plastic constitutive equations, similarly as in the GTPT. Thus, the transverse direct and shear stresses in this case are calculated directly without assumptions regarding the kinematics of deformations which provides a distribution of stresses close to reality.
In the present analysis, the nonlinear strain-displacement equations involving all nonlinear terms typically used for the geometrically nonlinear analysis with the application of a solid element are replaced by the von K arm an equations to retain consistency with the formulations for the plate model: " yy ¼ @v @y þ 1 2 @w @y 2 ; " zz ¼ @w @z ; xy ¼ @u @y þ @v @x þ @w @x @w @y ; Elastic-Plastic Buckling and Postbuckling Finite Element Analysis

Derivation of equilibrium equation
Nonlinear¯nite element equations in the incremental formulation are derived using the principle of virtual work which for increment t þ Át and iteration i þ 1 is given by where tþÁt iþ1 W ext is the virtual work of external forces for an increment t þ Át and where b is the vector of the body forces acting in volume V , p is the loading distributed over area , while u denotes the displacement functions dependent on the formulation corresponding to either the plate theory or the solid formulation. Virtual work of internal forces (second Piola-Kirchho® stresses) tþÁt iþ1 W int has the form as follows: Increments of the Green-Lagrange strains tþÁt iþ1 Á" assuming large displacements can be derived using von K arm an nonlinear strain-displacement relations.
Applying the¯nite element approximation and introducing strain-displacement matrices, B ð1Þ and B ð2Þ , the increments of the strains, are derived from the strains de¯ned by Eqs. (2), (7), and (8).
where d are nodal displacements. The strain increments and matrices B ð1Þ and B ð2Þ are given in appendix. The stress increments are evaluated using the constitutive relationship as follows: The principle of virtual work using Eqs. (11)- (14) takes the form as follows: is the sti®ness matrix dependent on displacements (including also the linear term) is the sti®ness matrix dependent on stresses, is the reference load vector, and is the internal force vector. The buckling stress (bifurcation point) is found from the condition The structural response in the postbuckling regime was analyzed applying the pathfollowing technique in the form of the Cris¯eld constant arc-length method by adopting a constraint condition in addition to the equation set. The constraint delimiting the displacement increment in each load step tþÁt iþ1 Ád incr is expressed by

Finite elements for plate models
In the present analyses, the FSDT and GTPT were employed. The¯nite element used for the plate is a sixteen-noded isoparametric Mindlin plate¯nite element of the Lagrange family with Lobatto numerical integration, which is free from shear locking 68 and thus also adequate for the analysis thin plates. Shape functions for the plate element are as follows: where N ðLÞ I ðÞ are the shape functions derived using the one-dimensional Lagrange polynomials, and taking the location of the integration points according to the Lobatto quadrature 1 ¼ À1 More details are given in Ref. 68 where the results of the patch tests are presented and it is shown that the element is free from shear locking thus being adequate for the analysis of thin plates.

Finite elements for solid model
The 27-noded and 64-noded hexahedral elements are used to set up models for the analysis of bending and buckling of plates. To achieve appropriate accuracy, it is necessary to keep the element aspect ratio as close as possible to a regular cubic form. Moreover, the number of elements in the thickness direction is important for the solution accuracy in thick plates. The 27-noded element was discussed in Ref. 69 while the 64-noded solid element was developed speci¯cally here for the purpose of comparative analysis with the third-order plate theory as an extension of the 16noded plate element in the third direction. The Lobatto integration technique was also employed for this element as the position of the Lobatto points allows to investigate the value of shear stress on the delimiting surfaces. The shape functions are a natural expansion of the functions de¯ned in Eqs. (21) and (22) are as follows: The formulation allows to investigate the stresses at the boundaries (in this case the delimiting top and bottom surfaces) to speci¯cally verify the value of the shear stress there.

Material Model
The material behavior is modeled by the Ramberg-Osgood elastoplastic stress-strain relationship assumed in the form where the¯rst term corresponds to the e®ective elastic strain " E e , while the second to the e®ective plastic strain " P e , e is the combined stress, Y is the nominal yield stress while k and n are nondimensional parameters de¯ning the shape of the stress-strain elastoplastic curve.
From Eq. (24), we can evaluate the strain hardening parameter as follows: We note that for small values of e the strain hardening has very large value with the in¯nity for e ¼ 0. When performing incremental analysis using the¯nite element method with no-initial-stress condition this singularity can be overcome by assuming the elastic relationship in the initial phase. In fact, in quite large range of stresses the di®erence of behavior described by the linear equation (the¯rst term of Eq. (22)) and the Ramberg-Osgood model is negligible.
In reference to elastic-plastic response the application of the¯rst-order theory is described in Refs. 9 and 15, following Ref. 10. From equations presented there it follows that the relationships given there are of purely elastic type for the shear stress and strains.
Using an appropriate higher-order deformation theory we can use the kinematic assumption allowing for distribution of shear stresses (true shear stress) and employ constitutive equations resulting from the plastic°ow theory. Nevertheless, the results for the¯rst-order shear theory will also be presented.
We note that the constitutive matrix in Eq. (13) is di®erent depending on the plate theory. For the GTPT where there are six stress and strain components it is formed directly from the equations de¯ning the Prandtl-Reuss theory of plasticity while for FSDT the condition of inextensibility in the thickness direction must be used as follows: which results in a transformation of the constitutive matrix to the form where For GTPT, where six stress and strain components are involved, no modi¯cation takes place so that

Veri¯cation of the formulationelastic bending
The GTPT will be veri¯ed comparing the results obtained using this theory with the results obtained using solid elements. Since many results referring to the elasticplastic behavior plates are obtained using the FSDT, the presented results will also be compared to this theory. Obviously, the principal di®erence between the FSDT which requires the corrective factor to account for shear stresses and higher order theories which do not require this factor is a possibility to include the shear stresses directly in the constitutive equations upon the condition that the distribution of these stresses correspond to real distribution. Here, it is believed that the \real" distribution of the stresses is provided by the models based on solid elements. Sixty-four-noded solid elements with the Gauss-Lobatto integration scheme as well as twenty-seven-noded elements are used for the 3D analysis. A simply supported square plate subject to the uniformly distributed loading over the plate area is investigated¯rst. Dimensions of the plate are b ¼ 800 mm, t ¼ 80 mm, thus b=t ¼ 10 which means a thick plate. The Young modulus of the plate material was assumed equal to 73776.5 N/mm 2 while Poisson ratio 0.32. The loading was 3:90625 Â 10 À3 N/mm 3 for the solid model and the equivalent value of 0.3125 N/mm 2

2150095-11
for the plate model. Both values are the reference values and the actual loading is the reference value multiplied by the loading factor. Meshing was 16 Â 16 for the plate models and 16 Â 16 Â 4 for the solid models, thus keeping the aspect ratio within reasonable limits. A quarter of the plate was actually modeled (Fig. 1) with the appropriate boundary conditions representing symmetry. We note that for the GTPT these are not only the in-plane displacements and rotations but also the second and third order variables. Regarding the simply supported edges, the kinematic conditions are speci¯ed as follows: wð0; yÞ ¼ 0; wðb; yÞ ¼ 0; wðx; 0Þ ¼ 0; wðx; bÞ ¼ 0; Bending moment along edges y= 0 and y ¼ b is evaluated as follows: while the normal stress can be found using constitutive relationship

2150095-12
Using the strains components resulting from the displacement¯eld (Eq. (2)), we obtain the following condition: Thus, to ful¯ll the condition expressed by Eq. (31) we implement for the edges y = 0 and y ¼ b: Similarly, we apply for the edges x= 0 and x ¼ b Analogous boundary conditions should be applied in the case of the solid element.
Since rotations do not appear in the 3D formulation, we insert the boundary conditions v ¼ 0 at x ¼ 0 and x ¼ a as well as u ¼ 0 at y ¼ 0 and y ¼ b.
As a next step, the loading is changed to cover only a 50 Â 50 mm part of the plate in the center of the plate (patch loading), see Fig. 2. The resultant force is identical as in the previous case which means that the value of loading is now 250 Â 10 À3 N/ mm 3 for the solid model and the equivalent value of 20 N/mm 2 for the plate model (again both being the reference values). This type of loading induces considerable shear forces and the purpose is to investigate the response of both plate theories to verify the GTPT via comparison with the solid model in the presence of shear forces more signi¯cant than in the case of the uniform loading.
Presented results concern nonlinear analysis, both geometrical and material, as the main purpose was the veri¯cation of the GTPT for this type of analysis. The results are given for loading factors ¼ 1:0, ¼ 40:0 and ¼ 55:0. The¯rst loading Fig. 2. Simply supported square plate subject to patch loading.

2150095-13
level slightly di®ers from elastic response with small stresses resulting in positioning the response in the initial part of the Ramberg-Osgood curve and with the de°ections not contributing signi¯cantly to the geometrically nonlinear response. Two latter loading multipliers are large enough to generate plastic strains as well as more signi¯cant geometric e®ects.
The results in terms of de°ections at the center of the plate are given in Fig. 3. We can observe a good agreement of results with the di®erence between SOLID 64 and GTPT at the level of 0.3%. We can also see a good agreement of results between the two solid models which also con¯rms correct formulation of the SOLID 64 element proposed in this paper.
In the case of the plate subject to the patch loading, the di®erence between the de°ections is also small, and even if greater than in the previous case it is still less than 1% - Fig. 4.

2150095-14
We also note that the de°ections obtained using the FSDT are close to the values obtained using the other formulation -0.1558 mm for the uniform loading and 0.1187 mm for the patch loading.
We also compare distributions of the in-plane displacements across the thickness for coordinates X ¼ Y ¼ 375 mm with the displacement in the center obviously equal to 0. We see the excellent agreement of the displacements for both uniform loading (Fig. 5) and patch loading (Fig. 6). In the case of the patch loading, we observe a slightly greater value of the displacement close to the outer surfaces for SOLID 27. However, the agreement of GTPT and SOLID 64 is very good.
Comparing the response of the plates subject to uniform and patch loading we also present the distributions of maximal normal stresses (for outer surface) along the diagonal of the plate. The results are presented in Figs. 7 and 8 in terms of stresses versus the X coordinate.

2150095-15
We observe that the stresses are practically identical in the case of the uniform loading. However, there is a remarkable di®erence of stresses between the FSDT and GTPT. The GTPT model better follows the stress distribution obtained using the solid model with the di®erence of maximal stress 1.5% while the FSDT di®ers by 15%.
We also present the distribution of the maximal midsurface shear stresses along the diagonal, comparing the GTPT to the SOLID 64 model. We can see a good representation of the shear stresses even though in the location of the largest stressat the border of the patchthe di®erence is greater, Fig. 9.

2150095-16
Distribution of direct stresses for simply supported square plate subject to uniform loading across thickness evaluated at plate center as well as distribution of shear stresses for simply supported square plate subject to uniform loading across thickness at the location of large shear strains (X ¼ 400 mm, Y ¼ 50 mm) are presented in Figs. 10 and 11, respectively. We can see the linear distribution of the direct stress corresponding to the elastic response. As it was already mentioned, the load factor ¼ 1:0 despite involving both types of nonlinearity corresponds in fact to the linear response.
An important test of a correct behavior of the plate model is the result of zero shear stresses at the delimiting surfaces. In the GTPT, there is no condition enforcing  Distributions of stresses for the patch loading are presented in Figs. 12 and 13 for direct and shear stress, respectively. The stresses are present for the same locations as in the case of the uniform loading. Similarly to the results referring to displacements, the distributions di®er only slightly. For the normal stress, the agreement is very good, for the shear stress, we observe 10% di®erence of the maximum values, yet the distribution is parabolic. It should be pointed out, however, that the results are given Fig. 11. Distribution of shear stresses for simply supported square plate subject to uniform loading across thickness evaluated at X ¼ 400 mm, Y ¼ 50 mm, ¼ 1:0.

2150095-18
for the locations with the very large stress gradient and in the remaining part of the plate the di®erence is de¯nitely less signi¯cant.
The conclusion from the elastic analysis is that the GTPT is a good candidate for a theory describing the behavior of thick plates. The elastic-plastic tests follow in the next chapter.

Veri¯cation of the formulationelastic-plastic bending
The geometry of the models and the loading applied for the elastic-plastic analysis are the same as used in the previous chapter. Also the Young modulus (E ¼ 73; 776:5 N/mm 2 ) and Poisson ratio ( ¼ 0:32) were identical. Nonlinear material behavior was modeled using the Ramberg-Osgood model with nominal yield stress 0 ¼ 423:353 N/mm 2 , n ¼ 20 and k ¼ 0:3485 what corresponds to the 14S-T6 aluminum alloy.
The results are presented for loading factors ¼ 55:0 (uniform loading) and ¼ 55:0 (patch loading). The di®erence is due to the fact that, even though the resultant forces are identical, the patch loading induces greater stresses. We begin with the presentation of distribution of de°ections for the plate center, Figs. 14 and 15.
We see that in the case of elastic-plastic response, the de°ections calculated using the GTPT are in a good agreement with the SOLID 64 results. The di®erence is 0.5% for the uniform loading, again in the case of the patch loading is greater. However, it does not exceed 1.5%.
From results presented in Figs. [16][17][18][19] where the distributions of direct and shear stresses are given for uniform and patch loading we can con¯rm acceptable agreement of results for both types of loading. Regarding the normal stress we observe very good agreement between the SOLID 64 and GTPT models. For the shear stress, Fig. 13. Distribution of shear stresses for simply supported square plate subject to patch loading across thickness evaluated at X ¼ 400 mm, Y ¼ 325 mm, ¼ 1:0.

2150095-19
we note very good behavior the of SOLID 64 model yielding the values practically zero at these surfaces with the acceptable values for the GTPT.

Veri¯cation of formulationbuckling
To verify the correctness of the formulations, we analyze the plate originally investigated by Shrivastava, 1 with the results given also by Wang et al. 9 We analyzed square plate subject to uniaxial compressive load. The material of the plate was modeled using the Ramberg-Osgood model with Young modulus equal to 73,776.5 N/mm 2 (originally 10,700 ksi), Poisson ratio ¼ 0:32, nominal yield stress  The results of the buckling stress, calculated for the FSDT and GTPT plate theories and the two solid elements: 27-noded and 64-noded, are presented in Table 1 in comparison to those presented by Wang et al. 9 according to the incremental theory.
It should be explained that in the case of the solid elements the boundary conditions described in Sec. 5.1. are automatically ful¯lled as the displacements at the side surfaces are identical for each y coordinate when the equivalent loading vector is applied evaluated from the pressure acting on the compressed surface of the model.

2150095-22
We can see a good agreement of results between the GTPT and the reference values of the buckling stress, while for the CPT and FSDT, in the case of thicker plates ðb=t < 25Þ the buckling stress is lower and the course of the curve representing relationship buckling stress versus slenderness is di®erent (Fig. 12). We also note slightly less buckling stress in the case of solid model, with the di®erence ranging from 1.3% to 3.7% comparing to the GTPT but we¯nd the good agreement between the buckling stress for both solid models.

Numerical examples
In order to investigate the in°uence of boundary conditions, type of loading as well as material characteristics we shall consider several numerical cases involving simply supported and clamped edges with respect to bending as well as restraining e®ect (blocking midsurface along an edge) and with respect to in-plane displacements, uniaxial and biaxial loading and various stress-strain curves, all in the framework of the Ramberg-Osgood model. We shall also present the nonlinear response of compressed plates in the postbuckling range for various square plate slenderness.
The material of the plate is essentially the same as used in the previous analyses, the Young modulus being equal to E ¼ 73; 776:5 N/mm 2 , Poisson ratio ¼ 0:32, nominal yield stress 0 ¼ 423:353 N/mm 2 , n ¼ 20 and k ¼ 0:3485. Dimensions of the plates are 800 Â 800 mm with thickness depending on the slenderness.
As the¯rst example we consider again the model presented in Sec. 5.3, the difference being that we take stockier plates exhibiting slenderness of b=t ¼ 14 to 20. The comparison of results in terms of the buckling stress is given in Table 2.
The variation of the buckling stress versus slenderness, including also the values from Table 1, is presented in Fig. 20.
We can observe a very good agreement of results for the FSDT and GTPT theories for the plates having slenderness greater than 22 and the di®erence at the level of up to 5% for the stockier plates. Since the di®erence between these two theories in the buckling analysis is the inextensibility of the normal transverse in the FSDT and its extensibility in the GTPT we can formulate the conclusion that this is an essential e®ect for the plates with b=t < 20.
In the next study, the in°uence of the boundary conditions on the buckling stress of the elastic-plastic plate is examined. We continue with the same model as

2150095-24
state of stresses. We note that despite the fact that the response of the plate is stronger, the buckling stresses evaluated using the bifurcation condition are less. The buckling stress decreases with the increase of the slenderness b=t, while for the unrestrained plates the critical stress stabilizes at the value exceeding 600 MPa. We should also note that for the assumed material model the nominal yield stress of 423.3 MPa corresponds to the plastic strain 0.002 while the stress of 470 MPa corresponds to the strain ten times greater so that the buckling stress greater than this value has only a theoretical signi¯cance.
We continue presenting our results considering variation of the material parameters. First we investigate the in°uence of the parameter n on the buckling stress. The Ramberg-Osgood stress-strain curves for values of the parameter ranging from 5 to 40 are shown in Fig. 22. Figure 23 displays the variation of the buckling stress for the uniaxially compressed plates for four values of the coe±cient n (Eq. (24)). They reveal that the buckling stress converges for the plate slenderness greater than 20, where the elastic response dominates e®ects.
Next the in°uence of the nominal yield stress, varying from 80% to 120% of the original value. The stress-strain curves for¯ve various nominal yield stresses are shown in Fig. 24. Here, unlike in the previous example, we can see identical shape of the curves in the initial range and a similar one in the part where the equivalent stress is greater than the nominal yield stress.
Variation of the buckling stress depending on the plate slenderness and the nominal yield stress is presented in Fig. 25. The value of the buckling stress decreases with the increase of the plate slenderness. The di®erence is more apparent for thicker plates where the buckling occurs in the elastic-plastic range while for thinner plates the results converge as the elastic parts of the stress dominate.

2150095-25
As a next example, we present nonlinear behavior of compressed plates in the postbuckling regime. We investigate the response of three plates having slenderness b=t ¼ 23, 26 and 40. The results are presented in Figs. 26-28, respectively in the form of the central de°ection versus averaged stress at the compressed edge.
For the plate having slenderness of 23, the buckling stress exceeds the nominal yield stress, so there is practically no reserve of strength and thus we observe the decreasing loading in the postbuckling range.
For the second case of plate with b=t ¼ 26 for which the buckling stress is lower than for the unrestrained plate we receive a slight increase of compressive stress. We also note di®erent shapes of the stress-strain curves for the two plates, with unrestrained and restrained in-plane boundary conditions, respectively.

2150095-26
Interesting behavior can be observed for the plate with slenderness of b=t ¼ 40. The buckling stress corresponds to the initial part of the Ramberg-Osgood stressstrain curve which corresponds to buckling.
The nonlinear response for plate with b=t ¼ 40 and the restrained in-plane boundary conditions can be explained by observing the de°ection modes at particular stages. First the buckling occurs in the mode typical of a square plate, that is one half-sine wave in each direction, see Fig. 29.

2150095-28
Increasing the loading in the postbuckling range, we observe transition from the elastic buckling mode to another deformation mode. For the maximum loading, we have de°ections as presented in Fig. 30, which should be associated with the ultimate load capacity of the plate.
When we continue the analysis which is possible due to a path-following solution algorithm we observe the change of deformations in the form of a decreasing de-°e ction of the central node, see Fig. 31. The analysis was terminated when the de°ection mode was as presented in Fig. 32. This example also illustrates the di®erence between the response of the plates with restrained and unrestrained boundary conditions, as shown in Fig. 21. The unrestrained plate exhibits greater buckling stress, however, this is a restrained plate which reaches greater ultimate capacity in the postbuckling range.

2150095-29
For the remaining examples referring to stockier plates for which the stress-central de°ection curves are presented in Figs. 26 and 27 as well as to the same plate (b=t ¼ 40) with unrestrained in-plane boundary conditions, the de°ection mode was typical of the buckling mode of a square plate and remained as such throughout the analysis.

Concluding Remarks
In this paper, the application of the third order plate deformation theory originally presented by Reddythe GTPTwas employed for the elastic-plastic analysis of plates using the incremental theory of plasticity with the Prandtl-Reuss constitutive equations and Ramberg-Osgood material model. The approach was veri¯ed in a series of numerical tests in which displacements and stress distributions were compared to the results obtained using solid models. The results show good agreement of results and the GTPT can be acknowledged to ful¯ll essential criteria for the application to the analysis of thick plates. To capture the more practical elastic-plastic behavior, various cases were analyzed with regard to bending of the plates subject to uniform loading and more concentrated loading.
The formulation was also veri¯ed by comparing to the results of buckling analysis for the cases where the buckling stresses were given. Again, a good agreement of results has been obtained.

Acknowledgments
The work has been performed under the project functionally graded materialsstatic and dynamic analyses,¯nanced by the Polish National Science Centre (NCN) under the contract 2018/29/B/ST8/02723. The support is gratefully acknowledged.

Appendix A
The strain increments, evaluated from as the di®erence between iteration i þ 1 and i, are for the FSDT: For the GTPT, the strain increments are as follows: We apply the¯nite element approximation of the displacement functions for the FSDT: ðA:4Þ For the GTPT, we also approximate the remaining displacement functions: