A Continuum Damage Mechanics Approach for Reliability Analysis of Composite Laminates Containing a Central Circular Hole

This paper discusses the topic of probabilistic analysis of composite laminates rectangular plates containing a central circular hole under static tensile load. First, According to the continuum damage mechanics (CDM) approach, structural analysis models, the definition of damage and the treatment of random variables will be explored. Then, the matrix cracking and fiber/matrix debonding damage is modeled and the material constitutive relationships are implemented in the ABAQUS software by the subroutine. The probabilistic methods, first order reliability method (FORM) and second order reliability method (SORM) have been used to analyze the composite plates’ system failure probability. The failure functions and random variables have been obtained according to the CDM approach. Issues arising out of the use of composite material structures, in applications, due to non-homogeneity and anisotropic characteristics, and manufacturing defects will be highlighted, and uncertainties related to the material properties, loads, and boundary conditions will be discussed through examples. It is shown that the scatter in parameter values has a significant effect on damage development in the model.


Introduction
Composite laminates materials reduce the structure's weight without reducing strength and are used in many engineering structures. The use of holes in these structures is inevitable and reduces the strength of composite laminates material. According to the loading, environmental, and boundary conditions, the local damage initiates in the early stages of loading and propagates over time.
Investigating the behavior of composite materials involves many uncertainties, such as variation of the fiber and matrix volume fraction, cavities in the matrix, incomplete bonding between components, cracks, fiber damage, load distribution, boundary conditions, geometry and structure of laminate, and variation in the stacking sequence, etc. These ambiguities lead to significant random changes in the mechanical performance of composite structures. As a result, probabilistic analysis is often used as an alternative to the conventional deterministic methods.
Gosling et al. [1] used the first order reliability method (FORM) to assess the shells' reliability based on the plates and shells equation. They confirmed their results by Monte Carlo Simulation (MCS). Using the Tsai-Hill failure criterion, Zhang et al. [2] specified the statistical correlation effect between the properties of plies on the reliability of composite structures. The effects of ply thickness uncertainty on panels made of symmetrical composite materials have been discussed by Zhang et al. [3] in uniaxial and multi-axis loads and calculated the probability of failure by using the MCS and the Markov chain according to the Tsai-Hill criterion. Zhou et al. [4] incorporated the FORM and SORM to assess composite laminates materials' reliability. They provided a method according to the stochastic finite element and obtained a random structural response. Nadjafi and Gholami [5] studied the reliability of rectangular plates containing a central circular hole under static tensile load using the CDM approach for ductile fracture. In this study, to investigate the initiation and evolution of damages, anisotropic damage expressed by second order damage tensor is used to derive constitutive equations. Yaghoubi and Gholami [6] developed an explicit formula for the computation of the reliability of a system with dependent component.
In recent decades, damage models relying on the finite element method (FEM) have found the most attention [7][8][9][10]. These models can directly reveal the damage and fracture of composite materials for any type of stacking and material. The CDM approach in the formulation of the FEM is increasingly developed [11][12][13]. The CDM models are based on thermodynamic potential (Gibbs free energy function) to express material constitutive equations. In recent decades, the CDM model suggested for composite laminates by Ladeve`zeet al. [14,15] has been of interest to researchers in the area of damage in composites. This model has been developed for various conditions, including compressive loading of composite laminates [16], three-dimensional damage analysis on a reinforced plate with a stringer [17], and has been used to calculate the behavior of damage in composite materials subjected to low-velocity impact [18,19]. Gholami et al. [20] proposed a piecewise fatigue damage model (PFDM) for the prediction of damage in composite laminates under cyclic loading based on the CDM model. In the other work, Gholami et al. [21] developed a stochastic fatigue damage model to predict the fatigue life of fiber-reinforced laminated composites by CDM-based damage plastic model.
The main objective of this research is to estimate the failure probability of the rectangular composite laminates containing a circular hole under uniaxial static tensile load using reliability methods taking into account the mechanical properties and model uncertainties. Namely, the FORM and SORM are used to assess reliability anaysiss. According to FORM and SORM, a real failure function is required to analyze reliability. For this purpose, the CDM approach is used. By determining the random variables and the failure function relying on the CDM approach, the probability of failure made of S2-Glass/Epoxy composite laminates will be investigated. Using this model, the laws of stiffness degradation of the plies can be calculated base on irreversible thermodynamics.

Reliability analysis
Traditional structural assessment involves safety analysis in the substantiation of a special design standard. It is a simple method to confirm that the failure probability of the structure is acceptably small. In these models, the problem variables are given a certain value, and safety factors are applied. However, in structural reliability methods, the probability of failure is estimated more accurately using statistical distributions.Reliability analysis is performed similar to the traditional method, with the difference that each parameter of the problem is determined by statistical distributions.
In this method, first, the various processes of structural failure, basic variables, and failure criteria must be determined. After determining the boundary conditions, the failure functions are given relying on the determinant problem and the statistical distribution of the fundamental problem variables. Accordingly, a possible method for estimating structural reliability has been selected.
In this study, the FORM and SORM are used to calculate the failure probability of composite laminates containing a central circular hole under in-plane tensile load using the CDM approach.

First Order Reliability Method (FORM)
The FORM is mostly used in structural reliability assessment. This approach includes the Taylor expansion of the failure function, in which the linearization of the failure function is performed at a point called the Most Probable Point (MPP). In this method, the performance function is first defined as [22]: where X is a random variable, Z and S represent resistance and load, respectively. Failure occurs when ( ) ( ) Z X S X , or ( ) 0 g X . This method starts with the transformation of original variables to standard normal variables with zero mean and unit variance according to Rosenblatt conversion. The goal is to find the MMP, the point that has the least distance from the failure surface to the origin in the normal space. The least distance between the limit state surface and the origin in the normal space is defined as the reliability index ( ). As a result, reliability is: where is a cumulative distribution function of standard normal variable, while Pf and R indicate the probability of failure and reliability, respectively.

Second Order Reliability Method (SORM)
For linear limit state functions, the FORM solution is accurate. For non-linear limit state functions, it is generally difficult to accurately calculate the probability of failure or reliability in mathematical calculations. Therefore, for non-linear functions, the use of the SORM is recommended. In this method, the second-order of Taylor expansion is used to approximate the probability of performance function failure [22]: where ki represents the ith main curve of the performance function.

Laminate damage theory
In composite material structures, there are many failure criteria for modeling structural failure under static plate stress conditions, each of which is not fully exact. In general, predictive failure models of a composite laminate includes micromechanical and macromechanical models. Micromechanical models use the fiber-matrix level to predict layer failure. On the other hand, macro-mechanical models often rely on the description of medium failure without considering the fracture details at the fiber-matrix level. The damages at the microscopic state are quite complex and intricate to use to assess a structure. While the macro-mechanical model is the same as the CDM approach, it is still a common method.
A composite laminate may be modeled as a system with a set of subsystems (plies), each determined by a specific limit state function. The laminate failure may be expressed by either the failure of the first plie (FPF) or the failure of all plies in the laminate, failure of the last plie (LPF). In the FPF conditions, the goal is to prevent any micro-cracks in the material. On the other hand, in the LPF, the purpose is to use the laminate in the best feasible way and optimize the laminate application.

Classification of laminate material damage mechanisms
Damage modes of composite laminates materials can be mainly divided into fiber breakage, matrix cracking, fiber/matrix debonding, and delamination. Each damage mode occurs when the stress or strain component in the principal material coordinates attains a critical value. In this study, two failure modes, i.e., matrix cracking, fiber/matrix debonding will be discussed, and fiber breakage and delamination are ignored.
In matrix cracking, the growth of microcracks in a ply subjected to transverse stress is under unstable load control, and the sudden failure starts with macrocracks along the fiber/matrix interface. However, if the 90degree ply is placed in a multidirectional laminate, its failure is no longer sudden under transverse stress, and transverse stiffness gradually decreases with strain. Debonding of the fiber/matrix interface is mainly due to in-plane shear stress. As the interface debonds, the discontinuous surfaces' relative slip causes the plastic behavior in the shear direction.

Continuum damage mechanics (CDM)
In the CDM approach, an irreversible thermodynamic theory is used as a logical base to formulate constitutive equations with damage. In this model, the free energy, which indicates the relationship between the variables of the internal variable and their conjugate forces, and the dissipation function, which describes the development of internal variables, must first be expressed.
This work's CDM approach is an elementary ply model suggested by Ladeve`ze and Le Dantec [7]. At the ply state, the damage is supposed to be in the form of matrix cracking, fiber/matrix debonding, and fiber breakage. Fiber breakage is shown as an elastic-brittle failure, while matrix cracking and fiber/matrix interface debonding are considered as ply degradation, which causes damage growth in the ply. This model describes damage by reducing the stiffness of the material.
In the following, the concept of free energy and the dissipation function will be explained. The formulas are expressed in the materials' principal coordinates so that 1 and 2 show the fiber direction and perpendicular direction to the fibers, respectively.

Elastic constitutive equation
The Young's modulus of the undamaged ply in the fiber direction and perpendicular direction is expressed by 0 11 E .
The in-plane shear modulus is represented by 0 12 G . In G , which respectively, indicates the fracture in fiber, microcracking in the matrix and the debonding of fiber/matrix in the ply, then the Gibbs free energy for the ply in damaged state is defined as [7,14] : The elastic constitutive equation is [7,14] : For the assumed elastic fibers here, D = 0, so the damage development is controlled by two conjugate forces Y2 and Y12 calculated from the partial derivatives of the Gibbs free energy, respectively, relative to D2 and D12 [7,14]: It should be noted that the coercive forces of damage Y2 and Y12are called the energy density release rate functions, which are equivalent to the strain energy release rate G in fracture mechanics.

Plastic constitutive equations
It is supposed that the fibers are elastic. Therefore, the ply's plastic deformation is affected by the formation of matrix cracking and the matrix/fiber interface debonding. Assuming isotropic hardening of the material, the yield function may be written as [7,14] : where a, , and denote material constants.By specifying the accumulated plastic strain rate as [7,14] (11) and using F p , the plastic constitutive equation can be obtained as follows [7,14] : 11

Damage evolution equations
Variation in the elastic modulus in the experimental curves indicates that a reduction in material stiffness characterizes the damage. As mentioned, the ply damage is caused by matrix cracking and the matrix/fiber interface debonding. These states of damage are shown by the damage variables D2 and D12. The equivalent damage variable Y (Y2, Y12) and the damage criterion are [14]: 12 2 where b is the material parameter. Introducing

Methodology
To investigate the probabilistic behavior under inplane static loading, as shown in Figure 1, a rectangular glass/epoxy composite plate containing a central hole with [+45, -45 ]2Slamination sequence is discussed. The stress components in the principal coordinate system are: 11 22 12 1 2 (15) It is assumed that the normal strain in the fiber direction is small compared to the perpendicular direction. So, the strain components in the principal coordinate system are [23]: The ply damage model discussed here was implemented using a subroutine in Abaqus commercial software. A subroutine was applied to characterized matrix cracking G . In this subroutine, Y2 and Y12 are calculated to see if it has reached the set threshold level. In the occurrence of any damage, the damage variables will be updated.
The input parameters of the materials are extracted from a set of characterization tests, and maximum damage of ply is governed by critical energy values Ys. Using the Equation (14), the final damage states could result from the transverse tension or in-plane shear stress or ply strain. At this stage, the ply is completely damaged, and the damage variables are 12 2 1 D D . To apply reliability assessment, a failure function appropriate to the composite laminate is needed. The damage criterion introduced in this theory is used to obtain damage in laminate based on the CDM model. Therefore, when an additional quantitative variable Ŷ exceeds the damage threshold S Y , the damage will occur. Knowing this, the failure function for this case is: According to CDM approach, the energy density release rate is obtained according to the matrix cracking and fiber/matrix interface debonding. The brittle-damage threshold can also be achieved with experimental results. Using Equations (17)  , , , Uncertainties about the properties of the material should be extracted from some experiments and statistical modeling. Table 1 shows the mean value and coefficient of variation of the random variable for the glass/epoxy laminate.

Results
In this study, system probability of failure of rectangular laminates containing a central circle hole with various diameters under uniaxial tensile loading has been studied. To investigate the accuracy of this method, the shear stress curve is plotted according to the shear strain for a [+45, -45 ]2S laminate and the result is compared with the experimental graph reported in [24]. For mentioned laminate, the specimen length is 250 mm, width 50 mm, thickness 1 mm, edge distance mm. The material properties and parameters used in the finite element calculations are listed in Table 1. Figure 2 shows the result of this comparison, which indicates the accuracy of this method.
The growth of matrix cracking (D2) and debonding of fiber/matrix interface (D12) using damage contours in Figure 3 and Figure4 are shown d 5 and 10 mm, respectively. In finite element model, symmetry boundary conditions were used to model one-quarter of the laminated plate and unit displacement is applied at free end of the specimen.
As can be seen, the location of the maximum damage is at the edge of the hole and propagates over time to the free edge.
In this section, based on the scatter of material properties, random input parameters of the problem are modeled to establish feasible records for the probability distribution function (PDF) and cumulative distribution function (CDF) from plate damage. Probabilistic development of matrix cracking D2 and debonding of fiber/matrix interface of laminate D12 with 5-and 10-mm diameters  with respect to PDF and CDF is shown in Figure 5 and Figure 6, respectively. Figure 5 shows that in both hole diameters of 5-and 10-mm, the peak of matrix cracking damage (D2) is approximately 30% higher than the fiber/matrix debonding (D12) peak. As shown in Figure 6, although the two failure modes have overlap, it can be said that the probability of failure of laminate due to matrix cracking is greater than the probability of failure because of interface debonding. As the hole diameter increases, the probability of failure increases. According to the above survey, the effect of material dispersion has a considerable effect on the maximum and minimum damage range. If the scatter in the input parameters is not considered, there is a possibility of underestimating the development of damage.
The probability of failure due to matrix cracking and debonding of the fiber/matrix interface on the composite laminate with hole diameters of 5 and 10 mm for the +45 ply, -45 ply, series and parallel systems are shown in Figure 7. As shown in the figure, the failure probability of the +45 ply and -45 ply is equal. Details ofFigureFigure 7 using FORM and SORM are given in Table 10. In the case of series system failure, because of the growth of both the microcrack in the matrix and the debonding of the fiber/matrix interface in plies 1 to 8, according to FPF model, failure of only one ply leads to system failure. In a parallel combination, according to LPF model, if all its subsystems fail, the system fails. Therefore, if all plies (1-8) fail, the whole laminate will be failed for the studied laminate system.   It is obvious that much is still needed in order to arrive at a comprehensive probabilistic design tool for composite structures that encompasses the full simulation circle from the micromechanics level to the level of the complete structures. Several methods can be used for the estimation of the probability of failure of structures. In this study, the FORM and SORM are used to assess reliability analysis. Traditionally, the comparison of results from the developed stochastic methods is performed against Monte Carlo Simulation (MCS) data. But in the first place, MCS is too expensive for use in the complete design cycle, and in the second place, the alternative probabilistic method and the MCS both share the same modelling assumption.
It is worth noting that the results obtained are only for glass/epoxy laminate with [+45, -45 ]2S lamination sequence and with a central hole under tensile load and taking into account the two modes of failure, matrix cracking and fiber/matrix interface debonding. The results will be different for other materials, load, and geometry conditions. Investigating the probability of failure of carbon/epoxy laminate with different stacking sequences, considering other failure modes (in addition to matrix cracking and debonding of fiber/matrix interface) using the CDM approach and other probabilistic methods and models such as response surface method (RSM) are also recommended. It is also possible to obtain the acceptable failure probability of the laminate system by optimizing the central hole's diameter by considering the scattering of material properties and other factors using optimization methods. Theoreticalexperimental investigation of failure probability of laminated composite with complex geometry under fatigue loading by the current probability model is also suggested for future studies, which is ongoing research by authors.

Conclusion
The main emphasis of this research and development has been to develop a method based on the CDM approach that can be used to find the failure probability of composite laminates in complex structures. In this study, FEM, CDM approach, and reliability method are used to analyze the failure of composite structures containing circular hole with random variable properties under static tensile load. Studies show that by ignoring the scatter of random variables, underestimating the development of damage and failure in the structure is possible. The probability of failure of the laminate system by considering matrix cracking and debonding of matrix/fiber has been investigated, and it is shown that the failure probability of parallel system under problem conditions is within an acceptable range while failure probability of the series system is not within the acceptable range.