#### Abstract

A new 3D constitutive model for progressive damage analyses of unidirectional composite materials is presented, in which several important damage phenomena for the composite materials, such as the interfiber crack orientation, coupling of fiber failure and interfiber failure under longitudinal loads, closure effect for interfiber cracks, and longitudinal compressive behaviors under transversal constraints, have been considered comprehensively. A modified maximum stress failure criterion has been used for the damage onset prediction and a linear damage model has been adopted to establish the evolution rules of different damage. Numerical analyses with the model proposed have been implemented by using the subroutine UMAT in commercial software ABAQUS. Progressive damage analyses and static tensile experiments of a group of double-lap composite bolted joints have been carried out to validate the model proposed. Good agreements between the numerical and experimental results have been obtained.

#### 1. Introduction

With the increasing application of composite materials in aircrafts, automobiles, ships, and so forth, there is an urgent requirement for composite structure design methodologies. During the past decades, the structural design has mainly relied on experimental data by using the building block approach [1], which results in extremely expensive cost due to a large number of tests at each structural level. Thus, initiatives bridging the tests and analyses to reduce the design cost have become a common view in the world [2].

Great efforts have been taken to develop reliable and accurate analysis methods for composite structures. Empirical methods [3–6] obtained from plenty of tests and experiences are usually only suitable for certain structure types or load cases, though they are convenient and beneficial to the structure design. Besides, series of tests are required to determine some parameters for the empirical methods when the materials, ply sequences, or structural configurations change, which also lead to expensive cost. Theoretical analysis methods can significantly reduce the composite structural design and analysis costs and always draw lots of attention in the engineering. Among these methods, some predict the structural strength only with failure criteria. These methods will result in too conservative predictions since they ignore the structural nonlinear damage propagation process and fail to reveal the structural failure mechanism conveniently [2, 7]. In the later of the 20th century, a progressive damage method, which builds nonlinear mechanics models for composite materials and has capacity of accurately simulating the structural failure process from initial damage to ultimate failure, attracts widespread attention in composite structure analyses [8–11]. A progressive damage model contains three parts: a stress analysis model to obtain structural stress distribution, a failure criterion to estimate the structural damage and failure, and a material degradation model to control the property changes of damaged materials. Furthermore, the last one has been a hotspot in the composite research during decades since it directly represents the nonlinear damage propagation.

The material degradation models can fall into two categories, sudden degradation models and gradual degradation models. In the sudden degradation models, the material properties degrade to zero or a certain proportion of original values, which is called a degradation ratio. A total discounting model in which all the material properties degrade to zero when the damage is predicted is the simplest sudden degradation model. Based on the total discounting model, conservative predictions of structural strength are always obtained. Additionally, a usually used ply discounting model is one of extreme editions [12–15]. Another average practice is that the properties are degraded with targeted selection according to the failure modes. For example, Engelstad et al. [16] only reduced one material modulus at a time corresponding to a failure mode dominated by a stress component. More commonly, several material properties are reduced simultaneously for a failure mode, which is called interaction models [17–19]. Camanho and Matthews [10] distinguished the material degradation model by the tensile and compressive stress for the first time. Obviously, the sudden degradation models neglect the damage accumulation during the composite material failure process since they simply treat the material as undamaged and totally damaged.

To establish more accurate models, gradual degradation rules are good choices to describe the postfailure behavior of composite materials. The micromechanics theories are usually used for establishing gradual degradation models. Chang et al. [20–27] adopted the fiber bundle theory, while Lee et al. [28–30] used the shear lag theory to build the degradation model for the fiber failure. In contrast to the micromechanics methods, more researchers set the material property degradation rules with a predefined function as a matter of experience. In Lin and coworker’s [31, 32] gradual degradation models, the material degradation models were chosen to make the stress reduce linearly after the damage initiation. Hwang et al. [33] adopted an exponential function to establish the composite material degradation model.

The continuum damage mechanics is another way to propose the gradual damage model for composite materials, in which the damage variable is used to represent the damage scale in the materials. With the increasing damage variable from zero to unit, the material properties degrade from the original value to zero. In a strict sense, the evolution laws for damage variables should obey the thermodynamics principles of irreversible processes. Since the continuum damage mechanics has a rigorous theory basis, establishing suitable 3D continuum damage model for composite materials has been the major goal of the Third World-Wide Failure Exercise [34]. In addition, the continuum damage mechanics based gradual degradation models are more appealing. However, due to the complexity and intricacy of composite material failure mechanisms, the developments of continuum damage model for composite materials are slow. Existing continuum damage mechanics based gradual degradation models commonly focused on partial characters of the failure. Maimí et al. [2] and Raimondo et al. [35] proposed 2D gradual degradation models for composite laminates separately, in which the crack direction and crack closure were considered. Pinho et al. [36] developed a 3D gradual degradation model for composite materials. The fiber failure was focused on, and the crack direction was taken into account. However, the crack closure was not mentioned.

For cases where the composite structures surfer from 3D complex stresses such as bolted composite joints, the more accurate 3D gradual degradation model that can consider more damage phenomena of composite materials is still required. In this paper, a 3D gradual degradation model that takes account of the matrix crack direction, matrix crack closure, interaction of the fiber damage, and matrix damage as well as mechanical behaviors of compressive fiber failure is proposed based on the physical damage phenomena of composite materials. Furthermore, the progressive damage analysis is implemented for composite structures by using subroutine UMAT in ABAQUS and validated by the experimental results of double-lap composite bolted joints.

#### 2. 3D Progressive Damage Model Using Gradual Degradation Model

The degradation model of unidirectional composite materials is developed in this paper based on the characteristics of damage mechanisms of different failure modes. The modified maximum stress failure criteria [37] are used to predict the damage modes, initiation, and crack angles of composite materials. The linear damage evolution law is chosen to define the softening relationship for all failure modes. Finally, the model is embedded into the software ABAQUS with UMAT.

##### 2.1. Failure Mechanisms of Composite Materials

The composite materials have complex damage mechanisms in the mesoscale such as fiber breakage, fiber buckling, matrix cracking, interface debonding, and fiber bridging [2, 38]. The propagation and connection of damage in the mesoscale lead to the damage in the macroscale. According to the loads inducing the damage, the composite damage in the macroscale include five modes: fiber tensile damage, fiber compressive damage, matrix tensile damage, matrix compressive damage, and fiber-matrix shearing damage. In virtue of the regular fiber distribution in the mesoscale, some characteristic of the damage in the macroscale can be tracked, which is helpful to establish the mechanics model of composite materials.

Under the longitudinal loading, fiber damage initiates by isolated fiber fractures in weak zones [38]. The localized fractures will induce the matrix cracking, interface debonding, and matrix shearing failure in the adjoining fibers on account of local stress concentrations caused by the localized fractures, as shown in Figure 1. With the increment of the longitudinal loading, more fibers fail and further lead to the ultimate failure of the material. For the longitudinal tensile case, a macroscopic crack perpendicular to the fiber will be formed [2], while for the longitudinal compressive case a damaged band as a result of fiber buckling, fiber kinking, matrix crack, and interface debonding is formed [39, 40]. The macroscopic damaged band is assumed to be perpendicular to the fiber for simplification.

The matrix tensile failure, matrix compressive failure, and fiber-matrix shearing failure will form a crack parallel to the fiber in the macroscale, called an interfiber crack [41, 42]. Under transverse tensile or fiber-matrix shearing loadings, the propagation and connection of internal flaws such as interface debonding, rich resin zones, and microvoids result in the interfiber crack perpendicular to the middle plane of the ply [2, 41, 42], as shown in Figures 2(a) and 2(c). For the transverse compressive failure, the interfiber crack of most engineering composite materials usually has an angle of 53° inclined to the normal of the ply, as shown in Figure 2(b), rather than the maximum shear stress plane [41, 42]. Under the combination of a transverse compressive loading and a fiber-matrix shearing loading, the crack angle gradually reduces to 40° with the increase of the shearing loading but suddenly decreases to 0° when the shearing loads exceed a certain level [43], which is the same as that of a pure shearing crack. Generally, the interfiber crack for the matrix tensile failure, matrix compressive failure, and fiber-matrix shearing failure occurs in the matrix and fiber-matrix interface; as a result, it does not influence the properties in the fiber direction.

**(a)**

**(b)**

**(c)**

##### 2.2. Constitutive Model of Composite Materials for 3D Progressive Damage Analysis

The progressive damage model proposed in the paper groups the unidirectional composite material damage into two classes: fiber cracks and interfiber cracks. The fiber cracks include the fiber tensile crack and fiber compressive crack, in which the latter is not an actual crack but a kink band. The interfiber cracks include the matrix tensile crack, matrix compressive crack, and fiber-matrix shearing crack.

The progressive damage model treats the fiber crack as a plane of symmetry because the fiber crack shown in Figure 3 is generally perpendicular to fibers. The interfiber cracks and parallel to fibers will increase the anisotropy in the directions perpendicular to fibers. When an interfiber crack occurs, the material can be assumed orthotropic with and the plane perpendicular to as a symmetric plane. In addition, the interfiber cracks form a pair of free planes, on which only symmetric stresses exist and further result in symmetric damage. Generally, the symmetric damage leads to a new interfiber crack marked as , which is perpendicular to the interfiber crack and fiber crack . As the damaged material is orthotropic and the normal of interfiber cracks may change when the stresses are different, it is very hard to establish the constitutive relation in the global coordinate directly for the damaged material with arbitrary interfiber cracks. Thus, in current work, a constitutive model under a local coordinate O-1′2′3′ fixed to the cracks, as shown in Figure 3, is provided and further deduced in the global coordinate.

Based on the continuum damage mechanics, three damage variables , , and are adopted to establish the model: is corresponding to the fiber crack, to the first interfiber crack, and to the second interfiber crack. The constitutive model in the local coordinate O-1′2′3′ is expressed with (1). and are the local strain and stress vectors, respectively, given by (2) and (3). is the flexibility matrix. and () are the effective elastic moduli. () are the effective Poisson’s ratios. The model has considered the damage effects on the elastic moduli and Poisson’s ratios, the interaction between the fiber crack and interfiber cracks, the crack closure effect, and the longitudinal compressive behaviors when the materials are laterally constrained:

###### 2.2.1. Damage Effect on Material Elastic Moduli and Poisson’s Ratios

Longitudinal stresses will not produce a smooth crack but a damage band. There are not only fiber fracture (tension) and fiber buckling or kinking (compression) in the damage band, but also the matrix crack and fiber-matrix debonding. As a result, the fiber failure will decrease all the material properties. In contrast, interfiber damage seldom influences the load carrying capability in the fiber direction. To simplify the expression of the constitutive model, the effect of fiber failure on the transverse properties is considered by using the damage variables and in the interaction model of the fiber crack and interfiber crack. The effective elastic moduli and Poisson’s ratios can be expressed as

###### 2.2.2. Interaction Model of the Fiber Failure and Interfiber Failure under Longitudinal Loads

Longitudinal loadings can induce not only the fiber crack but also the interfiber cracks, while transverse and shearing stresses can only result in the interfiber cracks. Thus, to simplify the constitutive model, this paper assumes that when of the fiber crack is produced, the interfiber cracks equaling to occur in the materials, which can be expressed as follows:

###### 2.2.3. Closure Effect for Interfiber Cracks

The normal stress reversal from the tension to compression will make the interfiber crack close as shown in Figure 4. This paper assumes the closed crack will make the material recover the normal load-bearing capacity, which means the effective damage variable of the interfiber crack for the normal stress becomes zero. The effective damage variables of the interfiber cracks and for the normal stresses can be provided by (7). denotes the McCauley operator, defined by :

**(a)**

**(b)**

On the other hand, this paper assumes the closed cracks cannot recover the shear load-bearing capacity, even though the frictions between the crack surfaces exist. Thus, the effective damage variables of the interfiber cracks and for the shear stresses are equal to and , respectively, given by

###### 2.2.4. Consideration on the Longitudinal Stress Reversion

It is assumed that the effect of the fiber crack on the material properties under the longitudinal tensile stress is the same as that when it turns to compressive one. The effective damage variables and of the fiber crack during the longitudinal stress reversion are given by

###### 2.2.5. Model for the Longitudinal Compressive Behaviors considering Lateral Constraints

Under a longitudinal compressive loading, the composite material failure starts with the fiber kinking or buckling. With the load increases, the fiber kinking or buckling increases and induces the matrix cracks and fiber-matrix debonding. Consequently, the compressive load makes the fibers break into pieces [25].

The fully damaged composite materials have liquidity like sands [25, 26]. Thus, the behaviors of the composite materials under longitudinal compressive loads are related to the lateral constraint conditions. It is supposed that the longitudinal load-carrying capability of the composite materials will continually decrease up to be zero approximately owing to the outflowing of damaged materials without lateral constraints, as the curve I shown in Figure 5. And a higher longitudinal load can be withstood for damage materials with lateral constraints, as the curve II.

Due to different mechanical behaviors of damaged composite materials under longitudinal tensile and compressive loads, it is hard to use the same model to depict longitudinal compressive behaviors as that to characterize longitudinal tensile behaviors. Inspired by Chang’s model for the composite material longitudinal compressive behaviors [25], an approximately impressible constitutive model is utilized for the composite materials under longitudinal compressive loads in current work. As a result, the materials have a high load-carrying capability when lateral constraints are applied and a low load-carrying capability for case without lateral constraints by setting a relative low elastic modulus. Meanwhile, the damage variable is used to realize the transformation from the damage mechanics model to the approximately impressible model by controlling the elastic modulus and Poisson’s ratios. The elastic modulus of impressible model is set as same as the transverse elastic modulus of composite materials, . The longitudinal elastic modulus changes from to and all Poisson’s ratios change to 0.5 controlled by during the longitudinal compressive loading, given by (10). Meanwhile, a term of is added to last item of (10) to include the effect of the intercrack closure. In addition, a value of near to 0.5 is set to Poisson’s ratios in the impressible model since the impressible constitutive relation cannot be used in general finite elements. Here, is set to be 0.49:

Based on the assumptions above, the constitutive model of composite materials is established. The material parameters of the flexibility matrix in the model are listed in Table 1, from which the stiffness matrix in the local coordinate O-1′2′3′ can be further deduced. Then the stiffness matrix in global coordinate could be obtained from and the interfiber crack angle.

##### 2.3. Failure Criteria for Damage Onset

The modified maximum criterion proposed by Zhao et al. [37, 44] is adopted here to predict the onset of fiber crack and interfiber crack.

###### 2.3.1. Fiber Crack Onset

If either of the fiber tensile failure criterion or the fiber compressive failure criterion is satisfied, the fiber crack onset is detected. For fiber tensile () and compressive () failure, the criterion can be written as where and are the longitudinal tensile and compressive strength, respectively.

###### 2.3.2. Interfiber Crack Onset

If any criterion among the matrix tensile failure criterion, matrix compressive failure criterion, and fiber-matrix shearing failure criterion is satisfied, the interfiber crack onset is detected. The criterion is expressed in (12) and (13), where , , and are the transverse tensile strength, transverse compressive strength, and fiber-matrix shearing strength, respectively. , , and are the intercrack angles corresponding to the matrix tensile failure, matrix compressive failure, and fiber-matrix shearing failure, respectively.

Matrix tensile () and compressive () failure iswhere , .

Fiber-matrix shearing failure iswhere .

##### 2.4. Damage Evolution

After damage onset, the damage variable increases from 0 to unit with the strains increase. A linear damage model is used to define the damage evolution as shown in Figure 6, in which the stress decreases linearly with the strain increase after the damage onset. denotes the stress at damage onset. represents the final failure strain. The damage variable in the linear damage model is defined as

The damage variable of the fiber crack could be provided as where and relate to the onset and final failure strains of the fiber tensile damage, respectively. and are the onset and final failure strains of the fiber compressive damage, respectively.

For interfiber cracks, the normal stress and two shear stresses and may lead to the damage propagation. Assume the most critical stress drives the interfiber crack without regard to the stress interaction; then the damage variable is expressed as where .

For a complete continuum damage model, the damage onset strain and final failure strain are key parameters. The former is often obtained from the damage onset stress which is generally referred to the interface strength. Since the interface strength is difficult to be obtained from common tests, appropriate values are often used [45, 46]. This paper determines the damage onset strain by using the modified failure criteria and engineering material strengths. Moreover, the final failure strain could be provided with (17) according to the cohesive model proposed by Bažant and Oh [47], which applies the fracture toughness of crack to the damage mechanics model. In (17), is the fracture toughness, and is the characteristic length of the finite element. For the proposed damage mechanics model, the for fiber crack under tensile stress, for fiber crack under compressive stress, for interfiber crack under tensile stress, for interfiber crack under longitudinal shear stress, and for interfiber crack under transverse shear stress are needed. can be examined by double cantilever beam tests. and are obtained by three-point bending tests or four-point bending tests. As measuring is difficult, let be equal to in this paper. There is a lack of standard tests for and yet. Pinho et al. [48] used compact tension tests to examine and compact compression tests to determine :

#### 3. Validation of the Progressive Damage Model

To validate the capability of predicting damage propagation in the composite materials under complex loading conditions, the 3D progressive damage model was embedded into the software ABAQUS and applied to a double-lap composite bolted joint under a tensile load. Meanwhile, the tensile experiments of composite bolted joints were conducted to validate the numerical results. On considering the main object and the limited length of this paper, the implements of the 3D progressive damage model in ABAQUS will be described in a following independent paper.

##### 3.1. Static Tensile Experiment of Composite Bolted Joints

A double-lap composite bolted joint was selected to verify the model proposed here. Its configuration and dimensions were designed as recommended in engineering and are shown in Figure 7. The laminate strips were made of T800 carbon/epoxy composites with lay-ups of [45/0/−45/0/90/0/45/0/−45/0]_{s}. Two outer strips were connected to the central one by using a titanium alloy protruding head bolt HST13-6-7 and a high locking collar HST1078 made of A280 high temperature alloy. Composite filler was added adhesively at the gripping section of the two outer strips.

**(a)**

**(b)**

The experiments were carried out by using a 250 kN Instron8803 testing machine. Five specimens were tested. The displacement load with rate of 1 mm/min was applied up to the catastrophic failure of the joint. The applied load and grip holder displacements were measured by sensors of the testing machine automatically, while the hole deformations were examined by an extensometer mounted on the specimens. More information about the experimental process can be found in [7].

##### 3.2. Numerical Simulations

The FE model of the double-lap composite joint is shown in Figure 8, in which the composite laminates, fastener, and collar were established by using 3D element C3D8. Each layer of the laminates was established with one element through the thickness direction. Meshes around the hole of the composite laminates were refined to obtain accurate numerical predictions. Three meshes having 69260, 103100, and 150380 elements, named Model I, Model II, and Model III, respectively, were used to verify the mesh independence on the prediction results. Contact pairs were set in the potential contact regions which were between the outer strips and central strip, the hole and bolt, and the strips and bolt. The direct “hard” method was applied to calculate the normal behaviors. A model without friction for tangential behaviors was used to make the numerical simulation converge more easily. A preload of 4.7 kN calculated from the assembly tightening moment of 4.5 N-m by using empirical formulae [49, 50] was applied to the internal section of the bolt shank with a bolt load method in ABAQUS. Equivalent boundary constraints and displacement load shown in Figure 8 were established to decrease the element amounts. The ends of outer strips were fully constrained, while a displacement along the longitudinal direction was applied to the other end with transverse degrees constrained.

The elastic properties and strengths of the composite material are GPa, GPa, GPa. , MPa, MPa, MPa, MPa, and MPa. The Poisson ratio is given by using Christensen’s [51] equation: . is calculated from the transverse compressive strength using the equation [36]. is the transverse compressive crack angle, usually 53°. and were tested as 0.368 and 1.480 N/mm, respectively. has the same value as . and were set as 55.00 and 35.00 N/mm after tentative analysis within the reported range [35, 36, 48]. The elastic modulus and Poisson’s ratio of the titanium alloy are 110 GPa and 0.29, while being 200 GPa and 0.30 with respect to the A280 high temperature alloy, respectively.

The gradual degradation model and failure criteria were applied to the FE model by developing a user defined subroutine UMAT for ABAQUS. At the beginning of each load increment, the UMAT will be called. The damage state will be checked using failure criteria. Once material damage is predicted, the gradual degradation model will be used to determine each damage variable according to the stress conditions and give the stiffness matrix of damaged material. This operation is repeated time and time again until the structure is fully destroyed.

##### 3.3. Results and Discussion

Figure 9(a) gives the predicted load-hole deformation curves of double-lap composite joints by using the models having different meshes. It shows that the three curves are almost the same below 15 kN and slightly separate after that. With the element number increases, the load-hole deformation curve is higher before they reach each peak point. However, the differences between three curves are relatively small. The ultimate loads of three models are 19.78, 20.50, and 20.49 kN, respectively. The numerical predictions indicate that adopting model II will slightly affect the prediction of the damage process but have safely ignorable influence on the ultimate failure loads. Further, Figure 9(b) compares the numerical load-hole deformation of the double-lap composite joint by using Model II and experimental ones. The numerical curve has the similar trend as the experimental results. Both of them indicate that the double-lap composite joints exhibit large nonlinear behaviors after the elastic stage, although the predicted final failure deformation is smaller than the experimental ones. The predicted ultimate load from Model II is 1.5% lower than experimental average load 20.82 kN, as listed in Table 2. It follows that the proposed progressive damage model has successfully predicted the stiffness and strength of the composite bolted joint under tensile loads.

**(a)**

**(b)**

Figure 10(a) shows the experimental failure morphology of the central plate in the double-lap composite bolted joint. Large bearing damage occurred near the central plate hole. The materials were compressed into fragments and plumped up in the bearing damage zone. Close to the central plate end, shear-out damage mode can be found. The inner plies were squeezed out of the plate end in the zone. Figure 10(b) shows the predicted failure morphology of the central plate. Large hole deformation and bearing damage were predicted. The laminates bulge from the cut view indicates that the model has successfully simulated the bearing behavior of damaged composite materials. The model also shows similar damage morphology with the experimental one at the right end of central strip, where shear-out damage was formed with inner plies squeezed out.

**(a)**

**(b)**

Good agreements on the joint strength, stiffness, and failure morphology between the experimental outcomes and numerical predictions give evidence that the proposed PDM has great potential on simulating the mechanical behaviors of composite materials under complex loading and constraints.

#### 4. Conclusion

A new 3D constitutive model for progressive damage analysis of unidirectional composite materials is proposed in this paper. The model has taken account of the possible fiber fracture planes caused by longitudinal tensile and compressive loads as well as interfiber fracture planes at any angles with respect to the ply thickness direction caused by transverse stresses and fiber-matrix shear stresses. Mathematical models are also presented to describe complex phenomena of damaged materials, such as crack closure due to load reversion, interaction of fiber damage and interfiber damage, effective elastic moduli, and Poisson’s ratios of damaged materials. Besides, an approximate model has been created for the longitudinal compressive behaviors. The damage onset has been predicted by using the modified maximum stress failure criterion [37] and a linear damage model has been adopted to establish the evolution rules of different damage.

The proposed model has been embedded into ABAQUS using the subroutine UMAT and applied to predict the progressive damage accumulation and collapse failure of a double-lap composite bolted joint under tensile loads. The load-hole deformation curve, strength and stiffness information, progressive damage process, and final failure morphology can be obtained from the progressive damage simulation with the model proposed. The numerical prediction is in good agreement with the static tensile experimental outcomes of the composite bolted joint, which gives evidence of the effectiveness of the model proposed.

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgments

The research work is supported by the National High Technology Research and Development Program of China (2012AA040209) and the National Science Foundation of China (11372020, 11202012, 11472027, and 10902004).