Mixed-mode damage into a CGRP cruciform subjected to biaxial loading

In this paper, a three-dimensional progressive damage model (PDM) is implemented within a chopped glass-reinforced polyester (CGRP) cruciform structure for modelling its damage under loading. Three different cruciform specimens subjected to biaxial tensile loading are studied. In order to simulate the computational behaviour of the composite, the constitutive model considers an initial elastic behaviour followed by strain-softening. The initiation criterion defined is based on the maximum principal stress of the composite and once this criterion is satisfied, stiffness degradation starts. For the computation of damage, the influence of the fibre and the matrix are taken into account within the damage rule. Realistic values of the energy dissipated during damage are computed. The computational results obtained by means of an explicit time marching solver are compared with experimental outcomes for validation purposes. Finally, it is concluded that the PDM is able to localise the damage effectively as well as predicting its initiation. In the best of authors’ knowledge, this is the first time a three-dimensional PDM is implemented into a composite cruciform structure subjected to biaxial loading. 2015 The Authors. Published by Elsevier Ltd. This is an openaccess article under the CCBY license (http:// creativecommons.org/licenses/by/4.0/).


Introduction
In the last decade the amount of engineering applications using composite materials has increased significantly. This huge increment is attributed, mainly, to the high strength-weight ratio that they provide, for instance, in aircraft structures. Computational simulation of failure for composites has been an area of interest during many years because the testing costs can be notably reduced. For instance, according to Cox and Yang [5], in airframes on which the human lives depend currently 10,000 experimental tests of material specimens are required. Taking the previous example into account, it can be noted how important it is to obtain more reliable computational tools for this kind of materials in order to save and/or use the resources more efficiently.
In some cases, for validating material models, an uniaxial loading test is considered. This practice is not very often close to reality due to the fact that many applications using composite structures are subjected to multi-axial loading (biaxial or triaxial). Thus, more knowledge about the response of this kind of materials under different loading cases is required [17,18,12,25,21,20,6,26]. In particular, two principal approaches to deal with biaxial loading in composite materials can be distinguished: the first one using tubular specimens e.g. [24] and the second one by means of cruciform specimens e.g. [27,16]. In this work, a CGRP cruciform structure is studied for modelling damage. The cruciform structure develops a biaxial stress state in its central zone under loading of its arms. To address the damage behaviour of the structure, a new stress damage curve has been proposed for the in-plane loading under study. Special care has been taken into account in order to model a realistic energy of damage during the degradation process of the CGRP composite. As a result of the analysis, the PDM here considered is able to localise properly the damage produced in the area of damage initiation of the structure for each geometry and loading case considered. In the damage process, it is taken into account the influence of fibre rupture and matrix cracking. The validation of the numerical results is accomplished by comparison with experimental tests. To sum up, for a relatively complex loading problem the PDM predicts successfully the damage produced.
This paper is organised as follows. Section 2 describes briefly the PDM used for the fibre-reinforced composite studied. Section 3 shows the material behaviour and the numerical model proposed. Section 4 shows the implementation of the PDM in a CGRP cruciform specimen under biaxial tensile loading and its validation by comparison with experimental tests. Finally, Section 5 provides the conclusions of the work. In this work, the PDM proposed by Curiel-Sosa et al. [9] is modified in order to model damage in a cruciform structure. The main characteristic of the initial model [7] is that the paths for computation of damage offer an effective localisation of the several damage modes i.e. fibre rupture, matrix cracking, etc. In addition, comparing this PDM with others from literature it is noticed a notable different. The initial PDM considers interaction between damage modes. Hence, a certain damage mode, for instance fibre rupture, is affected by others and, so that, this characteristic makes the model more realistic. In reality, during the fracture process of composite structures, damage modes are not independent. For the simulation of the CGRP cruciform structure the initial model proposed is changed. The principal adaptations of the PDM to the problem here presented are: -Due to the quasi-isotropic material behaviour, a single damage variable is considered. This damage variable takes into account the influence of matrix cracking and fibre rupture for computing damage. -The initial mesh-dependency of the PDM has been alleviated by means of a relation between material properties. Hence, a realistic value of energy is dissipated during the damage process.
The initial model [9] is framed into strain-space damage models according to continuum damage mechanics. The thermodynamic framework of an irreversible process for dissipative materials [4] is considered as well as the idea of effective stressr [3] (Eq. (1)). The main mathematical formulation of the model is presented in the following equations: 1. Relation between effective stress and nominal stresŝ 2. Diagonal second-order tensor composed by damage variables 3. Effective stress taking into account the strain equivalence principlê 4. Stress-strain relation 5. Normalised energy release rates 6. Stress damage surfaces f n ðr; gÞ ¼ r T Á F n ðgÞ Á r À 1 n ¼ 1; 2; . . . ; m 7. Strain damage surfaces g n ð; gÞ ¼ T Á G n ðgÞ Á À 1 n ¼ 1; 2; . . . ; m 9. Damage directors v ð1Þ ¼ ½k ð1Þ 10. Definition of growth functions U n Firstly, for clarity in the notation, it has to be noted that bold characters denote tensor and vector variables. In Eq. (1), DðgÞ is the second-order diagonal tensor formed by the internal damage variables, g ij , and r ¼ ½r x ; r y ; r z ; r xy ; r yz ; r zx T is the nominal stress. These variables, g ij , are responsible of the stiffness degradation due to the different damage modes in the composite material: fibre rupture, fibre buckling or kinking, matrix cracking and matrix crushing.
The damage tensor D is built as a diagonal tensor and contains all the damage variables associated with each damage mode, Eq.
(2). The effective stress tensorr, which takes into account the strain equivalence principle [15], can be expressed as in Eq. (3). C 0 is the second-order constitutive tensor containing all stiffness components of the undamaged material and is the strain tensor.
The stress-strain relation is given in Eq. (4) where CðgÞ is the nonsymmetric damaged stiffness tensor.
In terms of stress, according to the plasticity theory, the first question that must be answered is how the stress state in the material is. For this purpose, stress damage surfaces associated with each damage mode n are defined. Then, the undamaged domain is delimited by the stress damage surfaces. These stress damage surfaces are built taking into account the so-called normalised energy release rates (NERR) given by Eq. (5). Every damage mode is characterised by a certain combination of these NERRs where E is the Young modulus, G is the in-plane shear modulus, X is the tensile strength and S the shear strength. In the composite here considered, i and j directions correspond with the x and y axis. The form of the stress damage surfaces is presented in Eq. (6). In this equation, F n ðgÞ is a second-order tensor associated to the damage mode n and m is the total number of damage modes modelled.
Once f n is obtained, the strain damage surfaces g n can be calculated by mapping into the strain space given by Eq. (7). For computing the damage, the damage variables give a value of the damage that is occurring in the structure due to the different damage modes. The time variation of damage variables is given by the damage rule expressed by Eq. (8). In that equation U n are growth functions and v n are the unitary damage directors. The modelling of the damage directors is made according to the degradation of stiffness due to a damage mode. For instance, fibre rupture v ð1Þ affects to the degradation in direction xx, xy and zx (see Eq. (9)). Thus, following the references [8,7] the weights k ij represent the influence of the damage modes in the computation of a certain damage mode. In previous research, these weights k ij were defined in a qualitative manner. In this work, it is proposed a new relation based on the volume fraction and the properties of both components of the material (fibre and matrix). By means of this new relation, a realistic amount of energy is dissipated during damage.
The mathematical definition of the growth functions U n for each damage mode n can be made according to Eq. (10). In this equation the growth functions U n are defined as the non-negative inner product between the strain gradient of the damage surface in the strain space, r Á g n , and the strain rate _ . Note that if the strain increment vector is pointing to the interior of the damage surface, for a determined damage mode, then no progression of such damage mode occurs.
A flowchart of the explicit time integration [2] is presented in Fig. 1 and the constitutive law subroutine in which the progressive damage model is implemented. A detailed description of the Finite Element Method (FEM) is out of the scope of this work, so just a general description is given to the reader for a better understanding of the whole procedure used. Interested readers in numerical integration or assembly can consult FEM references such as [13].

CGRP cruciform structure
The composite that forms the structure tested has a polymeric matrix with 20% volume of glass fibre reinforcement. The polyester matrix (a synolite resin) is reinforced with 4 layers of chopped strand mat. Due to the random distribution of the fibres this material presents a quasi-isotropic behaviour and its material properties are defined in the Table 1. Additionally, because of the manufacture process, the composite is homogeneous through the thickness, i.e. there is not an interface between layers like traditional laminates.
It is important to notice that the cruciform geometries under analysis are subjected to different biaxial loadings. Those loadings provoked the appearance of a macro-crack throughout the central zone of the cruciform specimen. Then, for geometry A the loading which will lead to failure through the diagonal of its central zone is 1/1, which means that the same load is applied in horizontal and vertical arms. In geometry B, the loading condition is 1.5/1, the first term means that more load is applied in horizontal arm (50% more) compared to the vertical one. Geometry C has a loading case 0.5/1, so that, half of the vertical load is applied in the horizontal arm.

Numerical model
The composite under analysis presents a quasi-isotropic behaviour and homogeneity through the thickness as has been demonstrated experimentally in [22,23]. Based on this fact, computationally, the composite is treated as an isotropic and homogeneous material (see Fig. 2).
In this analysis, an explicit central-difference time integration rule is considered with an automatic time increment. The numerical simulations have been performed by means of the finite element (FE) software ABAQUS [11]. For provoking failure in the  structures, a displacement in the tip of the arms is defined. This displacement is applied incrementally and quasi-statically. For the FE discretization, two different mesh regions are distinguished with different element size in the structure (see Fig. 3). The first one, it is in the central zone with a 0.4 mm size (zone of interest) and the second one defined in the arms with 1.4 mm size. The finite element used is an eight noded hexahedral element with reduced integration. Hourglassing control has been considered in order to avoid spurious deformation in the FE mesh. A nonlinear explicit dynamic analysis has been performed. A 3D model is studied, considering a 1/8 of the structure due to the symmetry of the problem and also in order to save computational costs considering symmetry boundary conditions.
A damage variable g is defined in order to model damage in the cruciform structure. This damage variable degrades equally the x and y components of the stress tensor i.e. r x and r y , due to the quasi-isotropic behaviour of the composite under analysis. The damage rate considers the influence of fibre rupture and matrix cracking as presented in Eq. (11): where _ g represents the damage variable rate. Note that, this variable is defined as a linear combination of damage growth for fibre rupture U f and matrix cracking U m . So, the increment of damage is a contribution of the fibre and matrix breakage. k f and k m quantify the influence of fibre rupture and matrix cracking respectively into the damage growth.
The definition of k f and k m is based on the rule of mixtures within the two components of the CGRP composite: fibre and matrix. The percentage in volume of the fibre is v fibre ¼ 20%, and its elastic modulus is E fibre ¼ 70 GPa. The remaining percentage of volume corresponds to the matrix (v matrix ¼ 80%) and its elastic modulus is E matrix ¼ 2 GPa. Tacking into account those material properties, the stiffness and the percentage of volume of the fibres is described as a function of the matrix properties in Eqs. (12) and (13).
Based on the above equations and considering the stiffness and the percentage of volume of each component separately it is possible to define two ratios denoted as a and b for fibre and matrix respectively. Those ratios are defined as: Notice that the term 3 8 in Eq. (14) has been considered to address the randomness of fibres throughout the CGRP composite [14]. Finally, the values of k f and k m are described as follows: By means of Eqs. (16) and (17), a more physical definition has been considered for k f and k m . As shown in Eq. (11), those scalar variables affect the computation of damage. In particular, k f and k m control the amount of energy dissipated during strain-softening i.e. the area under the curve stress-strain. The goodness of the proposed definition is verified in the next Section where a realistic value of energy is dissipated during simulations. The PDM here adopted is a local damage model in which damage variables depend on the strain state of the element under consideration and the numerical simulations exhibit mesh dependency. Therefore, it has been noticed that during mesh refinement the energy dissipated tends to zero when damage occurs. To overcome this limitation, some authors had considered different approaches to deal with mesh dependency. For example, regularisation techniques [10], nonlocal damage models (the damage variables depend on the strain state of the neighbourhood giving a characteristic length [19]) or the crack band model [1]. Using this last approach particularly, the strength limit is not kept constant in order to preserve the fracture energy constant. In this work, this strategy is not considered in order to preserve the material strength during simulations and therefore give a more realistic model. Hence, for this application, the strategy followed is based on preserve the energy dissipated during the fracture process in order to provide realistic results without changes in the material strength.
The energy released during an uniaxial test has been compared with the computational energy released in the rounded zone of the cruciform where the macro-crack is initiated. The value of energy estimated for the uniaxial test is G C I ¼ 6210 Pa [21], being mode I the dominant fracture mode. A fragile fracture process has been produced, considering that the strain at the moment of failure f is considered 1% higher that the yield strain Y ¼ 0:0138.
In the other hand, the expected strain-softening provided by simulations reached ultimate failure when the failure strain is 8% higher than yield strain Y . This deviation of energy released compared with mode I (uniaxial test) is attributed to the fact that in the rounded zone of the specimen a notable shear stress component is observed during loading [20] that consequently induces mode II of failure.
During simulations, an initiation criterion based on the maximum principal stress r max into CGRP cruciform is defined.
Consequently, damage is initiated once the maximum principal stress r max within the model reached the material strength i.e.
90 MPa. For the post-pick behaviour in the constitutive model, a damage variable is defined considering the influence of matrix cracking and fibre rupture. This consideration is realistic because the fibres are randomly embedded into the matrix so both damage processes occurred at the same time and aim failure.
A second-order tensor is defined into the PDM. The components of this second-order tensor are in-plane components. F damage corresponds to the damage tensor associated to the material degradation of the composite. By mapping those tensors into the stress space it is possible to define the stress damage curve, f damage ðr; gÞ ¼ r T Á F 1 ðgÞ Á r À 1. This stress curve defines the undamaged domain of the CGRP composite. Once damage grows this curve is degraded reducing the undamaged domain because of the material degradation. This behaviour is illustrated in Fig. 4. Looking at that figure, when the material is undamaged i.e. g ¼ 0, the undamaged domain is defined by the solid line and once the damage is developed (dashed line in Fig. 4), for instance g ¼ 0:6, the undamaged domain is reduced.

Results and validation
In this section, the results from simulations are compared with the experimental tests in order to validate the proposed numerical approach for a biaxial loading context. Several simulations ( Table 2) have been performed with three different geometries and loading conditions. Experimental results for geometry A are presented in Fig. 5. In that figure, the macro-crack is initiated in the rounded zone and propagated throughout the central zone. The same pattern of failure is observed during experiments for geometries B and C under its corresponding biaxial loading cases shown on Table 2. The CGRP composite owns random fibre distribution within the matrix so any crack observed will consequently provoke matrix cracking and fibre rupture. Due to the manufacture process, the composite is homogeneous through the thickness as In-plane stress space curves for the undamaged CGRP material, g = 0, and the damaged material, g = 0.6. Experimentally, for the cruciform specimen A, a loading case 1/ 1 yields to an adequate failure across the central zone and the crack is initiated in the rounded zone. Considering the last fact, the PDM was embedded in geometry A. Hence, in Fig. 6(b), it is illustrated the map of damage predicted by the numerical approach. In that Figure, a damage variable called SDV1 is defined. It is important to notice that this variable is ranged between 0 (undamaged material) and 1 (total failure) and represents the percentage of damage as a contribution of fibre and matrix breakage. The localisation of damage in geometry A is properly addressed by the PDM and higher values of damage are found in the rounded zone as depicted. In Fig. 6(a), the map of principal stresses is presented. This map of stresses corresponds with the initiation of degradation in the area of stress concentration.
It has to be highlighted that experimentally for a loading case different to 1/1 in geometry A failure appears in the arm with bigger load applied so non-adequate failure is achieved. This fact is predicted by the proposed technique and it is shown in Fig. 7(b). In that figure, higher values of damage are located in the vertical arm where double load is applied compared to the load applied in the horizontal arm.
The PDM was also implemented in geometries B and C. In these geometries the biaxial loading case considered are 1.5/1 and 0.5/1 respectively. In Figs. 8(b) and 9(b), the corresponding map of damage for geometries B and C are addressed. For both geometries, the localisation of the damage is properly localised when comparing with experimental tests. Hence, the crack is initiated in the      rounded zone where the PDM predicts higher damage for all geometries.

Conclusions
The implementation of a PDM within a cruciform structure has been addressed. Mixed damage modes have been considered being matrix cracking and fibre rupture the agents that provoke the macro-crack initiation in the rounded zone of the cruciform. Special care has been taken into account in order to solve the mesh dependency issue and provide realistic results during simulations. Thus, the PDM has been able to: -Accurate localise the damage compared with experimental tests. This localisation depends on the two damage modes (fibre rupture and matrix cracking) that are involved on the failure. -Address the initiation of damage based on a maximum principal stress criterion. Then, when the value of maximum principal stresses in the rounded zone was 90 MPa the damage process started. This fact was observed in experimental observations.
Based on the computational results obtained and the comparison with experimental tests, it is possible to conclude that the PDM successfully predicts damage initiation in complex loading cases such as biaxial. Additionally, during the damage process, a realistic amount of energy is dissipated. The strategy adopted takes material properties of the composite such as the stiffness/volume ratio.