Progressive Stiffness Loss Analysis of Symmetric Laminated Plates due to Transverse Cracks Using the MLGFM

With the purpose to create high strength advanced structures, new materials are being developed presenting favorable characteristics for specific applications. Composite Materials are examples of these developments. They can be formed by high strength long fibers, conveniently oriented in a matrix, to form a lamina of composite material. The lamina presents high strength in the fiber direction, but, since it is slender, does not have enough rigidity, what makes impossible the use of an isolated lamina. Piling up and gluing a set of laminas, a laminate is formed which one presents better characteristics than original and isolated materials. The main strength of each lamina is oriented according to the fiber directions. Thus, micro cracks can be produced if sufficient tension is applied in the transverse direction of the fibers, as shown in figure 1, since the resistance of the lamina in these directions depends only on the matrix material. The rise of several transverse cracks produces loss of stiffness in the laminate.


Introduction
With the purpose to create high strength advanced structures, new materials are being developed presenting favorable characteristics for specific applications.Composite Materials are examples of these developments.They can be formed by high strength long fibers, conveniently oriented in a matrix, to form a lamina of composite material.The lamina presents high strength in the fiber direction, but, since it is slender, does not have enough rigidity, what makes impossible the use of an isolated lamina.Piling up and gluing a set of laminas, a laminate is formed which one presents better characteristics than original and isolated materials.The main strength of each lamina is oriented according to the fiber directions.Thus, micro cracks can be produced if sufficient tension is applied in the transverse direction of the fibers, as shown in figure 1, since the resistance of the lamina in these directions depends only on the matrix material.The rise of several transverse cracks produces loss of stiffness in the laminate.
Several papers are found in technical literature dealing with the behavior of composite materials with transverse cracks.Vejen & Pyrz (2002) investigated the transverse crack growth in long fiber composites using the finite element method.Three criteria concerning pure matrix growth, fiber/matrix interface growth and crack kinking out of a fiber/matrix interface were implemented to form a software package for crack propagation calculus.Cain and colleagues (2003) have studied unidirectional graphite bismaleimide composites to determine the effect of the matrix dominant properties on the failure of the material.The authors showed that the final fracture was caused by the development of a dominant matrix shear crack parallel to the fibers.They also concluded that the decrease in shear modulus of the composite was the most sensitive and best represented by damage evolutions.Ogihara et al (1998) have proposed a two-dimensional model which considers that, in the case of displacements and stress fields in the interlaminar cross-play laminates, there is a prevalence of plane-strain case, even in the presence of transverse cracks.They have also commented that the failure process of cross-ply laminates is due an accumulation of transverse cracks and delamination.An analytical model based on the principle of minimum potential energy was developed by Ji et al (1998) and applied to determine the two-dimensional thermoelastic stress state in cross-ply composite laminates containing multiple equally spaced transverse cracks in the 90 o plies subjected to tensile loading in the longitudinal direction.The criterion of strain energy release rate was employed to evaluate the critical applied stresses for two of the possible fracture modes.After some numerical experiments, they have concluded that the formation of new cracks never takes place until pre-existing cracks extends through the entire thickness of the 90 o plies.Wada et al (1999) have presented a damage mechanics model to predict the nonlinear behavior of laminated composites due to crack evolution.A new concept of cracking layer is proposed by a technique based on uniform work-softening layer.With this concept, the constitutive equations for a cracking layer are constructed according to modern plasticity theory.So, the lamina damage surface is defined in the stress space and the constitutive equations for a cracking layer are constructed by applying the defined damage surface to the associated flow rule.One of the first damages that occurs in composite laminate are the transverse cracks, as mentioned by Allen and colleagues (Allen et al 1987 a,b) and Tay and Lim (Tay & Lim, 1996).The cracks appear in a layer where the highest stress values act transversally to the fiber, exceeding the matrix resistance.With loading increment, the increase on the number of transverse cracks may happen in a diffuse way reducing the structural rigidity.The accumulation of this damage can accelerate the beginning of delamination, changing the natural frequency of the composite structure and causing a greater degradation in severe environment, jeopardizing its service life.
After the initiation and development of micro cracks, there is a process of accumulation of damage that reduces the structural stiffness.The tolerance to the damage is related to the stiffness of the structure which, in turn, is affected by the accumulation of micro defects during loading.The process of damage evolution in composite laminate is generally very complex due to the multiplicity of failure modes such as transverse cracks, delamination, decoupling fiber-matrix interface, and fiber breakage.The characterization of this process is generally possible when single cases are analyzed, where each failure mode can be separated and studied individually.The use of Fracture Mechanics, especially in terms of linear elastic fracture, has presented good results for isotropic material because, in this case, can be adequately characterized by a single parameter (the stress intensity factor).However, attempts to apply this method in composite laminates, whose behavior is orthotropic, have met unsatisfactory results, mainly when transverse cracks in the matrix are studied.Therefore, to determine changes of the mechanical properties in a laminate, the total number of cracks formed in the transverse layers must be taken to account, or, under a generalized crack distribution, the most appropriated methodology is based on Damage Mechanics.
Many researchers have developed studies to evaluate the properties of laminates subjected to generalized cracks in the matrix.Among these ones, can be cited the papers of Allen et. al. (1987 a,b), Hashin (1987), Talreja (1984) and Lim & Tay.
The present paper has the objective to apply the Continuum Damage Mechanics Theory to long fiber laminate composites.The transverse cracks appearance in the matrix implies in a rigidity loss due to damage accumulation.The increase of the load is considerate monotonically.Several failure criterions are presented and implemented such as, the Maximum Stress Criterion, the Maximum Strain Criterion, Tsai-Hill and Tsai-Wu Criterion.
The proposed methodology is restrict to the case of symmetric laminate and it is evaluated by a numerical approximation technique known as Modified Local Green's Function Method (MLGFM), which one will be briefly descript on this paper.

Constitutive relations
The models developed by Talreja & Boehler (1990), Allen et. al. (1987 a,b) and Lim and Tay (1996) to describe the damaged composite laminates were based on the Continuum Damage Mechanics using internal state variables.In the presente paper, the model proposed by Allen et al (1987 a,b) will be used, which describes the damage through a set of internal state variables.The final result of the distributed damage is built in the constitutive equations through these variables.Thus, the stress-strain relationship of the representative volume of a damaged material at the level of a lamina is assumed as: ijkl ijkl However, it is important to emphasize that equation ( 1) does not provide any information on how the damage state has been attained, that is, the history of damage accumulation.Thus, it is necessary to turn to Fracture Mechanics in search of a suitable criterion to evaluate the damage growth.Thereby, equation ( 1) is sufficiently general to permit the use of Classical Laminate Theory to determine the composite laminate constitutive relations with transverse cracks in the matrix.
Supposing the representation of the laminated plate by plane elements located in its middle surface, the loads in a certain point inside this surface can be evaluated by the following expressions: where {N}e {M} are, respectively, the force and moment resultants vectors,  x ,  y ,  xy are the stresses in the plane of the lamina and t is the thickness of the laminate.Taking to account that { 0 } e {κ 0 } are the strain and bending vectors in the middle surface of the plate, [A], [B] and [D] are the laminate extensional stiffness matrix, coupling stiffness matrix and bending stiffness matrix, respectively, {D N }e {D M } are the damage vectors related to the force and moment resultants, the expressions (3) and ( 4) can be transformed to: Assuming that z k-1 e z k are the corresponding distances from the middle surface to the inner and outer surfaces of the k th lamina, respectively, [] k Q and {} k  are the transformed reduced material stiffness matrix and the transformed vector of the internal state variables (expressed in global coordinates), respectively.The matrix and vectors presents in equations ( 5) and ( 6) can be expressed by: The internal state variables vector has two components and is expressed by: where P is the total number of damage models being considered and  22 and  12 are the internal state variable of the problem.

Determination of internal variables
In spite of the random characteristic, as can be found in the work of Silberschmidt ( 2005), the transverse cracks are assumed to be uniformly distributed.In this way, the laminate behavior can be adequately represented by a representative unit volume of material containing a transverse crack, as shown in figure 2. In the particular case of symmetric laminates, the damage models are simplified and incorporate only two types of fracture, namely, Mode I (crack opening) and Mode I coupled with Mode III (shear out of plane).They are represented respectively by the internal variables  22 e  12 .As only symmetric laminates are analyzed in this paper, just the α 22 variable will be developed.
The internal variable, equation ( 15 where u 2 is the crack opening displacement, n 2 is the unitary vector normal to the crack surface, V is the representative element volume and S c is the crack surface.Lim & Tay, 1996) Considering that t 1 and t 2 are the thickness of the 0° and 90° plies, respectively, t is the total thickness of the laminate, l is the distance between two adjacent transverse cracks, ρ is the non-dimensionalized crack density (the quantity of cracks per unit of length), is the nondimensionalized maximum crack opening displacement, u 2 is the maximum crack opening displacement, ψ is a normalized function of the crack opening profile and  is the normalized distance between the cracks center, the following expressions can be defined: The maximum crack opening displacement u 2 can be determined by a simple finite element analysis, considering the boundary conditions specified in figure 3.An arbitrary displacement is imposed.The value of is determined by the equation ( 18), and the displacement u 2 is determined by MEF.The non-dimensionalized maximum crack opening displacement can be obtained using ρ, as shown in the expression (20): The constants a 1 , b 1 , c 1 , c 2 e c 3 in the expression (20) depends on the type of the used material.
The table 1 exposes the value of these constants for the laminate glass/epoxi (Gl/Ep).

Material
Formulation in terms of ρ Glass/Epóxi 1.03 -0.81 2.28E-2 0.94 1.00 Table 1.Coefficients for the expression (20) (Lim & Tay, 1996) As the internal variable used in this problem depends on the maximum crack opening displacement according to equation ( 15

Approximation by computational methods
In this paper, the maximum crack opening displacement u 2 is determined by the expressions ( 18) and ( 20).The crack density ρ depends on the distance between two adjacent transverse cracks l, and its values are successively modified by the verification of the composite material rigidity loss.
Considering a conventional structural approximation by conventional Finite Element Method, the problem can be expressed by a system of algebraic equations, representing a typical element: where K ij are the element stiffness matrix, (d x , d y , d z ) are the components of the element displacement vector and {F a } and {F d } are the applied force vector and the element damage force vector.
The procedure used in this paper to obtain the expected results is a little different because it uses a different computational method known as Modified Local Green's Function Method (MLGFM), in witch the system defined in expression ( 22) is not directly applied in a conventional FEM.A detailed explanation of this procedure can be found in Barbieri et al (1998a,b) and Machado et al (2008).The MLGFM is an integral method that determines the unknowns on the boundary, similarly to the Boundary Element Method, but the fundamental solutions are generated automatically by projections of the Green's Functions developed from de field, as in the Finite Element Method.The matrix and the vectors indicated in ( 22) will be used to produce values at the boundary, as explained in the next topic.

The Modified Local Green's Function Method -MLGFM
The Modified Local Green's Function Method (MLGFM) is an integral technique that associates the Finite Element Method and the Boundary Element Method, solving the problem through an integral equations system.Unlike to the BEM and the Trefitz Methods, the MLGFM does not use a fundamental solution and/or a Green's function.The term "Local" indicates that the calculation of the GF projections can be done locally, that is, for each element.
Essentially, the MLGFM uses a transverse integration technique and reciprocity relations to determine, at a local level, the Green's Function, transforming the partial differential operator in an ordinary partial operator (Barcelos & Silva, 1987).The MLGFM uses finite elements at the domain with the purpose to create discrete projections of the Green's Function, corresponding to fundamental solutions, that are used later in the integral equations system associated to the boundary approximation.To analyze a continuum mechanics problem through the MLGFM, such as the plate indicated in figure 4, two meshes are necessary, one for the domain and other for the boundary.With domain elements, the method generates a set of domain equations, which are used to generate automatically the domain Green's projections.Later, a set of boundary equations are also generated and the boundary Green's projections can be determined with the domain projections developed before.At the end, the system is solved only for boundary equations, where the main variables are calculated.Domain values may be obtained once the boundary values are known after the solution of the boundary equation system.The most important steps of the MLGFM are detailed in the work sof Barbieri et al (1998a) and Machado et al (2008).It is possible to show that through the MLGFM two sets of equations are formed, the first one in the domain (equation ( 23)) and the other one on the boundary (equation ( 24)):  [G T (P,q) a(P)]d +   [G T (p,q) f(p)]d ; P  ; p, q  (24) where Q, P are two points in the domain; q, p are other two points on the boundary; a(P) is the vector of independent terms for the original problem; f(p) is the vector associated to the fluxes on the boundary; G(i,j) are the Green's functions which may be understood as the generalized displacement in the point i in the direction of an unitary vector n i , when a generalized force is applied over the point j, in the direction of a unitary vector n j .
Equations ( 23) and ( 24) describe completely the problem.Since these equations involve domain and boundary integrals, two types of meshes are necessary, one in the domain and the other on the boundary, using FE and BE methods, respectively.The FE domain approximation is also used to develop the Green's functions which are associated to the matrices G(P,Q), G(p,Q), G(P,q) and G(p,q).
The approximation shape functions are the same as in the conventional FE and BE methods.
For the present work, quadratic shape functions were employed to construct nine nodes lagrangean finite elements and three nodes boundary elements.
Developing discrete equations from the nodal values, two sets of linear equations are determined: A u  = B f + C a (in the domain) ( 25) where u  and u  are the domain and the boundary displacements, respectively, a and f are the independent and the fluxes variables vectors.The matrices A, B, C, D, E and F can be written as: where (Q) and (q) are matrices with the shape functions in the domain and on the boundary, respectively, and G  (Q), G  (Q), G  (q), G  (q) are the Green's function projections over the boundary  and the domain , evaluated on the points Q and q.The Green's projections can be written as: In order to determine the Green's functions automatically, it must be considered the following functional F, which one depends on G  or G  : As in the variational approach of the FEM, the minimization of functional F(G  ,G  ) in Equation ( 37) results in a linear equation system which can be solved to determine the Green's projections where [K] is the global stiffness matrix, evaluated in the same way of the conventional finite element stiffness matrix; A and D are the matrices of Equations ( 27) and (30), respectively.In this way, the Green's projections are determined directly from Equation ( 41), and can be applied in Equations ( 28), ( 29), ( 31) and (32) to complete the matrices of equations ( 25) and ( 26), which are the main system of the MLGFM.

Damage evolution
With purpose to quantify the damage accumulation due to a monotonic load increment, some failure criterion will be used.Generically, failure criteria can be considered as: where F a is the failure criterion, ij  are the local stress, X ij are the principal material strength and Z is the failure value characteristic to each criterion.The rigidity degradation of a component occurs due to a progressive process during its serviceable life.It is important to note that the evolution of rigidity loss in a structure can be characterized by a single crack or by the occurrence of generalized cracks.Combinations of failure modes can act together causing changes in the material properties and in its local stress distribution.In this way, the main difficulty in this kind of analysis is the adoption of a failure criterion, a F , that conveniently describes the damage evolution due to a failure mode.
The theories introduced to prevent the failure of an orthotropic laminate are adaptations of failure criterion for isotropic materials, modified for biaxial stress cases, such as, Maximum Stress Criterion, Maximum Strain Criterion, Tsai-Hill Criterion and Tsai-Wu Criterion (Reddy, 1997;Vasiliev & Morozov, 2001;Mendonça, 2005).
A new criterion is presented, based on the strain energy release rate, to evaluate the formation of a new micro crack (Anderssen et al., 1998;Ji et al., 1998;Kobayashi et al.,2000).The released energy is used because it is practically independent from the crack length (Anderssen et al., 1998).Some of these criterions are presented here.

Maximum stress criterion
According to the Maximum Stress Criterion, for orthotropic materials, while the stresses in the principal directions of the material are lower than strength of the material in this direction, there are no fails, which means: where t X is the longitudinal tensile strength, t Y is the transverse tensile strength, c X is the longitudinal compressive strength, c Y is the transverse compressive strength and C is the shear strength of the lamina.

Maximum strain criterion
This theory is analogous to the Maximum Stress Criterion, but the fail criterion is controlled by deformation limits in the principal directions of material.In this theory the material will fail when one of the following limits are reached: where X, X, Y, Y are the maximum deformation deformation values in the principal directions 1 and 2, for tensile and compressive loading, C  , is the maximum angular distortion in the plane 1-2.

Tsai-Hill Criterion
An adaptation made by Tsai in the Hill Criterion for transverse orthotropic laminate at plane stress condition, resulted in the expression (43):

Strain Energy Release Rate Criterion
The Strain Energy Release Rate Criterion to describe the damage evolution was applied by Lim and Tay (Lim & Tay, 1996).Consider a symmetric laminated composite of width b and length L, with the configuration [0°l/90°m] s , where l and m are integers.When the laminated is loaded uniaxially in tension, the stress-strain curve is linear until the failure criterion is reached for the first time, at point A (figure 5).A transverse crack is introduced in 90° layer.are the strain in x direction, y direction and xy plane of the lamina i, respectively.Therefore, the energy i U , necessary to form a new crack, can be defined as: Where t is the thickness of the laminated, t 2 is the thickness of the 90° plies an L is the length of the laminated.
In this way, a transverse crack is assumed to be formed when: Where Ic G is the mode I energy release rate for the formation of a transverse crack.The process of determining the transverse crack density is repeated for each successive micro crack, using the same value for Ic G .As seen in figure 5, a series of points (A, B, C, D, …, E) can be generated until the limit OI is reached.From this limit, matrix cracking in the 90º layers no longer influences significantly the laminate stress-strain behavior.It must be observed that in practice, the intervals between the points are very small, turning the curve smooth, rather than the curve shown in figure 5.

Analysis of laminated plates by the MFLGM
The first application refers to the analysis of a laminate plate, whose materials of its lamina are defined in table 2. The aim of this application is to determine the stiffness loss E/Eo due to the improvement of crack density  for the Gr/Ep [0 o ./90o ]s laminated, using () formulation  The loss of stiffness is observed in figure 6 for different meshes and compared with the results obtained by Lim & Tay (1996) and experimental results.As the crack density grows up, the stiffness diminishes.The results are better with finest meshes, but even with coarse meshes, the approximation is good.Figure 7 shows the loss E/Eo versus crack density ζ for the case Gl/Ep [0°/90°] s -using () formulation.The same considerations are made for this case.

Progressive stiffness loss of laminate
To evaluate the stiffness loss of laminated plates due to micro-crack accumulation under increasing monotonic loading using the MLGFM, the following conditions were considered: a. Stress-strain relations of a thin orthotropic laminate are considered in plain stress state; b.Dimensions of the squared plate are 2,0 m x 2,0 m, but only a ¼ was modeled due to its double symmetry: 2 {( , ) :(0 1,0;0 1,0} xy R x y     ; c. Axial tension loading in "x" direction; d.The properties of the material used are listed in the tables 2 and 3; e.The value of Ic G adopted is 250 J/m 2 for the glass/epoxi laminate (Tay & Lim, 1993).(Highsmith & Reifsnider, 1982) In order to compare the failure criterion, a [0 o /90 o 3 ] s glass/epoxy symmetric laminated with total thickness of 1,624mm was used.All layers on the laminated have the same thickness.
The results are presented in figure 8.All criterions were implemented in the same program to facilitate the comparison.To demonstrate the model capability to prevent the laminate stiffness loss, a comparison between the results obtained by the strain energy criterion implemented in this paper and the results obtained by Talreja (Talreja,1984) and Tay (Tay & Lim, 1993) was made.This analysis also used the [0 o /90 o 3 ] s glass/epoxy symmetric laminated with total thickness of 1,624mm.The results are presented in figure 9.

Conclusion
The present paper deals with damage composite laminate with transverse cracks in the matrix applying Continuous Damage Mechanics Theory, which was initially proposed by Kachanov (Kachanov, 1958) and than adapted by Allen (Allen et al., 1987a,b) for orthotropic laminated composites.This theory was also applied by Lim and Tay (Lim & Tay,1996) in laminates with transverse cracks to describe the stiffness loss of the structure.The adapted Damage Theory considers the mechanism associated to the transverse cracks through the internal state variables inside the constitutive relations based on the Continuous Damage Mechanics.
The theoretical model was implemented in a computational program, developed in FORTRAN language, based on the Modified Local Green's Function Method (MLGFM).The approximated solution was obtained by the MLGFM.The damage evolution model, originally developed for FEM, can be applied also to MLGFM without substantial changes in the original code.
In the presented results, it can be observed that the conventional criterions catch only the moment when the 90º layers no longer influences the stiffness of the laminated.Most of the criterions were able to determine the loss of stiffness.The strain energy criterion is able to evaluate the damage evolution, identifying the moment when the transverse cracks starts to affect the laminated rigidity.However, during the strain increase, the efficacy of the method to evaluate the stiffness loss decreases.Even so, as shown in the figure 8, the implemented

Fig. 1 .
Fig. 1.(a) A composite [0 o /90 o /0 o ] s laminate plate with transverse micro-cracks in the matrix; (b) extension of pre-micro cracks; (c) formation of new micro crack.
is the applied stress tensor, C ijkl is the constitutive relation tensor of the undamaged material, kl  are the strain tensor, ijkl I  are the elements of the damage matrix, kl   are the internal state variables, and  = 1, 2, 3, …, refers to the damage modes.As suggested by Allen et al (1987 a,b), a first simplification can be made considering that the tensor I ijkl is the actual tensor of constitutive relationships, as shown in Equation (2).
), proposed by Allen(Allen et al., 1987 a,b) can be determined by a computational analysis based on Finite Element Method.The representative volume is modeled, as shown in figure3, for the symmetric laminate [0 o /90 o /0 o ].A uniform displacement is imposed in one side of the element to determine the opening of the crack.The size and shape of the representative volume depend on the thickness of the different plies and the crack density (the number of cracks per unit of volume).Then, the internal variable can be determined by:

Fig. 4 .
Fig. 4. Symmetric boundary conditions and plate for a 2x2 mesh (4 finite elements and 8 contour elements) by the MLGFM.
) where G -corresponds to G  or G  depending on the case of interest; B -is a bilinear form, developed to G  or G ;  and  are constants whose values are:  = 1 and  = 0 to determine G   = 0 and  = 1 to determine G  B 1 , B 2 , B 3 -are bilinear forms which can be written as: was proposed by Tsai-Wu, changing the Tsai-Hill Criterion in equation (43).When the tensile and compression strength are similar, the expression (44) becomes the Tsai

Fig. 5 .
Fig. 5. Stiffness loss in composite laminates [0°l/90°m] s -Lim & Tay (1995).The result is a reduction in the effective stiffness of the laminate in the loading direction, and this is represented by the OB segment in figure5.Upon further loading, this reduction is verified by the segment BC.This process is repeated until line OF is reached.Note that

Table 4 .
Strength limits for glass/epoxy laminate