Energy based collocation method to predict progressive damage behavior of imperfect composite plates under compression

A new collocation methodology is presented to predict failure and progressive damage behavior of composite plates in this paper. The present work deals with composite plates containing initial geometric imperfections and different boundary conditions under uniaxial in-plane compressive load. In the present study, the domain is discretized with Legendre-Gauss-Lobatto nodes and the approximation of displacement fields is performed by Le-gendre Basis Functions (LBFs). The onset of damage and damage evolution are predicted by Hashin’s failure criteria and by proposed material degradation models. Three geometric degradation models are also assumed to estimate the degradation zone around the failure location which are named complete, region and node degradation models.


INTRODUCTION
Over the last decade thin-walled structures have extensively used in various lightweight structural components.If these structures be subjected to in-plane compressive loads, they are prone to buckling.Composite plates are one of these structures that can be optimized to carry in-plane compressive loads beyond their buckling resistance.They combine high strength with low weight, which makes them ideal for use in many industries such as aerospace engineering.Several laminate theories have been exploited to explain composite laminated plates' behavior.The Classical Laminated Plate Theory (CLPT) can be considered as an extension of classical plate theory based on the Kirchhoff hypothesis for isotropic plates and can be applied, if thickness of the laminate is small and the effects of transverse shear stresses are neglected (Dong et al. 1962;Yang et al.1966;Ambart ͡ sumi an 1970).An improvement on the CLPT leads to creation of the First order Shear Deformation Theory (FSDT) which accounts for the transverse shear effects (Whitney 1969;Pagano 1970) in relatively thick laminates.
For moderately thick composite plates with geometric imperfection, linear, nonlinear and post-buckling analyses without taking account of damage and also material degradations are very conservative analyses and give uncertain estimates of critical loads.Hence, damage analyses and ultimate strength of such structures have been of considerable research interest.In the field of buckling and post-buckling analyses, many researchers have extensively investigated these behaviors for composite beams and plates without considering the damage effects.Turvey and Marshall (1995) and Argyris and Tenek (1997) published excellent reviews on past studies of buckling and post-buckling of structures using different methods and recently, Pagani and Carrera (2017) have analyzed the large deflection and post-buckling behaviors of laminated beams by Carrera unified formulation.
Finite Element Method (FEM) is one the most common methods for analyzing the buckling and post-buckling of composite laminated plates.This method is based on dividing the plate's domain into a finite number of simple subdomains called elements.The behavior of each element is described by some specified functions.Finally, the solution of the whole plate as an assembly of its elements can be obtained according to the procedures applicable to the standard discrete problems such as those described by Zienkiewicz (1977) and Bathe (2006).Although FEM is an intensely powerful technique, it has some disadvantages.It needs large number of elements to reach acceptable stress results and in aspect of computational costs, it takes much longer computing time and requires larger amount of core storage.
Finite Strip Method (FSM) is another universally applicable method for buckling and post-buckling analyses of plates and plate structures.It can be considered as a particular kind of simplified finite element method in which a special element called strip is used.Finite strip method is based on discretization of the domain into longitudinal strips and interpolates the behavior in the longitudinal direction by different functions and in the transverse direction by polynomial functions.Cheung (1976) may be considered as the pioneer who first proposed the concept of FSM.Cheung established FSM for the analysis of simply supported plates.The studies presented by Smith and Sridharan (1978) proposed FSM for buckling of isotropic plate under edge loading.Recently, Ovesy, Ghannadpour and their co-workers (Ovesy et al. 2005) have made a contribution by introducing two different versions of finite strip methods, namely the full-energy and semi-energy finite strip approaches.Two other different versions of finite strip method, namely spline and semi-analytical methods are also developed by them for predicting the response of rectangular laminates with non-symmetric and symmetric forms of initial imperfection.They used both formulations to predict the non-linear response of channel sections when subjected to uniform end-shortening in their plane (Ovesy et al. 2006).An exact finite strip is introduced by Ghannadpour and Ovesy (2008) to investigate the exact relative post-buckling stiffness of I-section struts.To extend their works, they developed a high accuracy finite strip for the buckling and post-buckling analyses of moderately thick symmetric cross-ply composite plates based on FSDT (Ovesy et al. 2016).
Furthermore, many other methods employed to investigate composite plate's behavior subjected to compressive loads such as new class of numerical methods called meshless methods (Liu, 2009).These methods have attracted the attention of many researchers in recent years due to the fact that meshless methods do not require a mesh to discretize the problem domain as in the finite element method and require only a scattered set of nodes to model the domain of interest.Meshless or mesh-free methods have been proposed in multiple varieties such as Generalized Finite Difference method (GFD) (Liszka, 1984), which partial differential equations have been solved by applying irregular grids or clouds of points, the Smooth Particle Hydrodynamics (SPH) (Monaghan, 1988), a fully Lagrangian meshless method which at first was used for simulating astrophysical phenomena, the Diffuse Element Method (DEM) (Nayroles et al., 1992), which is the first meshless method employed moving leastsquares approximation, the Element-Free Galerkin method (EFG) (Belytschko et al., 1996) and Reproducing Kernel Particle Method (RKPM) (Liu et al., 1995).A corrected collocation method to apply essential boundary conditions in the meshless methods have been developed by Wagner and Liu (2000).
In the field of studies that investigate the plates behavior using mesh-free methods, it can be pointed out to the research conducted by Krysl and Belytschko (1995).They have implemented EFG to analyze the thin plates bending behavior.Lin and Jen (2005) proposed Chebyshev collocation method to solve the governing differential equations of a laminated anisotropic plate.A high order collocation method have been developed by Ferreira et al. (2009) for the static and vibration analyses of composite plates.Buckling of laminated composite plates subjected to various mechanical and thermal loads using meshless collocations have been studies by Singh et al. (2013).Free vibrations of beams were investigated by Carrera et al. (2013) using Radial Basis Functions (RBFs).They have analyzed the free vibration characteristics of beams on the basis of unified formulation.Liew and Huang (2003) and Liew et al. (2006) studied buckling and post-buckling of laminates using the moving least-squares differential quadrature and mesh-free kp-Ritz method.Also, Liew et al. (2011c) published a review paper for meshless methods and their applications in the buckling analysis of laminated or functionally graded plates.Recently, Ghannadpour and Barekati (2016) investigated the post-buckling behavior of composite plates using Chebyshev techniques by considering initial imperfection effects.More recently, Ghannadpour et al. (2017) have developed a high accuracy mesh-less analysis as a collocation method with Legendre Basis Functions (LBF) for nonlinear analysis of thin and moderately thick composite plates.
As emphasized earlier, all above mentioned research works, investigate the linear and nonlinear behaviors of plates without considering the effects of damage or failure and therefore the obtained results may be unrealistic.On the other hand, as was observed, all mentioned studies have been developed by different numerical or semianalytical methods.However, strength analyses of metal plates and damage analyses of composite plates are usually carried out using FEM in which the storage required is extremely large, and the computational time is too lengthy.Only a limited number of such analyses have been implemented by other numerical or semi-analytical methods.
Several simplified semi-analytical methods have been purposed by Brubak and Hellesland (2007a,b, 2008, 2011) and Brubak et al. (2007c) for investigating the strength of stiffened and unstiffened metal plates with imperfection and cutout.Accelerated analysis techniques were proposed by Orifici et al. (2008) using FEM tools such as ABAQUS to model the degradation and fracture mechanics of composite stiffened fuselage panel in a postbuckling regime.Some experimental works with comparison by FEM analysis can be found on ultimate strength of composite plates such as study investigated by Hayman et al. (2011).These studies are included a parametric study of ultimate strength analysis of rectangular simply supported composite plates with different geometrical imperfections.Yang and Hayman (2015a,b) and Yang et al. (2013) recently carried out a valuable studies established simplified methods to predict ultimate strength of composite plates with different initial imperfection values under in-plane compressive loads.
In the present paper, a simplified progressive damage methodology is developed to predict the ultimate strength of imperfect composite plates with different boundary conditions under compression.The simplified methodology is referred to the assumptions of small deflection theory by which the computational time is greatly reduced.A new methodology is presented based on collocation method in which the interested domain is discretized with Legendre-Gauss-Lobatto nodes and therefore this approach do not require a mesh to discretize the problem domain.However, with less computational efforts, an acceptable field of stress can be obtained which is necessary for damage analyses.The formulations are based on the concept of the principle of minimum potential energy and the approximation of displacement fields is performed by Legendre Basis Functions (LBFs).In this new collocation approach which is based on the total potential energy, there is no need to enforce the natural boundary conditions, and the essential boundary conditions are satisfied very easily.Therefore, it is easy to apply different boundary conditions.The structural model is based on the first order shear deformation theory.The onset of damage and damage evolution are predicted by Hashin's failure criteria (Hashin and Rotem, 1973) and by proposed material degradation models.Three types of degradation models are assumed to estimate the degradation zone around the failure location which are named complete, region and node degradation models.With these assumed models, the accurate results can be extracted by reducing the area of failure around a node, and this is done very easily and quickly due to the proposed formulations.The properties of damaged materials instantaneously reduced to a 1% of primary value.Some examples involving various boundary conditions, initial imperfections and different thickness to width ratios are investigated to demonstrate the validity and capability of the proposed method.The accuracy of the present work is examined by comparing the numerical results with the previous studies.

Theoretical formulations
A typical rectangular laminated plate of dimensions of a b  and total thickness h with an initial geometric imperfection i w at the center of plate and in the z-direction is shown in Figure 1.The laminate is subjected to inplane compressive load on the edge / 2 x a  in the x-direction called x N .Two different types of boundary conditions are considered that they will be mentioned later.The laminates are assumed to be moderately thick, thus the formulations are based on the first order shear deformation theory (FSDT).To described the deformation of a laminated plate, displacement vector and expressed by equation (1).
Whose components of displacement vector are given by equation (2). Where . According to the FSDT (Reddy, 2004),   , , u v w are in-plane and out-of-plane displacement of mid-plane, and x  and y  denote the rotations of a transverse normal about axes parallel to the x and y axes, respectively and w is initial geometric imperfection of the laminate in the form of a single half sine-wave in both directions.Due to moderately large displacements assumption, the vector of Green's strainse in a total Lagrangian formulation is written as equation (3).
Where p e and n e are in-plane and out-of-plane component vectors with subscripts p and n , respectively.By substituting geometric nonlinearity in the von Karman sense, the strain-displacement relations are given as: where ⊗ denotes the Kronecker product.Also equation ( 4) defines the plate strain vectors that are the in-plane strains vector 0  , the curvatures vector  and the shear strains vector  .Since the structural model in this study is based on the assumption of small deflection theory, therefore the nonlinear and imperfection strains vectors   are neglected henceforth.Differential operators p  and n  can be defined by equation ( 5).
Latin American Journal of Solids and Structures, 2018, 15(4), e35  5/25 The mechanical state is described by the stress vector  , which can be written as equation ( 6).
On the assumption that the plate is in a state of plane stress 0 zz   , the constitutive equations for the ori- ented k th orthotropic ply of the laminate are written as equation ( 7). ; The plate internal actions are including, namely, the in-plane stress resultants N , the moment resultants M and the transverse force resultants T .All these stress resultants are measured per unit length and can be obtained by equation ( 8).
In which p A , B , D and n A are the generalized stiffness matrices, which are extensional stiffness matrix, ex- tensional-bending stiffness matrix, bending stiffness and interlaminar shear stiffness matrices, respectively.They can be computed by the following equations.
Where K is the shear correction factor.By applying the external force vector on the plate, the total potential energy Π is consists of the strain energy of the plate  and potential energy of external forces  .Therefore, the total potential energy can be obtained by Latin American Journal of Solids and Structures, 2018, 15(4), e35 6/25 Where  is the domain occupied by the plate mid-plane.It is emphasized that in this paper the plate is only subjected to in-plane compressive load x N .
   However, as it can be seen, the labelling schemes are included in these figures to assign the related boundary conditions.The letter S refers to simply supported boundary condition and the letter C refers to clamped boundary condition on the specified edge.Therefore, the two aforementioned boundary conditions are different from the viewpoint of out-of-plane boundary conditions.The details of edges conditions for both Type A and Type B are expressed in Table 1.It is noted that the words "Free" and "Held" denote whether the edge is free to move or restricted against any movement, respectively.Also, the word "straight" denotes the edge can move longitudinally (for u ) or laterally (for v ) but remain straight.
Latin American Journal of Solids and Structures, 2018, 15(4), e35 7/25 Since the Rayleigh-Ritz technique is used in this study, approximated displacement fields have to satisfy only the aforementioned essential boundary conditions.The approximation of displacement fields is performed by Legendre basis functions P .Legendre polynomials are one of the most powerful mathematical series for numerical methods.Legendre basis functions or Legendre polynomials are solutions to the following Legendre differential equation is written as equation ( 12).
Also Legendre polynomials satisfy the three-term recursion expressed by equation ( 13). Where Therefore, the displacement fields of the problem can be approximated by equation ( 14).
Where   , , , , The so-called boundary function is also chosen to ensure the fulfillment of the essential boundary conditions mentioned in Table 1.It can be defined as equation ( 16).
Where  denotes the edge number and the exponents    can take the value 0 for free condition and the value 1 according to the conditions of held or straight for each displacement field   , , , , Specifically, as noted earlier, to the aim of the present work Legendre orthogonal polynomials have been chosen.It is worth mentioning that different types of polynomials can be used, whose characteristic features have been discussed in the literature.
As can be seen in right-hand side of equation ( 14), each displacement field

 
, , , , 4 Equilibrium equations and solution procedure As mentioned before, the equilibrium equations are obtained based on the concept of the principle of minimum potential energy.Therefore, the elements in equation ( 11) should be written as described in previous section.For this purpose, the vectorsu ,  and u can be rewritten in compact matricial form as equation ( 18).
where , , , , U U and U are defined as equation ( 19).
Accordingly, the strain vectors defined in equation ( 4) with the assumption of small deflection theory can also be rewritten as equations ( 20) and ( 21). Where And finally by using equations ( 20), the total potential energy  can be shown as: Latin American Journal of Solids and Structures, 2018, 15(4), e35 9/25 Since in the current study, the domain of the plates is discretized by a set of nodes therefore the above continuous integrals should be replaced by summations where they can be calculated over all nodes.For this purpose, Legendre-Gauss-Lobatto nodes are established here and can be obtained by solving equation ( 23).In order to achieve better accuracy and also to avoid excessive number of nodes to reduce computational costs, an appropriate weight coefficient can be considered for each node.Calculation of the weight coefficients for nodes is performed by taking idea from Gauss-Lobatto rules and therefore, the continuous integral of total potential energy (i.e.equation ( 22)) after eliminating the constant factors is then converted to the following relation:

   
Latin American Journal of Solids and Structures, 2018, 15(4), e35 10/25 Where ,   indicates the th  node along x direction and th  node along y direction as represented in Figure 3.
The coefficients   and   are weight coefficients of nodes along x and y directions, respectively and they can be computed by equation ( 25).
To obtain the equilibrium equations of the problem using the principle of minimum potential energy, the discretized form of the total potential energy equation ( 24) should be minimized with respect to the unknown primary variables u U and  U (i.e. ).Therefore, the final set of equilibrium equations can be written as equation ( 26) and ( 27).
, , 5 Progressive damage model methodology In this section, the methodology of progressive damage analysis including the failure criteria, material degradation model and ply geometric degradation models is described in details.

Failure criteria
In order to determine the failure load and the corresponding failure mode, a proper failure criterion must be established.The selected failure criteria applied in the present study are proposed by Hashin and Rotem (1973).
The Hashin failure criteria used herein include four different damage functions which correspond to the different modes of failure namely fiber tension, fiber compression, matrix tension, and matrix compression.These criteria can be presented as equation ( 28).
Fiber failure in tension 1 0   : Latin American Journal of Solids and Structures, 2018, 15(4), e35 11/25 Matrix failure in tension 2 0   : Here T X and C X denote tensile and compressive strengths of fiber and T Y and C Y are tensile and com- pressive strengths of matrix, respectively.Failure happens when any of these modes reaches unity.

Material degradation model
As it is known, when damage occurs in a composite structure, the effective material properties are reduced.This reduction can be modeled in this study by the following matrix.
  Where 1 E , 2 E , 12 G , 12  and 21  in equation ( 29) are undamaged material properties.The parameters f d and m d are the damage factors in fiber and matrix directions, respectively.Note that for the Hashin criterion, because the shear failure component is associated with both the fibre and matrix modes of failure, the shear damage factor is assumed to be not independent and given by the remaining damage factors.
According to degradation model presented by Hayman et al. (2011) and also in order to better comparison with the results obtained by Yang et al. (2013), the transverse shear stiffness matrix n Q is not degraded during the damage analysis.
In the current progressive damage model, the reduction of material properties is considered to happen instantaneously.When failure is detected in specific zone (complete, region or node), its properties are instantaneously reduced to %1 of their initial and undamaged values (i.e.f d or 0.01 m d  ).

Ply geometric degradation models
In this study, three geometric degradation models are assumed to estimate the degradation zone around the failure location which are named complete, region and node degradation models.When failure occurs in a location, the material properties of its zone should be changed.6 Numerical results and discussion In this study, in order to implement the proposed formulations for analyzing the progressive damage of composite plates, a computer program is developed based on Fortran 77 software pack-age.To reduce the execution time on multicore processors, a parallel programming technique is used using Open Multi-Processing (OpenMP) interface.The program is executed on a computer with TYAN FT48-B8812 mainboard, 4 AMD CPU by 2.20 GHz frequency (4×16 cores) and 128.00 GB RAM.
The mentioned proposed formulations should be verified through a number of comparisons and this is done by comparing the results with those obtained by Hayman et al. (2011) and Yang et al. (2013).For this purpose, the plates which are considered in this study are ( ) with two types of boundary conditions as described before.Also, a uniform compressive load, x N , is applied in the x-direction at / 2 x a  as mentioned earlier.Initial imperfection is considered to be a single half sin wave in both directions by maximum values of 0.1%, 1%, 2% and 3% of the length a in center of the plates.The composite plates under consideration are modeled as , 0 45 90 45 distinct layers whose mechanical properties are listed in Table 2 as taken from Yang et al. (2013).The parameter X takes the values of 2, 3 and 4 in this study.The thickness of each ply is considered to be equal to 1 mm and therefore, the composite laminated plates have total thickness of 16, 24 and 32 mm.Since in this study, one of the geometric degradation models is region degradation model (RDM) in which a plate is divided into 9 regions, so specifying the size of each region is necessary.Therefore, in order to better comparison between the results, the size of the regions has been chosen as those considered by Yang et al. (2013): regions dimensions 1, 3, 7 and 9 are 160×160mm.Regions 2 and 8 are 180×160mm each.Regions 4 and 6 are 160×180mm and region 5 is 180×180mm.
To do the convergence analysis, two composite plates with total thickness of 16 and 32 mm are selected.The boundary conditions are imposed as described for Type A and the plates have the maximum initial imperfection of 3% and 2%. Figure 5 shows the convergence study for both plates with the number of terms t N in the dis- placement fields.
Latin American Journal of Solids and Structures, 2018, 15(4), e35 13/25 Figure 5: convergence study of 16 mm and 32 mm composite plates with 3% and 2% initial imperfection respectively As can be seen, for the plates under consideration the convergence studies with regard to the number of terms have revealed that 8 terms are sufficient to obtain converged results.However, of 9 terms is used to ensure accurate convergence in all analyses.Thus, the total number of unknown coefficients is 407.Also, similar analyzes were conducted with regards to the number of nodes for all three geometric degradation models and it was concluded that the total number of 169 nodes (m×n=13×13) sufficient to obtain converged results.

Results for CDM.
To investigate the results for complete degradation model (CDM), two types of boundary conditions are considered as shown in Figure 2. The results for composite plates with geometrical and mechanical properties mentioned in previous section and for CDM are given in Tables 3 and 4 for different maximum initial imperfections and total thicknesses.Other given data in these tables are first ply failure (FPF) stress, number and angle of that layer and the location failure, last ply failure (LPF) stress and also number and angle of last failed ply.Total number of failed plies is also given in these tables.Ultimate strength of composite plates with boundary condition Type A is tabulated in Table 3 and they are compared by results reported by Yang et al. (2013) and the results computed for boundary condition Type B are presented in Table 4. Also, the load-center out-of-plane displacement history   x tc N w  and load-end shortening response   x c N u  are graphically represented for composite plates having total thickness of 16 and 24 mm and for two types of boundary conditions A and B in Figures 6-9.The parameter tc w is defined as the value of total out-of-plane displacement in the center of the plate (i.e.

   
0, 0 0, 0 w w  ).As can be seen in these figures, the results for different values of initial imperfection are also incor-porated in order to better comparison and observation.
Latin American Journal of Solids and Structures, 2018, 15(4), e35 14/25  in the lower-most ply.The only exceptions for both types are the two thicker plates with small initial imperfection (0.1%).For these special cases, the first ply failure occurred in ply with o 90 fiber orientation angle.Another important point is the location of first failure for both composite plates with different boundary conditions.For composite plates with all simply-supported edges (Type A), first failure usually occurs in the middle of the plates while for plates having one clamped edge (Type B) damage correspond to the first ply failure load mostly initiates from the mid-point of the clamped edge.
It is also seen that the ultimate strength of the composite plate is usually attained at the first incidence of fiber failure especially for plates having one clamped edge.Furthermore, the investigations show that for thin plates with small initial imperfections, there are little reserve strength beyond the first ply failure load while the ultimate strength of the thicker plates or the plates with larger initial imperfections have higher values than first ply failure loads.In these laminates, all plies usually have to experience matrix failure.
It is also seen from both tables and figures that by increasing the maximum value of initial imperfection of the plates, both FPF and LPF loads are decreased.Furthermore, as it can be observed from tables and figures, ultimate strength of composite plates having at least one clamped edge, is usually much greater than those plates with boundary condition Type A.

Results for RDM
In this section, the results for region degradation model (RDM) are presented and investigated.Tables 5 and  6 show the results for composite plates having boundary conditions Type A and Type B, respectively.Similar to the section of CDM results, first ply failure stresses and locations have been achieved.Also in both tables, last ply failure stresses and the number of matrix and fiber failed regions are included for each laminate thickness and initial imperfection.Furthermore, the values of ultimate strength obtained by CDM (Tables 3 and 4) are also retabulated here for better comparison.It can be observed that there is no difference between the first ply failure loads obtained by both CDM and RDM.However, the CDM model gives generally slightly lower ultimate loads than the RDM model.This difference is slightly higher for plates with boundary condition Type B.
In Figures 10 and 11, the applied load is plotted against the central out-of-plane displacement and end shortening displacement respectively, for plates having total thickness of 32 mm and 2% initial imperfection.In order to better comparison, the results for both types of boundary conditions A and B are also represented in these figures.It is seen that the results obtained for plates with boundary condition Type B have 10-20% upper ultimate loads with respect to the results obtained for Type A. This is due to the fact that the boundary conditions have significant effects on the response of imperfect composite plates and as mentioned earlier the last ply failure loads of plates having one clamped edge are greater than the corresponding values in plates with all simply-supported edges.
However, regarding the first ply failure, for plates with small initial imperfection, Type A gives lower values than the results of Type B but by increasing amount of initial imperfection, FPF loads take higher values than Type B. In general, the laminates with boundary condition Type B can carry higher loads after their first ply fails.
Moreover, considering the location of first ply failure, plates with boundary conditions Type A fail in the center of the plate as mentioned in previous section, while laminates with boundary condition Type B fail in near the center of clamped edge.

Results of NDM.
The results for node degradation model (NDM) are presented and discussed here.Only the results correspond to the plates with boundary conditions Type A are presented in this section.These results are tabulated in Table 7.As for previous models, first and last ply failure loads are reported in this table.In addition, the number of failed nodes in each fiber or matrix mode are also given.In order to better discussions, the results obtained by previous models (Tables 3 and 5) and those reported by Hayman et al. (2011)are tabulated in the separate columns.
Latin American Journal of Solids and Structures, 2018, 15(4), e35 20/25 The behaviors of load against central out-of-plane displacement and load against end-shortening displacement for plates with    As it can be observed, there is no difference between the values of first ply failure loads calculated by three degradation models however, the complete model gives generally slightly lower values of ultimate load than the two others.
With comparison between three degradation models in Table 7 it is seen that the node degradation model gives upper last ply failure load and this can be due to the effects of dividing the plate (or ply) into smaller areas.Therefore, if one wants to obtain better and more accurate results, the node degradation model should be used in which more computation time is also needed but in the early stages of structural design, the results taken from region model can be accepted with slightly lower accuracy.
However, as it is seen, the ultimate loads predicted by all geometric degradation models in this study and those reported by Yang et al. (2013) are still smaller than those calculated by FEM analysis.This is mainly due to neglecting of the nonlinear terms in strain-displacement relations and therefore post-buckling effects, which are particularly considerable for thin plates.But for thick plates even by neglecting of those effects, the results have less deviation from FEM.

Conclusions
A new methodology has been introduced in this study to investigate the progressive damage analysis of composite plates with initial imperfection under in-plane compressive load.The concept of the first order shear deformation theory and the assumption of small deflection have been established to drive the equilibrium equations.In presented method, the domain of the plate is discretized with Legendre-Gauss-Lobatto nodes.The onset of damage has been predicted by Hashin's failure criteria and material properties of damaged zone have been degraded by instantaneous material degradation model.Three geometric degradation models were assumed to estimate the degradation zone around the failure location.The results were compared with linear analysis reported by Yang et al. (2013) and also non-linear FE analysis presented by Hayman et al (2011).By comparison, it has been concluded that the obtained results by increasing the thickness of the plates and decreasing the area around the location of failure can take better accuracy.

Figure 1 :
Figure 1: A typical imperfect rectangular laminated plate 11) 3 Displacement fields and boundary conditions In order to approximate the displacement functions of the problem, it is necessary to describe the specified boundary conditions.The laminates under consideration have two different types of boundary conditions as represented in Figure 2. As shown in this figure, for both types of boundary conditions, the in-plane displacement inxcompletely restricted and in the opposite edges uniform movement are allowed and then all four edges of the plates are being held straight.

Figure 2 :
Figure 2: Two different types of boundary conditions


containing the Legendre and boundary functions and a column vector   contain- ing the corresponding Ritz unknown coefficients.With this definition and also in order to comply with other notations, the previously mentioned initial geometric imperfectionw that has a shape of single half sine-wave in both directions in this study, can be generally represented by equation (17).
m and n denote the number of nodes in both x and y directions, respectively and x  and y  are non-dimensional coordinates of th  and th  node in the x and y directions.Figure3represents a scattered set of Legendre-Gauss-Lobatto nodes (

Figure 6 :Figure 7 :Figure 8 :Figure 9 :
Figure 6: Response of Load vs maximum out-of-plane displacement for 16 mm composite plates (Type A)

Figure 10 : 25 Figure 11 :
Figure 10: Response of Load vs maximum out-of-plane displacement for 32 mm composite plates with 2% initial imperfection (Type A and Type B) 16 h  and 32 mm and for two different values of initial imperfection are depicted in Figures 12-15.In these figures, in addition to the results obtained by NDM, the results taken from two other models, CDM and RDM, are also represented.

Figure 12 :
Figure 12: Response of load vs maximum out-of-plane displacement, compression between degradation model for 16 mm composite plate with %1 initial imperfection (Type A)

Figure 13 :
Figure 13: Response of load vs end shortening, compression between degradation model for 16 mm composite plate with %1 initial imperfection (Type A)

Figure 15 :
Figure 15:Response of load vs end shortening, compression between degradation model for 32 mm composite plate with %3 initial imperfection (Type A)

Table 1 :
Details of edges conditions for both Type A and Type B and y are non-dimensional coordinates defined as 2x a and 2y a , respectively and t N is the number of terms in series expansion which is taken same for all displacement fields.The coefficients ij  and c   are the Ritz unknown coefficients and the latter is for satisfying the straight conditions mentioned in Table1for in-plane displacements    is defined as equation (15).

Table 2 :
Material properties of fiber and matrix (in MPa)

Table 3 :
Results for progressive damage of composite plates using CDM and for boundary condition Type A

Table 4 .
Results for progressive damage of composite plates using CDM and for boundary condition Type B

Table 5 .
Results for progressive damage of composite plates using RDM and for boundary condition Type A Latin American Journal of Solidsand Structures, 2018, 15(4), e3518/25

Table 6 .
Results for progressive damage of composite plates using RDM and for boundary condition Type B

Table 7 .
Results for progressive damage of composite plates using NDM and for boundary condition Type A max 