Large deflection and initial instability analysis of anisotropic plates by the generalized finite element method

This paper presents investigations laminated plates under moderately large transverse displacements and initial instability, through the Generalized Finite Element Methods - GFEM. The von Kármán plate hypothesis are used along with Kirchhoff and Reissner-Mindlin kinematic plate bending models to approximate transverse displacements and critical buckling loads. The generalized approximation functions are either C 0 or C k continuous functions, with k being arbitrarily large. It is well known that in GFEM, when both the partition of unity (PoU) and the enrichments functions are polynomials, the stiffness matrices are singular or ill conditioned, which causes additional difficulties in applications that requires the solution of algebraic eigenvalues problems, like in the determination of natural frequencies of vibration or the initial buckling loads. Some investigations regarding this problem are presently addressed and some aspects and advantages of using C k -GFEM are observed. In addition, comparisons are presented between the classical GFEM and the Stable-GFEM (SGFEM) with regard to the evaluation of the initial critical buckling loads. The numerical experiments use reference values from analytical and numerical results obtained in the open literature.


INTRODUCTION
Anisotropic laminated plates are amongst the most commonly used types of structural elements in aeronautical and naval industry.They are often subjected to in-plane compressive and shear loads which, when combined with their slenderness, make then prone to geometric instability by buckling, where various parameters such as thickness, layer stack, loading type, boundary conditions and geometry affect its ability to maintain structural integrity.Plate theories, through some kinematic hypotheses, allow the three-dimensional fields to be approximated by two-dimensional ones, thus reducing the complexity and the number of variables required to estimate its mechanical behavior.However, these simplifications lead to important mathematical consequences, such as: (a) First order kinematic models such as Kirchhoff and Mindlin, within the linear conditions, allow solutions to be exact with regard to the three-dimensional formulation, in the limit as the transverse displacement becomes small compared to the thickness.A traditional rule of thumb of application of these theories is that the errors involved are admissible in engineering applications if the transverse displacements are limited to some fraction of the thickness, e.g., a half.
(b) The simplest inclusion of nonlinear kinematic effects is provided by the von Kármán theory.The resultant strain-deformation relation contain a sub-set of the nonlinear terms present in the full relations.As a consequence, the range of applicability of the theory is usually referred to as displacements moderately large.This means it provides better accuracy in the same displacement range of validity of the linear models, but does not allow arbitrarily large displacements.A second aspect of the von Kármán theory is that it allows a first approximation of initial instability of the panel.This is obtained numerically by a standard linear eigenvalue problem which gives, in a relatively inexpensive way, an estimate for the mode shape and critical load of the structure.Even though the load estimate is only an upper bound (Brush and Almroth, 1975;Ventsel and Krauthammer, 2001), it is often very useful in the early stages of the design process.
With regard to finite element (FE) solutions, the plate bending theories available contain different continuity requirements on the function basis used.For instance, the Mindlin kinematic model requires 0 C continuous approximation functions and the Kirchhoff model requires 1 C functions for the transverse displacements.The Kirchhoff model is interesting by its reduced cost in computational analysis and absence of some of the pathologies which can suffer some of the FE formulations based on the Reissner-Mindlin model, like the shear locking.Frequently, the technologies developed to FE formulations on higher order kinematic models, like the third order shear deformation of Levinson or Reddy (Reddy, 1989(Reddy, , 1984)), are inspired in developments made for the Kirchhoff model by using, for example, non-conforming elements to alleviate the 1 C requirement on the function basis.
As for the 1 C requirement for the Kirchhoff finite element basis, it is virtually unknown in the literature a formulation of arbitrarily shaped plate elements with 1 C continuity everywhere in the mesh.Most of them have to be restricted to particular element shapes, or the continuity is restricted to discrete points, only enough to preclude numerical instability.Only recently a few approaches have been developed to fix this limitation.One of them is a variation of the GFEM (Generalized Finite Element Method), based on smooth k C continuous approximation functions, for arbitrary k (de Barcellos et al., 2009).This is one of the methods which will be tested in this paper, on modeling the nonlinear effects of anisotropic laminated plates by the von Kármán theory with the Kirchhoff and Mindlin kinematic models.
The theoretical and historical basis of the GFEM lies on early developments of meshless methods, particularly on the hp-Clouds method (Duarte andOden, 1996a, 1996b).These methods share some difficulties, for example on integration and on imposition of boundary conditions.Oden et al. (1998) and, in parallel, Babuška (Babuška and Melenk, 1997;Melenk and Babuška, 1996) proposed the use of clouds (patches of elements around each node) as support of approximation basis functions, with the contours conveniently defined by linear finite element meshes, but keeping from the hp-Cloud method its flexible way of enriching the basis with arbitrary functions defined in global coordinates.After this point the GFEM/XFEM developed exponentially, for example, (Belytschko et al., 2009;Belytschko and Black, 1999;Dolbow et al., 2000;Moës et al., 1999;Sukumar et al., 2000).
Most of GFEM, as well as FEM developments, are formulated with 0 C approximation function basis.However, in applications such as the Kirchhoff plate model the required approximation functions must be at least 1 C continuous.
The k C version of GFEM has been developed by the use the Edwards proposal (Edwards, 1996), which allows the construction of arbitrarily smooth weighting functions with support in convex clouds.These weighting functions are used to form a Shepard Partition of Unit (Duarte and Oden, 1995;Shepard, 1968).Duarte et al. (Duarte et al., 2006) extends Edwards proposal to non-convex clouds, by the use of a Boolean R-function.
A known characteristic of the GFEM formulation is that its set of approximating functions is linearly dependent, when it is formed by a polynomial partition of unity uniformly enriched with polynomial functions.As a result, the coefficient matrix (stiffness) contains more null eigenvalues than the number of rigid body motions in the model.When different combinations of types of partition of unity and enrichment functions are used, the rank deficiency can be eliminated but the stiffness matrix remains ill conditioned.Babuška (Babuška andBanerjee, 2012, 2011) proposed a modification in the enrichment functions, in which they are substituted by the difference between the original one and a linear interpolating function, defined in each element.This method, referred to as Stable-GFEM (SGFEM), was shown to improve the approximation space and also to alleviate problems due to non-homogeneous enrichment over the mesh.Gupta et al (2013) compared accuracy and conditioning between GFEM and SGFEM in a two-dimensional fracture mechanics problem, and showed that SGFEM improves the accuracy while keeping the conditioning similar to the one of FEM.Later, Gupta et al (2015) extended and improved the SGFEM formulation to the three-dimensional fracture mechanics.
A practical consequence of the rank deficiency in the stiffness matrix of 0 C -GFEM is that the number n of null eigenvalues is, a-priori, unknown, and can be relatively large when compared with the order N of the matrix.The standard procedure to evaluate the critical load consists in solving an eigenvalue matrix problem and compute the first, smallest, positive eigenvalue.Since N is large, an adequate method is used to compute not all N eigenvalues but only the first few of them.However, the difficulty is posed when one does not know, a-priori, the dimension of the kernel space.This forces the eigenvalue algorithm to work in a vector space of size, at least, 1 n  .On the other hand, the function basis formed by the smooth k C partition of unity uniformly enriched by polynomials is linearly independent (Mendonça et al., 2011).As a result, after imposition of sufficient Dirichlet boundary conditions, the stiffness matrix has a large condition number, but it is not singular.
The objectives of the present paper are the evaluation of the behavior of both types of GFEM, 0 C and k C in deformation response of Mindlin and Kirchhoff kinematic models under moderately large displacements using the von Kármán theory.Also, the response of both stable and non-stable 0 C -GFEM in the evaluation of buckling load is analyzed, with regard to the effects of the stiffness matrix rank on the computational effort to obtain the first useful eigenvalue.
The present paper is organized as follows.Section 2 summarizes a few key aspects of the GFEM, in both forms 0 C and k C , and its stable version; Section 3 shows the basic information about Kirchhoff and Mindlin plate theory with von Kármán hypothesis and its discretization for a generalized finite element procedure; Section 4 details aspects of the initial instability analysis by GFEM and Section 5 presents comparative results for deformation analysis and critical load evaluation in standard problems.

PARTITION OF UNITY AND APPROXIMATION FUNCTIONS
Let us consider a region V, defined by a thickness H and a planar middle surface  with boundary region  .The V region may be described, in Cartesian coordinates, as Let us consider a mesh of elements of domain e  , such that A Partition of Unit (PoU) is defined (Oden and Reddy, 1976) as a set of functions ( ) , such that: 1) ; 2)   is a function at least continuous over the domain, has having its cloud,   , as a compact support and has the required smoothness.Here, the partition of unity functions have at least k continuous derivative inside the cloud, that is, Once a PoU is available, the GFEM approximation functions ( ) The set of enrichment functions can be chosen in different ways.The simplest type is the set of monomials of a given degree.Also, it is very common the use a special function extracted from an asymptotic solution around a critical region of the model.In general, the enrichment functions are C  continuous over an infinite support.The product method Paulo de Tarso Rocha de Mendonça et  The usual GFEM is based on PoU functions which are the traditional FE shape functions, usually the linear tent functions, which naturally constitute a PoU.These functions have low computational cost and are easily integrated by numerical quadrature, although it is limited to 0 C continuity.

Shepard Partition of Unity
For problems requiring higher smoothness of the primal variables, like the transverse displacements in the Kirchhoff and higher order Reddy plate models, the linear tent functions used as PoU in the 0 C -GFEM are not adequate and a different approach is needed to generate approximation basis functions with the required level of continuity.Shepard (1968) The regularity of ( )   x depends only on the regularity of the weighting functions.Therefore, it remains to define the weight functions with the necessary regularity.In order to do that, a coordinate ( ) j  x normal to the straight edge j of the cloud is defined.It is computed as is the coordinate of the mid-point of the edge j , and , j  n is the unity vector normal to the edge pointing to the inside of the cloud.Therefore ( ) j  x is the smallest distance between x and the edge j .Cloud edge functions with continuity C  can be obtained using the following edge function (2.4) The parameter values are chosen as 0.6   and 0.3   according to Duarte et al. (2006).The weighting functions determine the influence of each node within a single element and are defined by where c  is a scale parameter that guarantees 1 W   in the node  x and M  is the number of edges on the boundary of the cloud.
The Shepard PoU functions constructed according to (2.3)-(2.5),as observed by Duarte et al. (2006), have C  differentiability for convex clouds.The regularity of the approximation functions is determined by the regularity of the PoU and the enrichment functions used.Since in the present paper only the polynomial enrichment is tested, the resultant approximation functions basis is C  .If non convex clouds were used, the result would be C  everywhere except at the concave nodes, where it would be k C , after using an appropriate Boolean product.In addition to being able to choose a basis that satisfy the differentiability required by any formulation, it is possible to use different regularities for each unknown field through the combined use of 0 C and k C functions (Mendonça et al., 2013).

Stable-GFEM
The first procedure of stabilization in GFEM was proposed by Babuška (Babuška andBanerjee, 2012, 2011) as a simple local modification in the enrichment functions employed in 0 C -GFEM.They consist in the difference between the original enrichment function and a linear interpolant function defined in each finite element: mod ( ) where mod i L  is the modified SGFEM enrichment function, i L  is the original enrichment function on the cloud   , ( ) is the interpolant function of the enrichment function defined at the element nodes, obtained with piecewise linear FE tent functions.
The approximation functions are obtained by the product of PoU functions with modified enrichment functions in the same way as in (2.2) (2.7)

PLATE MODELS AND GFEM DESCRETIZATION
The kinematic assumptions of Mindlin's model for plates can be summarized as the following, in Cartesian displacement components: and the Kirchhoff hypothesis result in / , where ε is the membrane strain, and κ is the curvature change due to bending.From the von Kármán hypothesis, the membrane deformations ε is divided into linear and non-linear parts 0 ε , NL ε , as The complete set of generalized deformations is the following: The transverse shear deformations s γ naturally results null in the Kirchhoff model.
The plate is considered an anisotropic laminated body, composed by N layers characterized by the following anisotropic linear stress-strain relations (generalized Hooke's law), such that, for each layer l k method Paulo de Tarso Rocha de Mendonça et al.
Latin American Journal of Solids and Structures, 2019, 29(3), e228  6 . The deformations ε are defined as The bilinear operator for the Mindlin model is where the last integral vanishes in the Kirchhoff case.Separating the linear and nonlinear membrane components, the bilinear form becomes where 0 0 and .

Discretization for Kirchhoff and Mindlin plate models
For the Kirchhoff model the generalized displacement fields can be organized as and for Mindlin model as . Each of these components is discretized with the enriched basis of GFEM over an arbitrary element e as where Nne is the number of nodes of element e, ij  are the PoU functions associated to the node j , i represent a component of u and k is the counter associated to the enrichment monomial functions.The limit ij m is the number of enrichment functions of the displacement , where 3 q  or 5 q  for Latin American Journal of Solids and Structures, 2019, 29(3), e228 7/24 Kirchhoff or Mindlin models, respectively, and dofe N is the number of degrees of freedom in the element.Therefore, the discretized generalized displacements are represented in the usual FEM or GFEM way as e e e  u x N x U (3.9) For future use, the transverse displacement is discretized as Next, a vector θ of the small rotations of the straight segments normal to the reference surface is defined along with its variation GU G U , and , where U is the global vector of displacement coefficients and  is the variation operator.The nonlinear strain vector in (3.7) can be decomposed as a product of two linear terms as follows Its first variation can be written in a more convenient form for computational application after a sequence of operations: in (3.14), we have the discretization of the deformation terms and their variations: K U U is the nodal vector of internal forces, the stiffness matrix The equilibrium algebraic system ( ) Following the usual procedure, the Newton-Raphson iterations are based on: where k U is a correction such that in k -th iteration the solution is approximated by U , which can be represented, in indicial form and using summation rule, and 0 ir K is the linear term, given by The first term comes from the differentiation of (3.18), which can be obtained taking The derivative of z A requires the differentiation of w, obtained from (3.10) as Latin American Journal of Solids and Structures, 2019, 29(3), e228 9/24 , 0 , 0 0 0 .0 , , 0 0 0 The symbol The influence of the transverse shear deformations of the Mindlin model appears in the linear stiffness 0 ir K , the second term in the integral in (3.22), which is null in case of Kirchhoff model.

PLATE MODELS AND GFEM DESCRETIZATION
Here we investigate the behavior of the Generalized Finite Element approach to deal with the initial stability problem of the laminated plate associated with the von Kármán's theory.The weak form is the standard one, given by Brush and Almroth (1975).
Here, { , , form a prescribed profile of in-plane resultant stresses acting in the plate, and one searches for the smallest value  which makes the equality holds for arbitrary admissible weighting functions.Hence, the plate is considered initially undeformed and loaded only by the prescribed in-plane stress resultants.In this way,

N N N 
gives an estimate for the initial bucking load of the plate, which is proven to be an upper bound of the real bucking load.After discretization, Eq. (4.1) results in the algebraic eigenvalue problem where 0 K is the linear stiffness matrix,  and U are eigenvalues and eigenvectors and G K is the geometric matrix, real, symmetric and, in general, singular.One searches for the smallest eigenvalue associated to the critical buckling load.Both in FEM and in GFEM, G K has non zero entries only on the positions associated to transverse displacement w , that is, all rows and columns associated with u , v (and x  , and y  in case of the Mindlin model) are zero.In addition, even rows associated with w can be null depending on the prescribed values of x N , y N and xy N .Since only one eigenvalue is needed, the smallest one, an interactive method is usually selected to search only for the first eigenpair, for example the method of subspace iteration or the Lanczos adapted to symmetric matrices.Then, it is sufficient to use a small penalty  on the zero diagonals of G K .These values generate a set of eigenvalues of order 1 /  , which are independent from the first ones, and will be the largest ones in the set.In case of FEM, both K and G K become positive-definite matrices after the imposition of adequate Dirichlet boundary conditions.
In case of the usual C 0 -GFEM, it is proven that the function basis formed by the linear partition of unity uniformly enriched by polynomial functions is linearly dependent.This generates a singular stiffness matrix 0 K .Besides the number of null eigenvalues associated with the number of rigid body motions of the model, there is an a-priori unknown number of null eigenvalues associated with the linear dependency of the basis.Considering the complete set of N eigenvalues of the algebraic system, the first non-null, the one of interest, will be n -th one in the ordered list.Usually, n can be as large as about a half of N .The singularity of 0 K is dealt with in this paper by the use of the , for arbitrary 0 s  ).Now, K becomes positive-definite for C 0 -GFEM.However, the iterative method has to solve for a number of eigenpairs large enough to include the first non-zero one.This defeats the whole purpose of eigenvalue methods like Lanczos or subspace iterations, because the size of the subspace becomes too large to render the method efficient.
On the other hand, function basis formed by the smooth C k partition of unity enriched by polynomial monomials is linearly independent (Mendonça et al., 2011).As a result, after imposition of sufficient Dirichlet boundary conditions, the stiffness matrix still has a large condition number but it is nonsingular.In the present paper the shift technique is used in both C 0 -GFEM and C k -GFEM.

NUMERICAL RESULTS
The present tests with the GFEM formulations were performed modeling the entire plate by a set of uniform meshes, defined by a mesh parameter M, some of which are illustrated in Figure 1.The C k -GFEM behavior in extremely distorted meshes has already been investigated in (Mendonça et al., 2011) and is not considered here.The same can be stated about benchmark problems with singularities, whose GFEM response have been already intensively investigated in the literature.These constraints are imposed by selectively restricting the enrichment coefficients.For example, for a node on a boundary .x const  , the coefficients of functions whose enrichment contains x are not constrained, since these functions are naturally zero on that part of the boundary.In this way, the enriched function is preserved throughout the domain.For arbitrarily shaped domains, where the boundary segments are neither straight nor parallel to a given coordinate axis, the procedure proposed in Garcia et al. ( 2009) can be used, where the partition of unity associated to the boundary node is substituted by it multiplied by a ramp function which is zero at the boundary.
The numerical integration is a special concern when the smooth approximation functions in k C -GFEM is used, because the partition of unity is not polynomial but of rational type.A study about the number of integration points required to integrate adequately the stiffness matrix is described in de Barcellos et al. (2009), using triangular and Gaussian integration rules.In all cases solved in the present paper, an indication of the number of integration points used, nip , is shown, and the triangular rule of Wandzura and Xiao ( 2003) is used.
In comparing results between C 0 -and C k -GFEM, we distinguish the degree p of uniform refinement and the degree b of the least complete polynomial the approximation basis is capable of reproducing.For C 0 -GFEM, the PoU is the linear tent and the enriched basis has a polynomial reproducibility of degree 1 b p   .For C k -GFEM, the largest polynomial the PoU can reproduce is the constant, and as a consequence, b p  .

Homogeneous, isotropic, thin plate problem
The first evaluation of the behavior of 0 C -GFEM and k C -GFEM is the simplest case of a homogeneous isotropic square plate, clamped on all edges, modeled by the Mindlin kinematic model, loaded in the nonlinear range.The problem is interesting for a first evaluation because there is an analytic solution available for the thin plate model by Levy (1942)  .This load generates maximum transverse displacement of about / 1, 7 w h  , which puts the response beyond the usual linear range ( / 0, 5 w h  ), but still inside the range of validity of the von Kármán model.For the results of k C -GFEM in Figure 2, for lower polynomial enrichments 1 p  and 2, the shape of the basis functions and their derivatives are considered simpler and easier to integrate, therefore we use 25 nip  integration points.For higher enrichment degrees, 3 p  and 4, we use 54 npi  points per element, in order to guarantee a proper integration.The results show a tendency for locking with enrichments 0 p  and 1, for 0 C -and k C -GFEM respectively, that is, for a basis of polynomial reproducibility 1 b  .For enrichments 2 p  and higher, 0 C -GFEM results show stabilization with mesh refinements, which can be caused by ill conditioning of the stiffness matrix.For the smooth GFEM, the curves show progressive convergence with h -refinement.The errors obtained when using mesh M1 are almost always high because this mesh is too crude, particularly in Figure 2 k C -GFEM case.
A more systematic evaluation of the effect of number of integration points ( nip ) in the response of k C -GFEM is described in Figure 3. Here, the relative error of transverse displacement at plate center, ,  Figure 4 shows the pointwise normalized transverse displacement w and the normal stress x  versus transverse load P , for 0 C -GFEM, with polynomial uniform enrichment of degree 2 p  and 25 nip  integration points.The displacements are at the plate center and the reference solution is taken from Reddy (Reddy, 2004).This reference uses a finite element solution and the stress is obtained at the point ( , , ) (6.25; 6.25;  / 2) x y z h   , chosen to coincide with one integration point of his model.That model consists of an uniform mesh of 16 nonconforming quadrilateral elements, Q9.The results obtained from GFEM uses a mesh M8, with polynomial enrichment of degree 2 p  .It must be observed that the results from Reddy (2004) are included here only to illustrate results available in the literature and to show a general trend.Both models have different number of nodes and types of approximation functions.Likely, the present model is more refined and accurate.The values of displacement and normal stresses are listed in Table 1 for the normalized values defined in the Figure 4 Figure 5 shows the evolution of the normalized in-plane stress x  along the line 0 y  , from the plate center 0 x  to the edge, at four different z levels, obtained with 0 C -GFEM with 2 p  , 250 P  and 25 nip  , and the mesh M8.Stresses obtained from the linear analysis are also shown, multiplied by 0.1.At this load level the nonlinear effect is pronounceable: the membrane stresses causes non zero stresses at the border / 2 x a  , differently from the usual zero value from the linear model; the laminate is symmetric and there is no coupling between bending and membrane behavior in the linear problem, and x  is zero for all z at the border in order to comply with the zero bending moment there.Along the range, from the center to near the border, the stress variation departs from the usual sine like shape of the linear case to an almost constant curve, also because of the nonlinearly generated membrane forces.Most of the transverse load is supported by the transverse component of the membrane stresses instead of the bending moment associated with the transverse gradient of x  .This can be seen by the difference between the values of x  at the plate center for the linear analysis, 76.4 x    , and nonlinear analysis, 1.69; 19.6) x   at the top and bottom surfaces respectively.

Cross-ply laminate problem
The problem considered here consists of a square laminated plate with sides 200 a b   mm aligned along x and y axis and thickness 2 h  mm, with three equal orthotropic layers with orientations [0 / 90 / 0 ]    .Each layer has the following orthotropic properties: Since the laminate is relatively thin, with aspect / 100 a h  , the behavior of the Kirchhoff model can be evaluated using conform Kirchhoff GFEM formulation with k C smooth functions.In this way, Figure 6 shows the transverse displacement w in the plate center versus the transverse load, using FSDT (Mindlin) with both 0 C and k C -GFEM and CPT (Kirchhoff) plate model using k C  GFEM.In all cases, the approximation basis has reproducibility polynomial degree 3 b  .C and k C -GFEM are very similar to each other, for the stresses at both surfaces.Figure 8 shows the normal stress x  across thickness at the center node, computed at two neighboring elements.Because the results are obtained with 0 C -GFEM, the discontinuity is expected, but it is very small in this case.The results for displacements and stresses for Mindlin and Kirchhoff are naturally similar because the laminate is thin.Also, the results from 0 C and k C -GFEM are similar for the same degree b of the approximation basis.For Mindlin model, a given b implies more degrees of freedom per node in the smooth than 0 C -GFEM, as seen in Table 2.
However, if the laminate is thin, k C -GFEM in Kirchhoff model requires equal or less dof's/node than 0 C -GFEM to attain a given degree b in the basis, for orders of 3 b  or higher.

Instability analysis
The problem considered here consists of a square laminated plate with equal sides 200 a  mm aligned along x and y axis, with thickness 2 H  mm, simply supported at all sides, with four equal orthotropic layers with orientations [0 / 90 / 90 / 0 ]     .Each layer has the following properties in their orthotropic directions: .Figure 9 and Figure 10 display results of relative errors of  versus number of degrees of freedom for Kirchhoff and Mindlin models respectively, with C k -GFEM.Each curve is associated to a given mesh, as indicated.The dots indicate the enrichment degree p of the basis, starting with 1 p  at the top, except for the first curve, mesh 1 M  , which starts from 2 p  .In this case, the discretization was too poor to allow a solution.We see the errors and the rates of convergence are similar for both kinematic models.All results were obtained with penalty  on G K equal to 2 10  times the smallest non zero value at the diagonal in G K .The shift factor was 10 s  .The eigenvalues were obtained with the subspace iteration method, with error tolerance   The abscissa scale in both Figure 9 and Figure 10 are the same.It can be seen a horizontal shift of the curves to the left in the Kirchhoff results, with respect to those of Mindlin model.This is natural because, for the first model, there is only one displacement component (w) and for the first order model there are three ( , x w  and y  ).
For the Kirchhoff model with enrichment 1 p  there is no convergence.This can be understood considering that only the transverse displacement is modeled.With 1 p  , the basis is too poor for allowing the adequate derivatives , x w and ,y w necessary to the model, although the space of approximation is still k C continuous, with arbitrary k .In the case of C 0 -GFEM, all displacement components, , x w  and y  , are modeled independently and convergence is observed.
For both, Mindlin and Kirchhoff models, Figure 9 and Figure 10 show very low convergence for enrichment degree 4 p  .This can be imputed to numerical difficulties due to the condition number of the stiffness matrix.
For the tests with C 0 -GFEM, as commented in the last section, the stiffness matrix contains an a-priori unknown number of spurious solutions with zero or very small eigenvalues.Figure 11   1.7.1 Comparison between GFEM, stable GFEM and C k -GFEM on the stability analysis The effects of the use of stabilization in C 0 -GFEM is evaluated using several meshes on the entire plate, in the same problem described at the beginning of this section, with  3 shows the estimated critical loads for the Mindlin model for polynomial enrichments 0,1, 2, 3 p  and 4, for both non-stabilized and stabilized C 0 -GFEM.The stabilization procedure (2.7), when applied to polynomial enrichment functions, automatically renders the stabilized linear function identically zero, such that remains in the approximation basis only the partition of unity function and the functions enriched by polynomials of degrees 2 p  and greater.Except for the very coarse mesh M1 (two elements and four nodes), C 0 -GFEM with 2 p  gives reasonable approximations for the critical load in all meshes, however, the first non-zero eigenvalue is the n -th one of the algebraic system.For , , 2, 3 p   the critical load is the 9th eigenvalue.In the stabilized case, C 0 -SGFEM, the critical load is the first eigenvalue for all p 's, and the more refined meshes, M4, M6 and M8.For the coarser meshes M1 and M2, the critical   The results in Table 3 are reorganized in Table 5 to show relative errors in critical load.Here the errors are given versus number of degrees of freedom for Mindlin model, using C 0 -GFEM, SGFEM and C k -GFEM, for several meshes.
To easy the interpretations, the results for two meshes, M2 and M8, were plotted in Figure 12.Every dot in a curve corresponds to a given degree p of polynomial enrichment, starting from 0, 1,2,3 or 4, as indicated in the figure.The relation between the results of C 0 and C k -GFEM follows as expected, because the difference between 1 b p   and b p  , respectively, generates a vertical translation between the curves, accompanied by a difference in convergence rates.However, it must be considered that the results for C k -GFEM were obtained computing only the first eigenvalue of the algebraic system, while for C 0 -GFEM the quantity s N of non-physical eigenvalues is unpredictable.Therefore, contrary to the tendencies of the convergence curves alone, C k -GFEM presents superior behavior in comparison to C 0 -GFEM for this type of application.
For C 0 -GFEM with mesh M8, it can be seen a lack of convergence between 3 p  and 4.This is imputed to the degradation of the condition number, which becomes more pronounced than in the mesh M2.The C k -GFEM results for M8 does not show an effect so intense, because, by construction, its approximation basis is linearly independent, which renders the coefficient matrix still ill conditioned but nonsingular.
The SGFEM aims at a stabilization of the condition number growth, by elimination of the linear dependency in the approximation basis by removing the linear content from the enrichment monomials.In most cases, as seen in n in Table 3, this is enough to renders a numerically nonsingular coefficient matrix, giving the critical load as the first eigenvalue.Even though, for coarse meshes M1 and M2, for 4 p  , it can be seen in the results a list of eigenvalues several orders smaller than the physical one.This happens as a result of the ill conditioning of the matrix, even after the stabilization procedure.
One side effect of the stabilization is the impoverishment of the approximation basis, as seen comparing the curves C 0 -GFEM and SGFEM in Figure 12.This has been observed with extensive numerical experimentation by Li (2014) in linear plane elasticity problems.It seems that the linear contents in the enrichment monomials, although generating linear dependency in the basis, provides also a useful enrichment, which is eliminated in the stabilization procedure.Judging by the curves in Figure 12, the resulting stabilized function basis is capable of spanning a space function similar to the one of smooth C k -GFEM: the curves are very similar in the range 2 p  and 3.

CONCLUSION
The current paper presents a comparison between two types of GFEM, 0 C and k C , in response of Mindlin and Kirchhoff kinematic models under moderately large displacements, modeled under the von Kármán theory.In addition, the response of both stable and non-stable 0 C -GFEM are evaluated in comparison with the smooth GFEM in the approximation of the buckling load, with regard to the effects of the stiffness matrix rank on the computational effort to obtain the first useful eigenvalue.The results allow the following observations and conclusions: 1) For Mindlin model, a given degree b of the polynomial reproducibility of the approximation basis implies more degrees of freedom per node in the smooth than 0 C -GFEM.However, if the laminate is thin, k C -GFEM allows the use of Kirchhoff model, which requires equal or less dof's/node than 0 C -GFEM to attain a given degree b in the basis, for orders 3 b  or higher, thus compensating the natural excess of dof's/node with respect to the 0 C version.
2) The convergence curves in pointwise relative errors of displacement for Mindlin model with 0 C -GFEM show stabilization of the curves for all enrichments, after some number of dof's, differently from the curves of smooth GFEM, which keep the convergence for all enrichment degrees, except for approximation basis of degree 1 b  .In fact, both forms of GFEM show locking in the nonlinear response for the Mindlin model with 1 b  .As for the stabilization of the curves for the 0 C -GFEM, it can be speculated it is related with bad conditioning of the stiffness matrix.The smooth version has the matrix naturally nonsingular, because the basis functions are linearly independent.method Paulo de Tarso Rocha de Mendonça et al.
Latin American Journal of Solids and Structures, 2019, 29(3), e228 22/24 3) As in the linear problem, the von Kármán response with the smooth GFEM requires a larger effort in the numerical integration than in 0 C -GFEM.This is due to the highly oscillatory aspect of the smooth enriched functions and their derivatives.
For the standard problem of the approximation of the upper bound of critical instability load by eigenvalue problem with k C -GFEM, we observe the following: 4) The convergence curves of the critical load by k C -GFEM and Kirchhoff model show lack of convergence for 1 p  , as expected because only the displacement component w is modeled with linear approximation functions and its first derivatives are too poor.For the Mindlin model, 1 p  generates a good convergence curve: the approximation basis is of degree 2 b  , applied over the three displacement components, w , x  and y  .
However, for thin laminates, the Kirchhoff model requires only one third of the dof's in the Mindlin model.
5) The behavior of the curves for enrichments 2 p  and 3 are excellent for both Kirchhoff and Mindlin models.For 4 p  , the response is not convergent, raising the hypothesis of ill conditioning in the coefficient matrix.
Considering the critical load obtained by 0 C -GFEM and SGFEM, we observe the following: 6) The critical loads obtained with 0 C -GFEM are, in general, as accurate as with smooth GFEM, but it corresponds to the n -th eigenvalue of the problem, where n is a priori unknown.This is different from the k C -GFEM, where all critical loads are the first eigenvalue of the model.
7) The stable form of GFEM, for non-coarse meshes, furnishes the critical load at the first eigenvalue, although with inferior accuracy.The processes of stabilization also impoverish the function basis.A positive aspect is that the number of dof's/node is reduced with respect to the non-stable GFEM because, in the process, the enrichment function of degree 1 p  naturally vanishes from the stabilized basis.However, with the impoverishment of the basis, the SGFEM requires more enrichment than GFEM to attain a given error.The increase in the number of dof's in SGFEM with respect to 0 C -GFEM, to obtain a given accuracy, is approximately the same as the increase in k C -GFEM, that is, both k C -GFEM and SGFEM show close convergence curves and produce the critical load at the first eigenvalue of the problem.

R
of elements)   associated to a node  as the union of all element domains that share the node  .In the present tests, the mid-surface of a laminated plate is modeled as an open bounded domain,   2 , Mindlin model involve five generalized displacements and Kirchhoff model three.For this type of equivalent layer kinematic model, the independent strains are the in-plane strains point of the plate.The hypothesis (3.1) allows the in-plane deformations to be represented in the separated in the form

E
are the elastic anisotropic material matrices of the layer, σ and s τ are in-plane and transverse shear stresses respectively.For the laminate, the constitutive relations are and e i u associated to node j , and jk L are the enrichment functions of this same node.The terms ij u are nodal values of e i u associated to the PoU functions and jk b are the coefficients of the enriched basis functions ij ijk L  .A compact notation can be obtained collecting the nodal parameters ij u and jk b in a vector e U and the corresponding basis functions in a matrix e N of dimensions dofe q N  representation (3.9) in the strain-deformations (3.2), one obtains the generalized linear deformations in discretized form where B is the deformation matrix associated with linear membrane deformation and change of curvature, and s B is the shear deformation matrix, used in case of Mindlin model.With regard to the nonlinear generalized deformation NL ε , here we follow the development introduced byZienkiewicz & Taylor (1991) for the von Kármán theory with Kirchhoff model, and perform the necessary adjustments to GFEM and the Mindlin model.First, one defines the matrices

Figure 1 -
Figure 1 -Examples of meshes used in the laminated plate numerical tests, with mesh indices M = 1, 2 and 4.

Figure 2 Figure 2 :
Figure 2 shows the pointwise relative errors of the transverse displacement at plate center, ( ) / r ex ex w w w w   , versus the number of dof, for 0 C -and 0 C -GFEM, for the different meshes M1 to M32, with uniform polynomial enrichment of degrees 0,..., 4 p  .The results are for the Mindlin model with shear parameter 1 s k  .The reference

rw
is represented versus Large deflection and initial instability analysis of anisotropic plates by the generalized finite element method Paulo de Tarso Rocha de Mendonça et al.Latin American Journal of Solids and Structures, 2019, 29(3are included, which are the smallest ones available in the integration rule of Wandzura.It is clear a tendency of convergence at 25 points for the present type of problems.

Figure 3 :
Figure 3: Relative error of the transverse displacement at plate center r w versus the number of integration nip.Mindlin model with Figure 4: (a) transverse displacement w and (b) normal stress x  versus transverse load P , 0 C -GFEM, enrichments 2 p  .and

Figure 5 :
Figure 5: Normal stress x  along the line 0 y  , at four different z levels at the thickness, for 0 C -GFEM with 2 p  and The laminate is thin, with aspect / 100 a h  .The boundary conditions are simply supported: at x and The tests are made with uniform transverse distributed loads ranging from 0.1 to 1 MPa.

Figure 6 :
Figure 6: Cross-ply laminate.Transverse displacement w in the plate center node, using FSDT and CPT plate models, 0 C -and k C -

Figure 7
Figure 7 shows the normal stress x  in the central node at top and bottom of the plate, versus transverse load, using FSDT and CPT plate models.The results are compared between the 0 C and k C -GFEM with 3 b  .The nonlinear results for FSDT model obtained with 0C and k C -GFEM are very similar to each other, for the stresses at both surfaces.Figure8shows the normal stress x  across thickness at the center node, computed at two neighboring

Figure 7 :
Figure 7: Normal stress x  in the central node at top and bottom of the plate, using FSDT and CPT plate models, 0 C -and k C - GFEM with b = 3. "b" and "t" refer to results at bottom and top surfaces of the laminate, respectively.

Figure 8 :
Figure 8: Normal stress x  through the thickness at the interfaces elements 1 first 2 eigenvalues.The size of the iteration subspace was set to 10 eigenvectors.

Figure 9 :
Figure 9: Relative error of  for Kirchhoff model versus number of degrees of freedom.

Figure 10 :
Figure 10: Relative error of  for Mindlin model with k C -GFEM versus number of degrees of freedom.
shows the number s N of zero eigenvalues obtained for each degree b of polynomial reproducibility of the basis, for Mindlin model.In all cases the 1 adequate approximation for the reference value of the critical load.The fraction of s N with the total number of degrees of freedom, / s N dof, grows approximately linearly with the degree b of the basis for all meshes, except for the coarsest one.

Figure 11 :
Figure 11: Number of null eigenvalues for Mindlin's model in 0 C -GFEM versus degree b of the polynomial basis. al.
proposed a procedure to generate smooth PoU by means of smooth weight functions.Let a function Paulo de Tarso Rocha deMendonça et al.
Latin American Journal of Solids and Structures, 2019, 29(3), e228 10/24 standard shift technique (solving the eigenproblem G . The data used here are the following: plate with sides Due to the symmetry, only one quarter of the plate is modeled.Considering the Cartesian axes 0xy with origin at the plate center, the boundary conditions for the Mindlin model are: at

Table 2 -
Number of degree of freedom per node in Mindlin and Kirchhoff models with 0 These geometric and material relations allow analytic solution of the critical load for the Kirchhoff model.Considering simply supported boundary conditions with

Table 4 -
Size of the basis nf for GFEM and SGFEM versus degree of enrichment p in two dimensions.

Table 5 -
Relative errors of critical loads obtained from 0 C -GFEM, SGFEM and k C -GFEM versus degree of enrichment p, for several meshes.Mindlin model.nsi is the number of subspace iterations, d1, d2 and d3 are the number of degree of freedom in the model.deflectionand initial instability analysis of anisotropic plates by the generalized finite element method Paulo de Tarso Rocha deMendonça et al.
Figure 12: Relative errors in critical load versus number of degree of freedom for Mindlin model, 0 C -GFEM and SGFEM and k C - GFEM.Every dot in a curve corresponds to a given degree p of polynomial enrichment, starting from 0, 1 or 2, as indicated.Large Latin American Journal of Solids and Structures, 2019, 29(3), e228 21/24