Statistical Analysis of the Progressive Failure Behavior for Fiber-reinforced Polymer Composites under Tensile Loading

An analytical approach with the help of numerical simulations based on the equivalent constraint model (ECM) was proposed to investigate the progressive failure behavior of symmetric fiber-reinforced composite laminates damaged by transverse ply cracking. A fracture criterion was developed to describe the initiation and propagation of the transverse ply cracking. This work was also concerned with a statistical distributions of the critical fracture toughness values with due consideration given to the scale size effect. The Monte Carlo simulation technique coupled with statistical analysis was applied to study the progressive cracking behaviors of composite structures, by considering the effects of lamina properties and lay-up configurations. The results deduced from the numerical procedure were in good agreement with the experimental results obtained for laminated composites formed by unidirectional fiber reinforced laminae with different orientations.


Introduction
Laminated polymer composites are of considerable interest in a spectrum of structural applications on account of their superior mechanical performance, such as light weight, high specific strength and OPEN ACCESS stiffness.To ensure the structural reliability of composite materials, it is an urgent task to comprehensively understand their failure mechanisms under thermo-mechanical loads and to develop a methodology that enables their failure behaviors to be predicted.
Performance in service of a composite structure is influenced by the progressive occurrence and interaction of some or all of many multi-mechanisms of damage, such as matrix cracking, delamination, fiber fracture and fiber/matrix debonding under loading and environment conditions.The first damage that occurs in laminates is transverse matrix cracking, the formation of intralaminar cracks running parallel to the fibers in the off-axis piles at a much lower stress [1].These matrix cracks develop in the fiber direction, and three-dimensional intersectional cracks across multi-layers and along multi-directions will be formed when the density of the matrix cracks increases with the increased applied load or fatigue cycles.Subsequently, interlaminar delaminations will be developed and extended at the tip of matrix cracks [2].Matrix cracking is usually referred to as transverse cracking, because the crack plane is commonly transverse to the laminate middle plane.This can result in considerable impairments of the performance of laminates and triggers the development of other more harmful damage mechanisms [3].Therefore, the stress in the primary load-bearing lamina will be redistributed, which will lead to the ultimate failure of composite structures [4].
The detailed mechanisms of transverse ply cracking in laminates have been a subject of active research efforts.A large variety of analyses have been conducted to investigate the damage phenomena of composite laminates: different modifications of the shear-lag model [5], a variational approach based on the minimization of complementary energy [6], a continuum damage mechanics approach [7], an efficient experimental technique [8] and numerical methods, for example, based on finite element analysis [9].The overwhelming majority of studies investigating the behavior and properties of composite laminates with matrix cracks assume that cracks are equally spaced [10].Such a deterministic approach, which ignores the fact that transverse cracking is a progressive damage mode, predicts the appearance of many transverse cracks simultaneously when the first transverse cracking strain is reached [11].In fact, the transverse matrix cracking is an inherently stochastic process due to the random variations of the local material properties of the plies caused by inner original defects [12].The random distribution of these defects is the major factor that affects both crack initiation and propagation.Prediction on the durability of these materials tends to be probabilistic in character.Modeling of transverse cracking in laminated composites should, therefore, consider the details of its construction and the use of statistical analysis.
Failure processes in real composites are of a great variety, which is why many papers investigating the strengths of particular types of composite laminates are of quite special interest.In order to present detailed reviews, a worldwide failure exercise was launched by Hinton et al. [13,14], where a great number of researchers had described their own theory and its intended application in order to compare the accuracy of these theories for predicting the degradation of laminate properties and failure, however, in which many involved models were valid only for a particular laminate lay-up, and most predictions were based on a ply-by-ply analysis with certain failure criteria [15][16][17].
In the previous papers [18,19], the reduction of the stiffness properties and the energy release rates due to transverse ply cracking for the balanced symmetric composite laminates were predicted by the equivalent constraint model (ECM).Within the framework of the ECM (Figure 1), the nearest neighboring plies of the damaged lamina under consideration (layer k) are called the two primary constraining layers.All the laminae below and above these piles, which are called the secondary constraining layers, are replaced with homogeneous layers (I and II) with the "layer properties" having the equivalent constraining effect.The nicest feature that the ECM contributes is introducing damage variables that depend on the strain state.The relation between the stress and strain field in the cracked lamina can be determined through the energy equivalent principle of the interstitial medium; thus, the in situ damage effective functions for the expression relating stiffness degradation to the transverse crack density can be derived for any general symmetric lay-up laminates.

Figure 1.
A schematic diagram of (a) a multilayer crack pattern in a multidirectional laminate; and (b) An equivalent constraint model (ECM) laminate for the matrix cracking in the k-th lamina (for details, see [18]).

(a) (b)
In the present study, we summarize our research work and propose an analytical approach that has the capacity of predicting the progressive failure of material structures damaged by transverse matrix cracking subjected to quasi-static loading, appropriate for the mechanical analysis of fiber reinforced composite laminates [20].The energy release rate associated with matrix cracking in the transverse layer of general symmetric composite laminates is presented using a 2D shear lag stress analysis and the ECM [2,21].An energy criterion, based on the fracture mechanics approach, is proposed to predict and describe the ignition and propagation in symmetrical composite laminates subjected to in-plane loading.This work is also concerned with a stochastic description of the inner defects with due consideration given to the scale size effect.The random properties are introduced with the use of a Weibull two-parameter probability density function.A numerical procedure is developed to investigate the multiplication of transverse cracking and to predict the statistical failure strengths of composite structures by taking into account the statistics of critical fracture toughness.The results are also compared with those obtained from experiments and other existing models.

Analytical Model
Micro-cracking has been a mechanism of failure extensively studied in many investigations.This mechanism occurs under transverse tensile or longitudinal shear loading conditions.A known analytical modeling of mechanical properties due to transverse cracking was performed by using the suggested ECM to determine the in-plane stress distributions in the damaged layers [22,23].The reduction in stiffness properties due to transverse ply cracking was described by two in situ damage effective functions introduced in [18]:  , called in situ damage effective functions (IDEFs), can be expressed as: where the matrix cracking damager parameter, D, can be written in the form of cracking spacing.ø i , Γ i and λ i are laminate constants and are functions of the lamina elastic properties and the geometrical parameters, such as the laminate stacking sequence.A detailed description of these parameters was provided in [21].Thus, these two damage parameters,    , are determined from the obtained stress field under the framework of shear-lag arguments.They characterize the in-plane stiffness reductions of the constrained, cracked lamina.The findings of some work show that a crack lamina behaves within a laminate in a different manner compared to an infinite effective medium containing many cracks [24].
Our simulation reproduces the initial crack and its propagation in the thickness direction using the energy criterion.Indeed, the potential energy method has been presented by Zhang et al. [21,23] for the evaluation of the energy release rate for damage initiation and propagation due to transverse ply cracking.The required energy release rate, G, for damage growth can be calculated using the energy approach and then predicting that damage will grow when G exceeds the critical energy release rate, G c , or the toughness of the material.
The energy release rate, G, associated with the appearance of a crack for a given stress state according to fracture mechanics is defined by the following expression: where PE stands for the potential energy of a damaged laminate element under the condition of the plane strain.A is the crack surface area.σ is the applied laminate stress.For the sake of convenience, the damage state of n cracks (c 1 , c 2 , …, c n ) is denoted by the parameter, D n , which does not necessarily mean that the damage state can be characterized by a single parameter.
For illustration purposes, a Cartesian coordinate system is introduced, in which the x-and z-axes are, respectively, parallel and perpendicular to the fibers in the lamina that fails first.The origin is located at the center of the lamina, as shown in Figure 2. The laminate specimen has a gauge length of 2L.2h is the thickness of the cracked lamina.Assuming that matrix cracking occurs randomly due to the statistics of fracture toughness, we define the system that has the ply cracks numbered one after another, where c j is the time coordinate of the j-th crack when appearing.b j is the x-coordinate of the j-th crack.The j-th crack appears right after the (j-1)-th crack has formed.By applying an opening displacement of the transverse crack (cracking opening displacement, COD) solution for a crack, one can obtain the total energy release rate due to matrix cracking, as follows: where   for -L ≤ c n ≤ b 1 and: for b n−1 ≤ c n ≤ L and: The parameters, χ and λ, are specified as laminate constants (see [25]).With the solution of the energy release rate, a damage criterion is proposed to predict the evolution of matrix cracking.The estimated quantity is compared with the associated critical value of the stain energy release rate: where G c denotes the fracture toughness of the material.

Modeling of Progressive Cracking
The origins of transverse cracking are the inherent material defects, such as the microcracks, voids, debonded fibers, areas of high fiber volume fraction, etc. [12].This causes the transverse layer to have a statistical nature of the fracture toughness along its length.Thus, the transverse cracking propagation can be investigated by assigning a random distribution of the fracture toughness along the transverse layers.
In the present paper, the crack multiplication can be simulated by dividing the initial gage length into equal elements along the direction of the length, 2L, and randomly assigning fracture toughness G c to each of them in accordance with a two-parameter Weibull distribution.
Based on the Weibull statistics, the fracture toughness is given as: where G 0 represents the scale parameter or the characteristic quantity of the material.The Weibull modulus, or shape parameter m, controls the degree of disorder in the distribution, experimentally found to describe a variety of materials.F in the range of zero to one is the probability that the fracture toughness of the material equals G c ; thus, the fracture toughness, G ci , of each segment can be obtained by generating uniform random number F in the range of [0,1].Figure 3 presents the effect of m on the shape of the probability density function when G 0 = 350 J/m 2 .This finding indicates that m variations significantly affect the distribution of fracture toughness, i.e., the smaller the shape parameter, the more scattered the facture toughness and vice versa.In the approach we use here, transverse cracking at some location takes place when the energy released by a cracking event becomes equal to the critical strain energy release rate, G c , at this position.In fiber-reinforced composites, the largest portion of the loads is resisted by the fibers.When matrix cracking occurs, the internal loads must redistribute to other areas of the structures and may cause a structural collapse [4].For every state of stress, a simple criterion, such as the Hoffman criterion, is expressed mathematically in the following fashion: in which σ 1 , σ 2 and τ 12 are longitudinal stress, transverse stress and shear stress, respectively.X t and X c are the tensile strength and compressive strength of the unidirectional layer parallel to the fiber direction.Y t and Y c are the tensile strength and compressive strength of the unidirectional layer transverse to the fiber direction.S is the shear strength of the unidirectional layer transverse and parallel to the fiber direction, respectively.It should be mentioned that in the present study, the ultimate failure of the materials takes place when the stress state of the primary load-bearing lamina satisfies the condition mentioned [26].

Numerical Algorithm
The simulation is conducted by controlling the stress loading of the model.Figure 4 shows the progressive failure analysis algorithm for estimating the damage growth of laminated composites under combined loading: Step 1: The ply material and ply orientation are selected to form a laminate.The laminate properties can be determined indirectly from mechanics analysis based on the classical laminate theory.The critical strain energy release rate, G c , is assigned to each element.Step 2: A virtual crack is introduced in all positions as possible sites for failure, and the work performed to close the crack surfaces [27] is calculated.New transverse cracking happens in any location once the energy released by a cracking event is equal to the critical value.The virtual cracking procedure is repeated unceasingly, until cracking is terminated under the same applied stress.
Step 3: The equivalent constraint model is employed to analyze the damaged lamina.In the ECM, all the laminae below and above the damaged lamina under consideration are replaced with Step 2 Based on ECM Degraded mechanical properties of damaged lamina Evaluation of the stress distribution in the primary load-bearing lamina Step 3

No
Failure of composites Step 4 homogeneous layers [18].The reduced stiffness properties of the damaged layer can be calculated by applying the laminated plate theory to the ECM.Thus, the in-plane microstress of the primary loading-bearing lamina can be obtained by the constitutive relationship.
Step 4: With the stresses calculated, they are substituted in failure criteria Equation (7) to check for failure.When the loads are increased monotonically, the matching strains are computed, and the resulting stresses are substituted in the failure criteria until they are satisfied.The ultimate load of the laminate is thus determined.

Results and Discussion
A detailed analysis is conducted to assess the predictive capabilities of the present methodology.The mechanical properties of each type of unidirectional laminate used in this work are given in Table 1.It should be noted that the Weibull parameters may be obtained by fitting experimental cracking data for laminates [28,29].In order to give a comparison between predicted values and experimental ones, the simulation of the progressive cracking behavior is carried out based upon the assumed Weibull parameters derived from the references.
Six different laminate configurations (denoted by A-F) are studied as shown in Table 2.Each set of Monte Carlo simulation consists of 300 data points.The element length of 2.5 × 10 −2 mm is chosen for the sake of efficiency [11].

Prediction of Stiffness Reduction
Transverse matrix cracking is a progressive damage mode that evolves with the increase in the applied strain, and a measure of accumulated damage is the crack density.Thus, the reduced elastic properties of a laminate, such as a longitudinal Young's modulus and the major Poisson ratio, are related to the number of cracks formed in the transverse layer.The comparison between the present predictions and experimental data [30] is shown in Figure 5 as plots of normalized axial modulus against crack density (number of cracks per centimeter) for the [0 2 /90 4 ] s glass/epoxy laminate.
Figures 6 and 7 show the reduction of the Poisson's ratio for the AS4/3501-6 and CFRP [0/90] s lay-up laminate as a function of matrix crack density.It should be mentioned that the numerical result of each case is the average value of 300 data.The results reveal that the theoretical predictions provided by the analytical approach agree well with the experimental values [31,32].Indeed, the probabilistic model for characterizing the transverse matrix cracking can emphasize the statistical nature of damage accumulation [33].

Analysis of Progressive Cracking
For the purpose of verifying the accuracy of the present model and demonstrating its predictability for the progression of matrix cracking, Figures 8 and 9 show the results obtained for the crack density as a function of the stress applied to the glass-epoxy laminates with a stacking sequence of [±θ/90 4 ] s , where θ = 15° and 30°.G 0 = 750 J/m 2 and m = 10 are used as the scale parameters in the Weibull distribution.The comparison between the deterministic model and the probabilistic model shows that the deterministic model predicts a very rapid accumulation of ply cracks succeeding the first cracking, whereas the probabilistic model presents a smooth increase of the crack density with the applied stress.The latter is closer to the observed values in the experiments [11].From these figures, we can see that the transverse matrix cracking procedure can be well described by the fracture toughness distribution.The statistical model can effectively characterize the initiation and propagation of matrix cracking.Moreover, the model is able to predict the constraining effect, namely, the cracking stress required for a given crack density gets smaller with an increase of the lay-up angle, θ, of [±θ/90 4 ] laminates [25].

Estimation of Ultimate Strength
Figures 10 and 11 show the uniaxial strengths as a function of the angle-ply, θ, expressed in degrees predicted for the T300/5208 [0/±θ/90] s and [0/90/±θ] s laminates, respectively.The test results from [34] are used to validate the present model.Due to the effect of the constraining action of the adjacent layer decreases as the degree θ increases, the ultimate strengths for the laminates become weaker with an increase of θ.Furthermore, a comparison of the experimental results with the current analysis predictions is made, which exhibits a good agreement.It is interesting to note that locating the 90° ply adjacent to the laminate mid-plane will lower the laminate strength.The reason for this is almost certain to arise from the fact that in such a case, the length of the ply cracks in the through-thickness direction is double that for the other laminates, leading to a larger stress concentration in the 0° plies and to fiber failures at lower stresses [35].The observed offset of the ultimate strength is due to the lack of the micro-delamination mechanism in this approach.Especially for thicker layers, delamination at the interface will occur earlier during the loading.One feature worthy of note, which is demonstrated by the deterministic model predictions [26], is that the current theory correctly takes account of the effects of the variability in fracture toughness due to microflaws distributed randomly along the length.

Conclusions
The objective of this exercise is to approach the modeling of progressive damage and failure in composite laminates.A statistical model based on the computation of the strain energy release rate associated with matrix cracking is proposed to study the stiffness reduction and to predict the ultimate strength of in symmetric laminated composite.The statistical distribution of the critical energy release rate, G c , in the transverse layer is described by a two-parameter Weibull function in this article.Several sets of application examples and comparison with experimental results show that the present numerical model is able to reproduce the mechanical behavior of the laminated composite formed by unidirectional fiber reinforced laminae with different orientations.

Q 22 
are the stiffness of an intact lamina and Q (k) are the degraded stiffness matrices of the k-th damaged lamina.  k and   k 66

Figure 2 .
Figure 2. Schematic of the lamina containing matrix cracking.

0 2 σ
x is the x-axis normal stress in the xy-plane and the superscript, 0, is used to indicate the quantities that belongs to the undamaged state. 0ij Q are the in-plane stiffness components of an undamaged lamina.g(L, D) is the normalized energy release rate, which can be expressed as:

Figure 3 .
Figure 3.The Weibull distribution of fracture toughness.

10 Figure 4 .
Figure 4. Flow chart for the analysis of the progression of transverse cracking.

1
Input data: -Dimensions of composite laminates -Mechanical properties of each lamina Start Weibull statistics Distributions of the fracture toughness of the transverse plies along their length Step Calculation of the energy release rate of the damaged layer Introduce virtual crack G

Figure 5 .
Figure 5. Normalized axial modulus as a function of crack density in the 90-degree layer.

Figure 6 .Figure 7 .
Figure 6.Normalized Poisson's ratio as a function of crack density in the 90-degree layer.

Figure 8 .Figure 9 .
Figure 8. Crack density in the 90-degree layer as a function of the stress applied to the composite laminates.

Figure 10 .
Figure 10.Comparison between the theoretical and experimental results angles for degrees.

Figure 11 .
Figure 11.Comparison between the theoretical and experimental results angles for θ.

Table 1 .
Basic properties of unidirectional composite material systems.

Table 2 .
Laminate configurations of composite structures.