# A Three-Dimensional Ply Failure Model for Composite Structures

Authors: Maurício V. Donadon; Sérgio Frascino M. de Almeida; Mariano A. Arbelo; Alfredo R. de Faria
Journal: International Journal of Aerospace Engineering
 (2009)

---

## Abstract

A fully 3D failure model to predict damage in composite structures subjected to multiaxial loading is presented in this paper. The formulation incorporates shear nonlinearities effects, irreversible strains, damage and strain rate effects by using a viscoplastic damageable constitutive law. The proposed formulation enables the prediction of failure initiation and failure propagation by combining stress-based, damage mechanics and fracture mechanics approaches within an unified energy based context. An objectivity algorithm has been embedded into the formulation to avoid problems associated with strain localization and mesh dependence. The proposed model has been implemented into ABAQUS/Explicit FE code within brick elements as a userdefined material model. Numerical predictions for standard uniaxial tests at element and coupon levels are
presented and discussed.


---

## Body

## 1. Introduction

Damage in composite structures is a very complex phenomenon which can occur through a different number of failure mechanisms such as fibre breakage, fibre buckling, matrix cracking, fibre-matrix debonding and delamination either combined or individually. The increasing computational resources have allowed reliable prediction of such phenomenon with a certain degree of accuracy by using the finite element method. However, there is still a lot of work to be done in this field to better understand the physics of failure in order to improve numerical failure models for composite materials. The damage modelling in composites can be broadly divided into four different approaches:(i)
failure criteria approach,(ii)
fracture mechanics approach,(iii)
plasticity approach,(iv)
damage mechanics approach.Failure criteria approaches were initially developed for unidirectional materials and restricted to the static regime. They are divided into two categories: interactive and noninteractive criteria [1]. Noninteractive criteria assume that the failure modes are decoupled and specific expressions are used to identify each failure mechanism. In the stress based criteria, for example, each and every one of the stresses in the principal material coordinates must be less than the respective strengths; otherwise, fracture is said to have occurred. In a similar fashion, the maximum strain failure criteria states that the failure occurs if one of the strains in the principal material coordinates exceed its respective failure strain.On the other hand, interactive criteria assume an interaction between two or more failure mechanisms and they describe the failure surface in the stress or strain space. Usually stress or strain polynomial expressions are used to describe the boundaries for the failure surface or envelope. Any point inside the envelope shows no failure in the material. Several interactive failure criteria can be found in the open literature, such as the Tsai and Wu [2] criterion, Tsai-Hill criterion, Hoffman criteria among others [1]. Nevertheless, the disadvantage in using such polynomial criteria is that they do not say anything about the damage mechanisms themselves, therefore modified versions have been used to distinguish between failure modes. Hashin [3] proposed a three-dimensional failure criteria for unidirectional composites. In his model four distinct failure modes associated with fibre failure in tension or compression and matrix cracking in tension or compression are modelled separately. Engblom and Havelka [4] proposed the use of a combination of the Hashin [3] and Lee [5] failure criteria. They used the Hashin criterion to detect in-plane failure and the Lee criteria to predict delaminations. The degradation was performed by reducing the stresses associated with each failure mechanism to zero in the composite material constitutive law.Shivakumar et al. [6] used the Tsai-Wu failure criterion and maximum stress criterion to model low-velocity impact damage in composites. In this approach the Tsai-Wu criterion was used in order to determine whether the damage has occurred and the maximum stress criteria was used to identify the failure mode. Good agreement was obtained between computed and experimental damage areas.F.-K. Chang and K.-Y. Chang  [7, 8] proposed a progressive in-plane damage model for predicting the residual strength of notched laminated composites. Three in-plane failure modes are considered: matrix cracking, fibre-matrix shearing and fibre breakage. Their model is currently available in the material model library of the LS-DYNA3D explicit finite element code. In their model the shear stress-shear strain relation is assumed to be nonlinear and an expression proposed by Hahn and Tsain [9, 10] is used to represent composite behaviour in shear. For fibre breakage and/or fibre-matrix shearing the degree of property degradation within the damaged area depends on the size of the damage predicted by the fibre failure criterion. F.-K. Chang and K.-Y. Chang  [7] proposed a property reduction model for fibre failure based on the micromechanics approach for fibre bundle failure. It is postulated that for fibre failure both Ey and νyx are reduced to zero, but Ex and the shear modulus Gxy degenerate according to the Weibull distribution.Choi et al. [11] investigated the low velocity impacts on composite plates using a line-nose impactor. They used the Chang-Chang failure criterion to detect the initiation of matrix cracking and delamination. Good agreement between experimental and numerical results was achieved.Davies  and Zhang [12] studied the low-velocity impact damage in carbon/epoxy composite plates impacted by an hemispherical impactor. The plates tested were quasi-isotropic having a symmetrical lay-up. Different impact energy levels, thickness, dimensions and boundary conditions were considered. The finite element modelling was carried out using plate elements.The Chang-Chang failure criteria were used for the in-plane damage predictions and in general good agreement was obtained between the computed force histories and experiments.The Chang-Chang failure criteria have shown a reasonable performance for in plane damage predictions in composite structures and therefore, their applications have been widely reported in the literature [13–18].The prediction of the onset of delamination depends on the interlaminar stress state and interlaminar strength of the laminate. Kim and Soni [19] used the distribution of interlaminar normal stress, and averaged that stress along a ply thickness distance in all cases they studied. They assumed that failure occured when the average of the normal interlaminar stress value over the fixed distance reached the interlaminar tensile strength.Jen et al. [20] developed a model based on boundary layer theory to predict initiation and propagation of delamination in a composite laminate containing a central circular hole. The Hashin-Rotem failure criterion was adopted in their model to predict the loading and location at which the initiation of delamination occurs.Brewer and Lagace [21] proposed a quadratic stress criterion for initiation of delamination using the approach suggested by Kim and Soni [19]. According to their criterion only out-of-plane stresses contribute to delamination and they assume that the predicted stresses should be independent of the sign of the interlaminar shear stress. A similar criteria was also proposed by Liu et al. [22] to predict delamination in composite laminates.The disadvantage in using stress based criteria for composite materials is that the scale effects relating to the length of cracks subject to the same stress field cannot be modelled correctly [13, 23]. In the failure criteria approach either the position and size of the cracks are unknown. For these reasons the fracture mechanics approach may be more attractive. Fracture mechanics considers the strain energy at the front of a crack of a known size and compares the energy with critical quantities such as critical strain energy release rate.The fracture mechanics approach has been used to predict compression after impact strength of composite laminates [24–27]. In such models the damaged area is replaced by an equivalent hole and the inelastic deformation associated with fibre microbuckling that develops near the hole edge is replaced with a equivalent crack loaded on its faces by a bridging traction which is linearly reduced with the crack closing displacement. The diameter of the hole is obtained from X-radiographs and/or ultrasonic C-scan images. The results showed good correlation between analytical and experimental values.Another potential application of the fracture mechanics approach is its indirect use to predict progressive delamination in composites [28–31]. In such models stress-displacement constitutive laws describe the interfacial material behaviour and fracture mechanics concepts are used. The area defined by the constitutive relationship is equal to the fracture energy or energy release rate and once the stresses have been reduced to zero the fracture energy has been consumed and the crack propagates. Linear and quadratic interaction relationships were assumed to describe the crack propagation in mixed-mode delamination. Comparisons were made with experimental and closed-form results and good agreement was obtained.However, the fracture mechanics approach cannot be easily incorporated into a progressive failure methodology because its application requires an initial flaw. A possible solution is to use ahybrid approach by using a stress or strain-based criterion for the failure initiation and a fracture mechanics approach for the failure propagation.The plasticity approach is suitable for composites that exhibit ductile behaviour such as boron/aluminium, graphite/PEEK and other thermoplastic composites. Vaziri et al. [32] proposed an orthotropic plane stress material model that combines the classical flow theory of plasticity with a failure criterion. In their work the material constitutive law is assumed to be elastic-plastic and it has two stages. The first stage is the post-yield and pre-failure where an orthotropic plasticity model is used to model the nonlinear material behaviour. The second stage is the postfailure where brittle or ductile failure modes start to occur. Favourable agreement is obtained between experiment and the model results.The damage mechanics approach has been investigated by many researchers in recent years and its application to damage modelling in composites has shown to be efficient. The method was originally developed by Kachanov [33] and Rabotnov [34] and it has the potential to predict different composite failure modes such as matrix cracking, fibre fracture and delamination.Ladeveze and Dantec [35] proposed  an in-plane model based on damage mechanics to predict matrix microcracking and fibre/matrix debonding in unidirectional composites. Two internal damage variables were used to degrade the ply material properties, one of them associated with the transverse modulus and another with the in-plane shear modulus. A linear elastic-damage behaviour was assumed for tensile and compressive stresses and a plasticity model was developed to account for the inelastic strains in shear. Damage evolution laws associated with each failure mechanism were introduced which relate the damage variables to strain energy release rates in the ply. Tension, compression and cyclic shear tests were performed to determine the constants required in the damage-development laws. Comparisons with experiments were performed by the authors and a good correlation between numerical and experimental results was obtained.Johnson [36] applied the model suggested by Ladeveze and Dantec for the prediction of the in-plane damage response of fibre reinforced composite structures during crash and impact events. The damage model was implemented for shell elements into the PAM-CRASH explicit finite element code. An experimental programme was carried out in order to validate the model. Low velocity impact tests were performed using a drop test rig. The plates were simply supported on a square steel frame and they were impacted at their centre using a hemispherical impactor with 50 mm diameter. The mass of the impactor was 21 kg. The plates were fabricated from 16 plies of carbon/epoxy with a quasi-isotropic layup. Different energy levels were considered to give a range of different failure modes from rebound to full penetration. Based on force history comparisons between experiments and numerical results the model overpredicted the peak load. According to Johnson, the discrepancy between measured and predicted peak load is explained by the neglected delamination in the model. This fact was demonstrated by the author in his subsequent work [37] where the delamination was included in the damage modelling using contact interface conditions.Williams and Vaziri [38] implemented an in-plane damage model based on continuum damage mechanics (CDM) into LS-DYNA3D for impact damage simulation in composite laminates. The model was originally developed by Matzenmiller et al. [39]. The model considers three damage parameters: fibre failure damage parameter, matrix failure damage parameter and an extra damage parameter to account for the effect of damage in shear response. Individual stress based criteria for each failure mechanism are used for the damage initiation. The failure criterion defines certain regions in stress (or strain) space where the damage state does not change. The damage growth law adopted by Matzenmiller is a function of the strain and it assumes an exponential form. The stress/strain curve predicted by this damage function is a Weibull distribution that can be derived from statistical analysis of the probability of the failure of a bundle of fibres with initial defects. Impact simulations with different energy levels were performed and the performance of the proposed model was checked against the Chang-Chang failure criteria and experimental results. Some limitations of the model were pointed out by the authors. Firstly, the response is predicted by a single equation and, as a result the loading and postfailure responses cannot be separated, thus restricting the versatility of the model. Also, the dependence of the damage growth on the loading rate as well as mesh size.In his recent work Williams et al.  [40] proposed a plane-stress continuum damage mechanics model for composite materials. The model was implemented for shell elements into LSDYNA3D explicit finite element code. Also, the model is an extension of their previous work [41] and it deals with some issues related to the limitation of the model suggested by Matzenmiller et al. [39]. The approach adopted by the authors uses the sub-structuring concept where one integration point is used for each sub-laminate and composite laminate theory is used to obtain its effective material properties. The model assumes a bilinear damage growth law and the damage process has two phases, one associated with matrix/delamination damage and another with fibre breakage. The driving force for damage growth is assumed to be a strain potential function and the threshold values for each damage phase are experimentally determined. High-velocity and low-velocity impact simulations were performed in order to assess the performance of the model. Also, the results were compared with the model proposed by Matzenmiller et al. [39] and the existing Chang-Chang failure criteria. Pinho et al. [42] proposed a three-dimensional failure model to predict failure in composites. Their model is an extension of the plane-stress failure model proposed by Davila and Camanho [43] and it accounts for shear nonlinearities effects. The model was implemented into LS-DYNA3D finite element code and the authors obtained good correlation between numerical predictions and experimental results for static standard in-plane tests. However, neither strain rate effects nor strain localization problems associated with irregular mapped meshes were addressed by the authors.This paper presents a formulation for a three-dimensional ply failure model for composite laminates. The proposed model is an extension of the authors previous work [44]. The new version of the model incorporates strain rate effects in shear by means of a semiempirical viscoplastic constitutive law and combines a quadratic stress based criterion with the Mohr-Coulomb failure criterion to predict interfibre-failure (IFF) without knowing a priori the orientation of the fracture plane. The new version of the model also uses the Second Piola-Kirchhoff stress tensor which potentially avoids material deformation problems. High-order damage evolution laws are proposed to avoid both numerical instabilities and artificial stress waves propagation effects commonly observed in the numerical response of Finite Element Codes based on Explicit Time Integration Schemes. These laws compared to the widely used bilinear law ensure smoothness at damage initiation and fully damaged stress onsets leading to a more stable numerical response. Iterative energy based criteria have been also added to the previous version [44] to better predict damage under multiaxial stress states. The model also incorporates a three dimensional objectivity algorithm [45] which accounts for orthotropic and crack directionality effects, making the model to be completely free of strain localization and mesh dependence problems.


## 2. Formulation

This section presents the failure model formulation. The formulation is based on the Continuum Damage Mechanics (CDM) approach and enables the control of the energy dissipation associated with each failure mode regardless of mesh refinement and fracture plane orientation by using a smeared cracking formulation. Internal thermodynamically irreversible damage variables were defined in order to quantify damage concentration associated with each possible failure mode and predict the gradual stiffness reduction during the fracture process. The material model has been implemented into ABAQUS explicit finite element code within brick elements as an user defined material model.


### 2.1.3-D Orthotropic Stress-Strain Relationship

The orthotropic material law that relates second Piola-Kirchhoff stressS to the Green-St. Venant strain e in the local material coordinate system is written as
(1)S=TtCTe,
where e=[e11,e22,e33,e12,e23,e31] is the Green-St. Venant strain vector and T is the transformation matrix given by
(2)T=[k2l20002kll2k2000-2kl001000000k-l0000lk0-klkl000k2-l2],
where k=cos(β) and l=sin(β). β is the inplane rotation angle around the Z-direction which defines the orientation of the material axes with respect to the reference coordinate system. The compliance matrix C-1 is defined in terms of the material axes as
(3)C-1=[1E11-ν21E22-ν31E33000-ν12E111E22-ν32E33000-ν13E11-ν23E221E330000001G120000001G230000001G31],
where the subscripts denote the material axes, that is,
(4)νij=νxi′xj′Eii=Exi′.
Since C is symmetric,
(5)νijEii=νjiEjj.
The Green-St. Venant strain tensor e is given by
(6)e=12(FtF-I),
where F is the deformation gradient and I is the identity matrix. The Cauchy stress tensor σij can be determined in terms of the second Piola-Kirchhoff stress tensor S as follows:
(7)σij=J-1FSFt,
where J is defined as the Jacobian determinant which is the determinant of the deformation gradient F. The present formulation will predict realistic material behaviour for finite displacements and rotation as long as the strains are small.


### 2.2. Degraded Stresses

The degraded or damaged stresses are defined as stresses transmitted across the damaged part of the cross-section in a representative volume element of the material. Based on the isotropic damage theory as originally proposed by Kachanov [33], an orthotropic relationship between local stresses acting on the damaged configuration can be written in terms of the local effective stresses in the undamaged configuration at ply level as follows:(8){σ1dσ2dσ3dτ12dτ23dτ13d}=[(1-d1(ϵ1))000000(1-dmt(ϵ̅mt))0000001000000(1-dmt(ϵ̅mt))(1-d12(γ12))000000(1-dmt(ϵ̅mt))0000001]{σ1σ2dmcσ3dmcτ12dmcτ23dmcτ13dmc},
where d1,dmt, and d12 are internal variables introduced to quantify the damage concentration within the Representative Volume Element (RVE) [46, 47]. A detailed description on the physical meaning and definition of these variables will be given in the following sections. The stress components with the superscript dmc are the stresses degraded locally on the action fracture plane. Details about the determination of the local action plane and the local degradation procedure is given in Section 2.4.The proposed orthotropic relationship between degraded and intact stresses ensures that the material stiffness matrix is positive defined during the degradation process. Moreover, it is physically based. By inspecting carefully (8) one can notice that transverse compression damage dmc degrades the through-the-thickness stresses σ3,τ23 and τ13.From the physical point of view, this relationship incorporates the the main features observed in the experiments, such as follows:(i)
It accounts for fibre damage effects due to tensile and compression loadings by means of the internal variabled1.(ii)
It provides the coupling between shear induced damage and matrix cracking during the degradation process using the internal damage variablesd12,dmtanddmc.(iii)
It accounts for damage coupling between the normal stressσ2,out-of-plane shear stress τ23 (transverse shear cracking, which leads to delaminations) and inplane shear stress τ12 for matrix cracking predictions.


### 2.3. Fibre Failure

The failure index to detect fibre failure in tension is given by(9)F1t(σ1)=σ1Xt≥1.In order to detect the catastrophic failure in compression related to the total instability of the fibres, the maximum stress criteria is used to detect damage initiation:(10)F1c(σ1)=|σ1|Xc≥1,
where Xt and Xc are the longitudinal strengths in tension and compression, respectively. When one of criteria given above is met, damage commences and grows according to the damage evolution law proposed by Donadon et al. [44]:
(11)d1(ϵ1)=d1t(ϵ1+)+d1c(ϵ1-)-d1t(ϵ1+)d1c(ϵ1-),
where d1t(ϵ1+) and d1c(ϵ1-) are the contributions of the irreversible damage due to fibre breakage in tension and fibre kinking in compression, respectively:
(12)d1t(ϵ1+)=1-ϵ1,0tϵ1+[1+κ1,t2(ϵ1+)(2κ1,t(ϵ1+)-3)],d1c(ϵ1-)=1-ϵ1,0cϵ1-[1+κ1,c2(ϵ1-)(2κ1,c(ϵ1-)-3)]
with
(13)κ1,t(ϵ1+)=ϵ1+-ϵ1,0tϵ1,ft-ϵ1,0t,κ1,c(ϵ1-)=ϵ1--ϵ1,0cϵ1,fc-ϵ1,0c,
where ϵ1,0t and ϵ1,0c are the failure strains in tension and compression, respectively. ϵ1+=max(ϵ1(t),ϵ1,0t) and ϵ1-=max(|ϵ1(t)|,ϵ1,0c) are the maximum achieved strains in the strain time history in tension (σ1>0) and compression (σ1<0), respectively. ϵ1,ft and ϵ1,fc are the final strains in tension and compression which are written as a function of the tensile fibre breakage and compression fibre kinking fracture toughnesses, respectively, as follows,
(14)ϵ1,ft=2GfibretXtl*,ϵ1,fc=2GfibrecXcl*
where Gfibret and Gfibrec are the intralaminar fracture toughnesses associated with fibre breakage in tension and compression, respectively and l* is characteristic element related to the size of the process zone. Experimental procedures and data reduction schemes to characterise the intralaminar fracture toughness for composites can be found in [48].The proposed damage evolution lawd1(ϵ1) results in a thermodynamically consistent degradation procedure, regularizing damage and combining stress, damage mechanics and fracture mechanics based approaches within an unified way similarly to the approach proposed by Miami et al. [49]


### 2.4. Interfibre Failure (IFF)

The interfibre failure modes consist of transverse matrix cracking either in tension or compression. Based on the Worldwide Failure Exercise experimental results [50] Pinho et al.  [42, 51] found that the failure envelope defined between the transverse stress σ2 and in-plane shear stress τ12 is accurately described by a quadratic interaction criteria. Thus, a failure index based on the interactive quadratic failure criterion given in [42, 51] has been used to predict tensile transverse matrix cracking. For tensile matrix cracking a failure index based on an interactive quadratic failure criterion written in terms of tensile and shear stresses is proposed in the following form:
(15)F2t(σ2,τ23,τ12)=(σ2Yt)2+(τ23S23)2+(τ12S12)2≥1.Once the criterion above is met, the proposed expression for damage growth due to tensile matrix cracking is given by(16)dmt(ϵ̅mt)=1-ϵ̅m,0tϵ̅mt[1+κm,t2(ϵ̅mt)(2κm,t(ϵ̅mt)-3)]
with
(17)κm,t(ϵ̅mt)=ϵ̅mt-ϵ̅m,0tϵ̅m,ft-ϵ̅m,0t,
where ϵ̅mt is the defined as resultant strain which is given by,
(18)ϵ̅mt=ϵ22+γ̅s,m2
with
(19)γ̅s,m=γ122+γ232,ϵ̅m,0t and σ̅m,0t are the damage onset resultant strain and stress, respectively, that is,
(20)ϵ̅m,0t=ϵ̅mt|F2t=1,σ̅m,0t=σ̅mt|F2t=1.In order to account for damage irreversibility effectsϵ̅mt=max(ϵ̅mt(t),ϵ̅m,0t) must be used in (16), where ϵ̅mt is the maximum achieved resultant strain in the strain time history.  The derivation of the resultant failure strain ϵ̅m,ft associated with tensile/shear matrix cracking is based on a power law criterion, which accounts for interactions between energies per unit of volume of damaged material within the RVE subjected to tensile and shear loadings. The power law energy criterion is given in the following form,
(21)(gmtgmct)λ+(gmsgmcs)λ=1
with λ=1 for UD laminates. The resultant stress in the transverse direction (or matrix direction) due to combined tensile and shear loadings is given by
(22)σ̅mt=σ22+τ122+τ232=σ22+(τ̅mt)2.The tensile and resultant shear stress components can be written in terms of the resultant stress as follows:(23)σ2t=σ̅mtcos(θ),τ̅t=σ̅mtsin(θ),
where θ is the angle defined between the resultant stress and tensile normal stress in the transverse direction, that is, θ=acos[max(0,σ2t)/σ̅mt]. In a similar way the tensile and resultant shear strain components are given as follows:
(24)ϵ2t=ϵ̅mtcos(θ),γ̅s,mt=ϵ̅mtsin(θ).For the damage evolution law given by(16) the specific fracture energies associated with tensile and shear stresses are respectively given by(25)gmt=∫0ϵ2fσ2tdϵ2=σ̅m,0tϵ̅m,ftcos2(θ)2,gms=∫0γ̅s,mfτ̅tdγ̅=σ̅m,0tϵ̅m,ftsin2(θ)2,
where the area under the stress-strain curves defined by the proposed polynomial damage evolution laws is identical to the one defined by the widely used bilinear softening law given in [42, 44]. Substituting (25) into (21) we obtain the following expression for the final strain due to the combined tensile and shear stress state,
(26)ϵ̅m,ft=2σ̅m,0t[(cos2(θ)gmct)λ+(sin2(θ)gmcs)λ]-1/λ
where gmct and gmcs are the critical specific fracture energies. These energies are related to the intralaminar fracture toughnesses. By using a smeared cracking formulation [52] and assuming that for UD laminates the values of intralaminar toughnesses associated with tensile matrix cracking and shear matrix cracking are comparable with mode I and mode II interlaminar fracture toughnesses, a relationship between specific critical fracture energies and intralaminar fracture toughnesses can be written as follows:
(27)gmct=GIcl*,gmcs=GIIcl*,
where l* is the characteristic length associated with the length of the process zone for each particular failure mode. A detailed description about the characteristic length calculation will be presented in the following section.The failure index to detect matrix cracking in compression failure is based on the criterion proposed by Puck and Schürmann [53, 54]. Their criterion is based on the Mohr-Coulomb theory and it enables the prediction of fracture planes for any given stress state related to Interfibre-Failure (IFF) Modes. This criterion is currently the state of the art to predict transverse compression response of composite laminates. The failure index to detect matrix cracking in compression based on the failure criterion proposed by Puck and Shürmann [53, 54] can be written as follows,
(28)F2c(σnt,σnl)=(σntS23A+μntσnn)2+(σnlS12+μnlσnn)2≥1,
here the subscripts n, l and t refer to the normal and tangential directions in respect to the fracture plane direction. S12 is the inplane shear strength and S23A is the transverse shear strength in the potential fracture plane (Action Plane), which is given by [55],
(29)S23A=Yc2(1-sin(ϕ)cos(ϕ))
with
(30)ϕ=2θf-900,
where Yc is the transverse compression strength. The fracture angle can be determined either experimentally or alternatively using (30), where θf maximises the failure criterion. Following the Mohr-Coulomb failure theory the friction coefficients can be determined as a function of the material friction angle as follows:
(31)μnt=tanϕ=tan(2θf-900).In the absence of experimental values an orthotropic relationship for the friction coefficients can be used [54],
(32)μntS23A=μnlS12.The stress components acting on the potential fracture plane are written in terms of angleθf which defines the orientation of the fracture plane in respect to through-the-thickness direction (direction-3 in the local material coordinate system):(33)σnn(θf)=σ2m2+σ3(1-m2)+2τ23mn,σnt(θf)=-σ2mn+σ3mn+τ23(2m2-1),σtt(θf)=σ2(1-m2)+σ3m2-2mnτ23,σnl(θf)=τ12m+τ13n,σlt(θf)=-τ12n+τ13m,
where m=cos(θf) and n=sin(θf). θf is the rotation around the local fibre direction (direction-1 in the local material coordinate system). It is clear that in order to apply (28) θf must be known. Many authors have defined θf≈530 based on experimental results for standard compression tests in UD laminates [42, 44, 51]. This is true and has been also confirmed by Puck and Schurmann [53, 54]. However, θf≈530 is only valid for uniaxial compression loading. This implies that θf changes for different stress states and alternatives are needed in order to handle such a problem. By examining (28) it is possible to see that the criterion has also the potential of predicting transverse intralaminar shear cracking for values of θf different from 530. The transverse intralaminar shear cracking is a very important failure mode because it leads to delamination between adjacent layers. In order tackle this problem we have used an iterative procedure to compute the fracture plane orientation for a given stress state. The procedure consists of incrementally varying θf within the interval [-900≤θf≤900] for a given stress state defined at ply level and check if the failure index for matrix cracking in compression being reached. Once the failure index is reached the local shear components σnl and σnt acting on the candidate fracture plane are degraded to zero according to the following damage evolution law:
(34)dmc(ϵ̅mc)=1-ϵ̅m,0cϵ̅mc[1+κm,c2(ϵ̅mc)(2κm,c(ϵ̅mc)-3)]
with
(35)κm,c(ϵ̅mc)=ϵ̅mc-ϵ̅m,0cϵ̅m,fc-ϵ̅m,0c,
where ϵ̅mc is the resultant shear strain on the action plane which is defined as
(36)ϵ̅mc=ϵnl2+ϵnt2,σ̅mc=σnl2+σnt2
with
(37)ϵnl(θf)=γ12m+γ13n,ϵnt(θf)=-ϵ2mn+ϵ3mn+γ23(2m2-1),
where ϵ̅m,0c and σ̅m,0c are the damage onset resultant strain and stress, respectively, that is,
(38)ϵ̅m,0c=ϵ̅mc|F2c=1,σ̅m,0c=σ̅mc|F2c=1.The resultant final strain for transverse compression failure is defined in terms of mode II interlaminar fracture toughness as follows:(39)ϵ̅m,fc=2Gmatrixcσ̅m,0cl*=2GIIcσ̅m,0cl*.In order to account for damage irreversibility effectsϵ̅mc=max(ϵ̅mc(t),ϵ̅m,0c) must be used in (34), where ϵ̅mc is the maximum achieved resultant strain on the action plane in the strain time history. After degrading the shear stresses acting on the potential fracture angle, the stresses are rotated back to the local material coordinate system using the following transformation:(40)σ2dmc=m2σnn-2mnσnt(1-dmc(ϵ̅mc))+σtt(1-m2),σ3dmc=σnn(1-m2)+2mnσnt(1-dmc(ϵ̅mc))+σttm2,τ12dmc=mσnl(1-dmc(ϵ̅mc))-nσlt,τ13dmc=mnσnn+(2m2-1)σnt(1-dmc(ϵ̅mc))-nmσtt,τ23dmc=nσnl(1-dmc(ϵ̅mc))+mσlt.


### 2.5. In-Plane Shear Failure

The observed behaviour of glass and carbon fibres laminates generally shows marked rate dependence in matrix-dominated shear failure modes and for this reason a rate dependent constitutive model has been used to model the in-plane shear behaviour. The constitutive model formulation is based on previous work carried out by Donadon et al. [44, 56] and it accounts for shear nonlineatities, irreversible strains and damage within the RVE. The stress-strain behaviour for in-plane shear failure is defined as follows,
(41)τ12=αG12γ12
with
(42)G12=G120+c1(e-c2γ12-1),
where G120 is the initial shear modulus and c1, c2 are material constants obtained from static/quasistatic in-plane shear tests. α is the strain-rate enhancement given by the following law,
(43)α=1+e(γ̇12/c3)
where c3 is another material constant obtained from dynamic in-plane shear tests. By decomposing the total shear-strain into inelastic γ12i and γ12e elastic components, the inelastic shear-strain can be written in terms of the elastic and total strain components as follows:
(44)γ12i=γ12-γ12e=γ12-τ12(γ12)G120.The failure index for in-plane shear failure is based on the maximum stress criterion and it is given by:(45)F12(τ12)=|τ12|S12≥1.The proposed damage evolution law for in-plane shear failure is given by(46)d12(γ12)=1-γ12,0-γ12,0iγ12-γ12,0i[1+κ122(γ12)(2κ12(γ12)-3)]
with
(47)κ12(γ12)=γ12-γ12,0-2γ12,0iγ12,0-γ12,0i-γ12,f,
where γ12,0 and γ12,0i are the total strain and total inelastic strain at failure (τ12=S12), respectively, that is,
(48)γ12,0=γ12|S12,γ12,0i=γ12i|S12,γ12,f is written in terms of the intralaminar toughness in shear:
(49)γ12,f=2GshearS12l*.In the absence of experimental results forGshear, is reasonable to assume Gshear=GIIc for UD plies, where GIIc  is the mode II interlaminar fracture toughness.


### 2.6. Objectivity Algorithm

The smeared formulation described in the previous sections relates the specific energy within a Representative Volume Element (RVE) with the fracture energy of the material for each particular failure mode. Since finite elements are volume based, mesh dependency problems will arise as a result of the mesh refinement. The correction of the postfailure softening slope according to the finite element size, as reported by Bazant [52] seems an attractive solution for the problem. However the approach has some limitations. Firstly, the crack growth direction must be parallel to one edge of the finite element, which is not the case for multidirectional composite laminates where layers can have arbitrary directions. Secondly, it cannot handle nonstructured meshes required in most of the complex finite element models with geometric discontinuities. In order to overcome such limitations and to ensure the objectivity of the model for generalized situations, a methodology originally developed by Oliver [57] has been used and extended to handle composite layers [58]. The dependence of the characteristic length on the fracture energy as well as its mathematical expression, were derived based in the work proposed by Oliver [57]. The method ensures a constant energy dissipation regardless of mesh refinement, crack growth direction and element topology so that, it is still applicable to nonstructured meshes.


#### 2.6.1. Crack Modelling in the Continuous Medium

Imagine a singular line in a two dimensional domain as a continuous material line, across which displacements are continuous but displacement gradients are discontinuous. The condition for a point belonging to a singular line with unit normaln at this point is that the determinant of the acoustic tensor in the n direction be zero, that is [57],
(50)det[niCijklnl]=0.
Nonpositive materials bifurcate, producing singular lines and the equation above permits their direction to be determined at each point. In the context of standard finite elements of C0 continuity, a singular line can be modelled only by the sides of the elements, these being the only points in the mesh where displacement gradient discontinuities can be obtained. However, a crack produces not only displacement gradient discontinuities but also displacement discontinuities. This latter kind of discontinuity cannot be modelled by a C0 finite element mesh for finite levels of discretization. However, a displacement discontinuity can be modeled as the limit of two parallel singular lines Γ- and Γ+ which tend to coincide with each other. The band delimited by these lines is known as singular band, and h is its width.By assuming an orthogonal curvilinear coordinate system (x′,y′) in the interior of the band, where y′ coordinates lines are parallel to the singular lines Γ- and Γ+, and x′ are the straight coordinates lines. Let u+(y′) and u-(y′) be the displacement vectors on Γ+ and Γ- the relative displacement vector can be written as,
(51)ω(y′)=u+(y′)-u-(y′)
as a vector representing the displacement “jump” between the two singular lines and
(52)δ(y′)=limh↦0ω(y′)=limh↦0(u+(y′)-u-(y′)),
If δ≠0, the singular band is modelling a discontinuous displacement field as the limit of a continuous one. This allows a crack to be idealised as a limit (with mesh refinement) of a band of finite elements where, by means of some numerical mechanisms, the condition δ≠0 is satisfied.


#### 2.6.2. Displacement and Traction Vectors in the Singular Band

Consider a singular band in the solid, with a width h according to Figure1.Analysis within the singular band.
(a)(b)Along a coordinate linex′, the displacement vector u can be expanded from its value in the Γ- line using Taylor's series as
(53)u(x′,y′)=u-(y′)+[∂u∂x′]-Δx′+O(h2),
and consenquently, for a point in the Γ+ line
(54)u+(y′)=u-(y′)+[∂u∂x′]-h(y′)+O(h2).
From the (51), (53), (54) we can write
(55)u(x′,y′)=u-(y′)+[∂u∂x′]-Δx′+O(h2),(56)u(x′,y′)=u-(y′)+Δx′hω(y′)+0(h2)≅u-(y′)+ϕ(x′,y′)ω(y′),
where ϕ is a function to be determined, which approximates Δx′/h when h↦0.From (51) and (56) it can be seen that
(57)ϕ=0⇒Γ-,ϕ=1⇒Γ+.
The equilibrium across the singular band will be enforced by assuming the traction vector t acting on the plane defined by the normal n:
(58)ti=σijnj,
which is constant in the x′ direction that is
(59)t(x′,y′)=t+(y′)=t-(y′).


#### 2.6.3. Energy Dissipation Within the Band

For a generic deformation process which takes place over a timeτ(0≤τ≤∞) the specific energy dissipation (energy per unit volume) within a closed domain Ω* (see Figure 1) is given by
(60)gf=∫0∞σij(x′,y′,τ)dϵij(x′,y′,τ)=∫0∞σijϵ̇ijdτ.For uniaxial deformation processgf would be, for a given point, the area under the stress-strain curve at that point. By taking the linearized geometric equations gf can be expressed as,(61)gf=∫0∞σij12(∂u̇i∂xj+∂u̇j∂xi)dτ=∫0∞σij∂u̇i∂xjdτ-∫0∞σij12(∂u̇i∂xj+∂u̇j∂xi)dτ.
The last integrand in (61) is zero, being the product of a symmetric and antisymmetric tensor, so that
(62)gf=∫0∞σij∂u̇i∂xjdτ=∫0∞∂∂xj(σiju̇i)dτ
where the Cauchy's equations for quasistatic processes and negligible body forces (∂σij/∂xi=0) have been considered. The total dissipated energy in the domain Ω* is
(63)W*=∫Ω*gfdΩ*=∫Ω*{∫0∞∂∂xj(σiju̇i)dτ}dΩ*.By applying the Gauss's theorem to (63) and using (58) we obtain,
(64)W*=∫Γ*{∫0∞σijnju̇idτ}dΓ*=∫Γ*{∫0∞tiu̇idτ}dΓ*Owing to the infitesimal width of the band, the curvilinear integral in (64) can be evaluated only on the lines Γ*+ and Γ*- (see Figure 1):
(65)W*=∫Γ*+∪Γ*-{∫0∞tiu̇idτ}dy′,
and taking into account (56) and (59) we obtain,
(66)W*=∫Γ*+∪Γ*-{∫0∞ti-(y′,τ)u̇i-(y′,τ)dτ}dy′+∫Γ*+∪Γ*-{∫0∞ti-(y′,τ)ϕ(x′,y′)ω̇(y′,τ)dτ}dy′.The first integral vanishes because the contributions onΓ*+ and Γ*- cancel each other out. Thus, the dissipated energy on Ω* is
(67)W*=∫Γ*+∪Γ*-{ϕ(x′,y′)∫0∞ti-(y′,τ)ω̇(y′,τ)dτ}dy′.Now, if the case whereΩ*=Ω is considered; that is, the whole band between points A and B in Figure 1, the total energy dissipated within the band between points A and B is
(68)W=∫ΩgfdΩ=∫yA′yB′{∫0∞ti-(y′,τ)ω̇(y′,τ)dτ}dy′.Equation (68) establishes that the energy dissipated within the idealized band can be written as a curvilinear integral along its length. The integrand of (68) represents the energy dissipated per unit of area, which in terms of fracture mechanics, is the fracture energy Gf:
(69)Gf(y′)=∫0∞ti-(y′,τ)ω̇(y′,τ)dτ.IfGf is assumed to be a material property independent of the spatial position of the point from (67) and (69) we obtain,
(70)W*=∫Γ*+∪Γ*-Gfϕ(x′,y′)dy′=Gf∫Γ*ϕ(x′,y′)dy′.By applying Green's theorem and taking into account (60), we can write the energy dissipated within the band as,
(71)W*=Gf∫Ω*∂ϕ∂x′dΩ*=∫Ω*gfdΩ*.The local form of (71) is
(72)gf=Gf∂ϕ∂x′=Gfl*,
where
(73)l*(x′,y′)=(∂ϕ∂x′)-1.The parameterl* plays the role of relating the specific energy (per unit of volume) and the fracture energy (per unit of area). l* is also identified as the characteristic length or crack band width used in existing cracking models [52]. For the unidimensional case and using the proposed Hermitian stress-strain softening law, the specific energy can be written as
(74)gf=σ0ϵf2,
where σ0 is the material strength and the degraded stress is given by
(75)σd=σ(1-d(ϵ)),
where the damage evolution law d(ϵ) is given in terms of the strain as follows:
(76)d(ϵ)=1-ϵ0ϵ[1+κ2(ϵ)(2κ(ϵ)-3)],
with
(77)κ(ϵ)=ϵ-ϵ0ϵf-ϵ0.Using  (74) the failure strain ϵf can be written in terms of the fracture energy and the characteristic length as it follows:
(78)ϵf=2Gfσ0l*.


#### 2.6.4. Determination of the Functionϕ and the Characteristic Length l*

In order to apply the theory presented in the previous section to the discretized medium, consider a mesh ofC0 continuous hexahedron solid finite elements (see Figure 2).Figure 2
Determination of the characteristic length for hexahedron elements.A set of cracked elements is determined by using failure criteria for detecting the crack initiation, the crack orientation depends on the fibre direction. The cracked plane is defined here as a normal vector which is parallel to the fibres for fibre failure and normal to the fibre direction for matrix failure. The algorithm described in this section for determining the characteristic length was proposed by Oliver [57] and it has the advantages of calculating the characteristic length for arbitrary crack directions and any finite element geometry.From (71) the function ϕ has to be continuous and derivable, satisfying (57).  A simple function defined in the isoparametric coordinates ξ and η which fulfils these requirements is
(79)ϕ(ξ,η)=∑i=1ncNi(ξ,η)ϕi,
where nc is the number of corner nodes of a virtual plane located at the midplane of the element (nc=4 for our case), Ni are the standard C0 shape functions of an element of nc virtual nodes in its midplane and ϕi is the value of the ϕ at corner i. If the crack location inside the element is known, ϕi takes the value +1 if the corner node i is ahead the crack, and 0 otherwise. The function defined by (79) fulfils the required condition of continuity within elements and takes the values +1 for the nodes on the boundaries ahead the crack and 0 for the nodes on the boundaries behind the crack (see Figure 3).Figure 3
Finite element band modelling.In general, however, the exact crack location is not known, and usually only some indication of the onset of cracking and the crack directions is obtained at the integration points. The following algorithm proposed by Oliver [57] has been used for determining the characteristic length at each integration point j, as shown in Figure 4,Figure 4
Computation ofϕ values at the virtual midplane of the element.(1)  A set of local cartesian axesx′, y′ is defined at the centre of the element, this being identified by the values of the isoparametric coordinates (ξ=0, η=0 and ζ=0). The direction of the local axis x′ is defined by the normal to the fracture plane, which is the fibre angle for fibre failure (θj=θf) and (θj=θf+900) for matrix failure. 										(80){xi′yi′}=[cos(θj)sin(θj)-sin(θj)cos(θj)]{xiyi}(2)  Values ofϕ at each corner node are established according to their position with respect to the local axis x′, y′ (ϕi=1 if xi′≥0, otherwise ϕi=0).(3)  The characteristic length, at the present integration pointj with isoparametric coordinates ξj and ηj and cracking angle θj is obtained as
(81)l*(ξj,ηj)=(∂ϕ(ξj,ηj)∂x′)-1=(∑i=1nc[∂Ni(ξj,ηj)∂xcos(θj)+∂Ni(ξj,ηj)∂ysin(θj)]ϕi)-1
where
(82){∂N∂x∂N∂y}=Jxy-1{∂N∂ξ∂N∂η},
and Jxy is defined as the Jacobian matrix given by
(83)Jxy=[∂x∂ξ∂y∂ξ∂x∂η∂y∂η],
where the partial derivatives with respect to the isoparametric coordinates are written as(84)∂x∂ξ=14(1+η)x1-14(1+η)x2-14(1-η)x3+14(1-η)x4,∂x∂η=14(1+ξ)x1+14(1-ξ)x2-14(1-ξ)x3-14(1+ξ)x4,∂y∂ξ=14(1+η)y1-14(1+η)y2-14(1-η)y3+14(1-η)y4,∂y∂η=14(1+ξ)y1+14(1-ξ)y2-14(1-ξ)y3-14(1+ξ)y4,
where the pairs (xi,yi) refers to the global coordinates of the virtual nodes defining the midplane of the element.For transverse compression failure a set of cracked elements is determined by using the stress based criterion defined by (28) and the cracked plane is defined by the fracture angle θf. The function ϕ is given b(85)ϕ(ξ,ζ)=∑i=1ncNi(ξ,ζ)ϕi,
where nc is the number of corner nodes of a virtual plane located at the midplane of the element according to Figure 5, Ni are the linear shape functions defined previously and nc virtual nodes defining the virtual cracking midplane and ϕi is the value of the ϕ at corner i.Figure 5
Computation of the characteristic length for transverse compression.x′, z′ is an auxiliar coordinate system defined at the centre of the element, this being identified by the values of the isoparametric coordinates (ξ=0, η=0 and ζ=0) with the direction of the local axis x′ is defined by the normal to the fracture plane,
(86){xi′zi′}=[cos(θj)sin(θj)-sin(θj)cos(θj)]{xizi}.The values ofϕ at each corner node are established according to their position with respect to the local axis x′,z′ (ϕi=1 if xi′≥0, otherwise ϕi=0) in a similar way as the one described previously.The characteristic length associated with transverse compression at the present integration pointj with isoparametric coordinates ξj and ζj and fracture angle θf is given by
(87)l*(ξj,ζj)=(∂ϕ(ξj,ζj)∂x′)-1=(∑i=1nc[∂Ni(ξj,ζj)∂xcos(θf)+∂Ni(ξj,ζj)∂zsin(θf)]ϕi)-1
where
(88){∂N∂x∂N∂z}=Jxz-1{∂N∂ξ∂N∂ζ},
and Jxz is defined as the Jacobian matrix given by
(89)Jxz=[∂x∂ξ∂z∂ξ∂x∂ζ∂z∂ζ],
where the partial derivatives with respect to the isoparametric coordinates are written as(90)∂x∂ξ=14(1+ζ)x1-14(1+ζ)x2-14(1-ζ)x3+14(1-ζ)x4,∂x∂ζ=14(1+ξ)x1+14(1-ξ)x2-14(1-ξ)x3-14(1+ξ)x4,∂z∂ξ=14(1+ζ)z1-14(1+ζ)z2-14(1-ζ)z3+14(1-ζ)z4,∂z∂ζ=14(1+ξ)z1+14(1-ξ)z2-14(1-ξ)z3-14(1+ξ)z4,
where the pairs (xi,zi) refer to the global coordinates of the virtual nodes defining the midplane of the element. In-plane shear cracking is strongly dependent on the fibre orientarion within the element therefore, the characteristic length associated with in-plane shear failure has been assumed to be the same as the one defined for fibre failure or failure in the warp direction. For out-of-plane shear failure modes, cracks are assumed to be smeared over the thickness of the element with a crack band defined between upper and lower faces of the element which is equivalent to assume θf=900 in (87),
(91)l*(ξj,ζj)=(∂ϕ(ξj,ζj)∂x′)-1=(∑i=1nc[∂Ni(ξj,ζj)∂z]ϕi)-1


## 2.1.3-D Orthotropic Stress-Strain Relationship

The orthotropic material law that relates second Piola-Kirchhoff stressS to the Green-St. Venant strain e in the local material coordinate system is written as
(1)S=TtCTe,
where e=[e11,e22,e33,e12,e23,e31] is the Green-St. Venant strain vector and T is the transformation matrix given by
(2)T=[k2l20002kll2k2000-2kl001000000k-l0000lk0-klkl000k2-l2],
where k=cos(β) and l=sin(β). β is the inplane rotation angle around the Z-direction which defines the orientation of the material axes with respect to the reference coordinate system. The compliance matrix C-1 is defined in terms of the material axes as
(3)C-1=[1E11-ν21E22-ν31E33000-ν12E111E22-ν32E33000-ν13E11-ν23E221E330000001G120000001G230000001G31],
where the subscripts denote the material axes, that is,
(4)νij=νxi′xj′Eii=Exi′.
Since C is symmetric,
(5)νijEii=νjiEjj.
The Green-St. Venant strain tensor e is given by
(6)e=12(FtF-I),
where F is the deformation gradient and I is the identity matrix. The Cauchy stress tensor σij can be determined in terms of the second Piola-Kirchhoff stress tensor S as follows:
(7)σij=J-1FSFt,
where J is defined as the Jacobian determinant which is the determinant of the deformation gradient F. The present formulation will predict realistic material behaviour for finite displacements and rotation as long as the strains are small.


## 2.2. Degraded Stresses

The degraded or damaged stresses are defined as stresses transmitted across the damaged part of the cross-section in a representative volume element of the material. Based on the isotropic damage theory as originally proposed by Kachanov [33], an orthotropic relationship between local stresses acting on the damaged configuration can be written in terms of the local effective stresses in the undamaged configuration at ply level as follows:(8){σ1dσ2dσ3dτ12dτ23dτ13d}=[(1-d1(ϵ1))000000(1-dmt(ϵ̅mt))0000001000000(1-dmt(ϵ̅mt))(1-d12(γ12))000000(1-dmt(ϵ̅mt))0000001]{σ1σ2dmcσ3dmcτ12dmcτ23dmcτ13dmc},
where d1,dmt, and d12 are internal variables introduced to quantify the damage concentration within the Representative Volume Element (RVE) [46, 47]. A detailed description on the physical meaning and definition of these variables will be given in the following sections. The stress components with the superscript dmc are the stresses degraded locally on the action fracture plane. Details about the determination of the local action plane and the local degradation procedure is given in Section 2.4.The proposed orthotropic relationship between degraded and intact stresses ensures that the material stiffness matrix is positive defined during the degradation process. Moreover, it is physically based. By inspecting carefully (8) one can notice that transverse compression damage dmc degrades the through-the-thickness stresses σ3,τ23 and τ13.From the physical point of view, this relationship incorporates the the main features observed in the experiments, such as follows:(i)
It accounts for fibre damage effects due to tensile and compression loadings by means of the internal variabled1.(ii)
It provides the coupling between shear induced damage and matrix cracking during the degradation process using the internal damage variablesd12,dmtanddmc.(iii)
It accounts for damage coupling between the normal stressσ2,out-of-plane shear stress τ23 (transverse shear cracking, which leads to delaminations) and inplane shear stress τ12 for matrix cracking predictions.


## 2.3. Fibre Failure

The failure index to detect fibre failure in tension is given by(9)F1t(σ1)=σ1Xt≥1.In order to detect the catastrophic failure in compression related to the total instability of the fibres, the maximum stress criteria is used to detect damage initiation:(10)F1c(σ1)=|σ1|Xc≥1,
where Xt and Xc are the longitudinal strengths in tension and compression, respectively. When one of criteria given above is met, damage commences and grows according to the damage evolution law proposed by Donadon et al. [44]:
(11)d1(ϵ1)=d1t(ϵ1+)+d1c(ϵ1-)-d1t(ϵ1+)d1c(ϵ1-),
where d1t(ϵ1+) and d1c(ϵ1-) are the contributions of the irreversible damage due to fibre breakage in tension and fibre kinking in compression, respectively:
(12)d1t(ϵ1+)=1-ϵ1,0tϵ1+[1+κ1,t2(ϵ1+)(2κ1,t(ϵ1+)-3)],d1c(ϵ1-)=1-ϵ1,0cϵ1-[1+κ1,c2(ϵ1-)(2κ1,c(ϵ1-)-3)]
with
(13)κ1,t(ϵ1+)=ϵ1+-ϵ1,0tϵ1,ft-ϵ1,0t,κ1,c(ϵ1-)=ϵ1--ϵ1,0cϵ1,fc-ϵ1,0c,
where ϵ1,0t and ϵ1,0c are the failure strains in tension and compression, respectively. ϵ1+=max(ϵ1(t),ϵ1,0t) and ϵ1-=max(|ϵ1(t)|,ϵ1,0c) are the maximum achieved strains in the strain time history in tension (σ1>0) and compression (σ1<0), respectively. ϵ1,ft and ϵ1,fc are the final strains in tension and compression which are written as a function of the tensile fibre breakage and compression fibre kinking fracture toughnesses, respectively, as follows,
(14)ϵ1,ft=2GfibretXtl*,ϵ1,fc=2GfibrecXcl*
where Gfibret and Gfibrec are the intralaminar fracture toughnesses associated with fibre breakage in tension and compression, respectively and l* is characteristic element related to the size of the process zone. Experimental procedures and data reduction schemes to characterise the intralaminar fracture toughness for composites can be found in [48].The proposed damage evolution lawd1(ϵ1) results in a thermodynamically consistent degradation procedure, regularizing damage and combining stress, damage mechanics and fracture mechanics based approaches within an unified way similarly to the approach proposed by Miami et al. [49]


## 2.4. Interfibre Failure (IFF)

The interfibre failure modes consist of transverse matrix cracking either in tension or compression. Based on the Worldwide Failure Exercise experimental results [50] Pinho et al.  [42, 51] found that the failure envelope defined between the transverse stress σ2 and in-plane shear stress τ12 is accurately described by a quadratic interaction criteria. Thus, a failure index based on the interactive quadratic failure criterion given in [42, 51] has been used to predict tensile transverse matrix cracking. For tensile matrix cracking a failure index based on an interactive quadratic failure criterion written in terms of tensile and shear stresses is proposed in the following form:
(15)F2t(σ2,τ23,τ12)=(σ2Yt)2+(τ23S23)2+(τ12S12)2≥1.Once the criterion above is met, the proposed expression for damage growth due to tensile matrix cracking is given by(16)dmt(ϵ̅mt)=1-ϵ̅m,0tϵ̅mt[1+κm,t2(ϵ̅mt)(2κm,t(ϵ̅mt)-3)]
with
(17)κm,t(ϵ̅mt)=ϵ̅mt-ϵ̅m,0tϵ̅m,ft-ϵ̅m,0t,
where ϵ̅mt is the defined as resultant strain which is given by,
(18)ϵ̅mt=ϵ22+γ̅s,m2
with
(19)γ̅s,m=γ122+γ232,ϵ̅m,0t and σ̅m,0t are the damage onset resultant strain and stress, respectively, that is,
(20)ϵ̅m,0t=ϵ̅mt|F2t=1,σ̅m,0t=σ̅mt|F2t=1.In order to account for damage irreversibility effectsϵ̅mt=max(ϵ̅mt(t),ϵ̅m,0t) must be used in (16), where ϵ̅mt is the maximum achieved resultant strain in the strain time history.  The derivation of the resultant failure strain ϵ̅m,ft associated with tensile/shear matrix cracking is based on a power law criterion, which accounts for interactions between energies per unit of volume of damaged material within the RVE subjected to tensile and shear loadings. The power law energy criterion is given in the following form,
(21)(gmtgmct)λ+(gmsgmcs)λ=1
with λ=1 for UD laminates. The resultant stress in the transverse direction (or matrix direction) due to combined tensile and shear loadings is given by
(22)σ̅mt=σ22+τ122+τ232=σ22+(τ̅mt)2.The tensile and resultant shear stress components can be written in terms of the resultant stress as follows:(23)σ2t=σ̅mtcos(θ),τ̅t=σ̅mtsin(θ),
where θ is the angle defined between the resultant stress and tensile normal stress in the transverse direction, that is, θ=acos[max(0,σ2t)/σ̅mt]. In a similar way the tensile and resultant shear strain components are given as follows:
(24)ϵ2t=ϵ̅mtcos(θ),γ̅s,mt=ϵ̅mtsin(θ).For the damage evolution law given by(16) the specific fracture energies associated with tensile and shear stresses are respectively given by(25)gmt=∫0ϵ2fσ2tdϵ2=σ̅m,0tϵ̅m,ftcos2(θ)2,gms=∫0γ̅s,mfτ̅tdγ̅=σ̅m,0tϵ̅m,ftsin2(θ)2,
where the area under the stress-strain curves defined by the proposed polynomial damage evolution laws is identical to the one defined by the widely used bilinear softening law given in [42, 44]. Substituting (25) into (21) we obtain the following expression for the final strain due to the combined tensile and shear stress state,
(26)ϵ̅m,ft=2σ̅m,0t[(cos2(θ)gmct)λ+(sin2(θ)gmcs)λ]-1/λ
where gmct and gmcs are the critical specific fracture energies. These energies are related to the intralaminar fracture toughnesses. By using a smeared cracking formulation [52] and assuming that for UD laminates the values of intralaminar toughnesses associated with tensile matrix cracking and shear matrix cracking are comparable with mode I and mode II interlaminar fracture toughnesses, a relationship between specific critical fracture energies and intralaminar fracture toughnesses can be written as follows:
(27)gmct=GIcl*,gmcs=GIIcl*,
where l* is the characteristic length associated with the length of the process zone for each particular failure mode. A detailed description about the characteristic length calculation will be presented in the following section.The failure index to detect matrix cracking in compression failure is based on the criterion proposed by Puck and Schürmann [53, 54]. Their criterion is based on the Mohr-Coulomb theory and it enables the prediction of fracture planes for any given stress state related to Interfibre-Failure (IFF) Modes. This criterion is currently the state of the art to predict transverse compression response of composite laminates. The failure index to detect matrix cracking in compression based on the failure criterion proposed by Puck and Shürmann [53, 54] can be written as follows,
(28)F2c(σnt,σnl)=(σntS23A+μntσnn)2+(σnlS12+μnlσnn)2≥1,
here the subscripts n, l and t refer to the normal and tangential directions in respect to the fracture plane direction. S12 is the inplane shear strength and S23A is the transverse shear strength in the potential fracture plane (Action Plane), which is given by [55],
(29)S23A=Yc2(1-sin(ϕ)cos(ϕ))
with
(30)ϕ=2θf-900,
where Yc is the transverse compression strength. The fracture angle can be determined either experimentally or alternatively using (30), where θf maximises the failure criterion. Following the Mohr-Coulomb failure theory the friction coefficients can be determined as a function of the material friction angle as follows:
(31)μnt=tanϕ=tan(2θf-900).In the absence of experimental values an orthotropic relationship for the friction coefficients can be used [54],
(32)μntS23A=μnlS12.The stress components acting on the potential fracture plane are written in terms of angleθf which defines the orientation of the fracture plane in respect to through-the-thickness direction (direction-3 in the local material coordinate system):(33)σnn(θf)=σ2m2+σ3(1-m2)+2τ23mn,σnt(θf)=-σ2mn+σ3mn+τ23(2m2-1),σtt(θf)=σ2(1-m2)+σ3m2-2mnτ23,σnl(θf)=τ12m+τ13n,σlt(θf)=-τ12n+τ13m,
where m=cos(θf) and n=sin(θf). θf is the rotation around the local fibre direction (direction-1 in the local material coordinate system). It is clear that in order to apply (28) θf must be known. Many authors have defined θf≈530 based on experimental results for standard compression tests in UD laminates [42, 44, 51]. This is true and has been also confirmed by Puck and Schurmann [53, 54]. However, θf≈530 is only valid for uniaxial compression loading. This implies that θf changes for different stress states and alternatives are needed in order to handle such a problem. By examining (28) it is possible to see that the criterion has also the potential of predicting transverse intralaminar shear cracking for values of θf different from 530. The transverse intralaminar shear cracking is a very important failure mode because it leads to delamination between adjacent layers. In order tackle this problem we have used an iterative procedure to compute the fracture plane orientation for a given stress state. The procedure consists of incrementally varying θf within the interval [-900≤θf≤900] for a given stress state defined at ply level and check if the failure index for matrix cracking in compression being reached. Once the failure index is reached the local shear components σnl and σnt acting on the candidate fracture plane are degraded to zero according to the following damage evolution law:
(34)dmc(ϵ̅mc)=1-ϵ̅m,0cϵ̅mc[1+κm,c2(ϵ̅mc)(2κm,c(ϵ̅mc)-3)]
with
(35)κm,c(ϵ̅mc)=ϵ̅mc-ϵ̅m,0cϵ̅m,fc-ϵ̅m,0c,
where ϵ̅mc is the resultant shear strain on the action plane which is defined as
(36)ϵ̅mc=ϵnl2+ϵnt2,σ̅mc=σnl2+σnt2
with
(37)ϵnl(θf)=γ12m+γ13n,ϵnt(θf)=-ϵ2mn+ϵ3mn+γ23(2m2-1),
where ϵ̅m,0c and σ̅m,0c are the damage onset resultant strain and stress, respectively, that is,
(38)ϵ̅m,0c=ϵ̅mc|F2c=1,σ̅m,0c=σ̅mc|F2c=1.The resultant final strain for transverse compression failure is defined in terms of mode II interlaminar fracture toughness as follows:(39)ϵ̅m,fc=2Gmatrixcσ̅m,0cl*=2GIIcσ̅m,0cl*.In order to account for damage irreversibility effectsϵ̅mc=max(ϵ̅mc(t),ϵ̅m,0c) must be used in (34), where ϵ̅mc is the maximum achieved resultant strain on the action plane in the strain time history. After degrading the shear stresses acting on the potential fracture angle, the stresses are rotated back to the local material coordinate system using the following transformation:(40)σ2dmc=m2σnn-2mnσnt(1-dmc(ϵ̅mc))+σtt(1-m2),σ3dmc=σnn(1-m2)+2mnσnt(1-dmc(ϵ̅mc))+σttm2,τ12dmc=mσnl(1-dmc(ϵ̅mc))-nσlt,τ13dmc=mnσnn+(2m2-1)σnt(1-dmc(ϵ̅mc))-nmσtt,τ23dmc=nσnl(1-dmc(ϵ̅mc))+mσlt.


## 2.5. In-Plane Shear Failure

The observed behaviour of glass and carbon fibres laminates generally shows marked rate dependence in matrix-dominated shear failure modes and for this reason a rate dependent constitutive model has been used to model the in-plane shear behaviour. The constitutive model formulation is based on previous work carried out by Donadon et al. [44, 56] and it accounts for shear nonlineatities, irreversible strains and damage within the RVE. The stress-strain behaviour for in-plane shear failure is defined as follows,
(41)τ12=αG12γ12
with
(42)G12=G120+c1(e-c2γ12-1),
where G120 is the initial shear modulus and c1, c2 are material constants obtained from static/quasistatic in-plane shear tests. α is the strain-rate enhancement given by the following law,
(43)α=1+e(γ̇12/c3)
where c3 is another material constant obtained from dynamic in-plane shear tests. By decomposing the total shear-strain into inelastic γ12i and γ12e elastic components, the inelastic shear-strain can be written in terms of the elastic and total strain components as follows:
(44)γ12i=γ12-γ12e=γ12-τ12(γ12)G120.The failure index for in-plane shear failure is based on the maximum stress criterion and it is given by:(45)F12(τ12)=|τ12|S12≥1.The proposed damage evolution law for in-plane shear failure is given by(46)d12(γ12)=1-γ12,0-γ12,0iγ12-γ12,0i[1+κ122(γ12)(2κ12(γ12)-3)]
with
(47)κ12(γ12)=γ12-γ12,0-2γ12,0iγ12,0-γ12,0i-γ12,f,
where γ12,0 and γ12,0i are the total strain and total inelastic strain at failure (τ12=S12), respectively, that is,
(48)γ12,0=γ12|S12,γ12,0i=γ12i|S12,γ12,f is written in terms of the intralaminar toughness in shear:
(49)γ12,f=2GshearS12l*.In the absence of experimental results forGshear, is reasonable to assume Gshear=GIIc for UD plies, where GIIc  is the mode II interlaminar fracture toughness.


## 2.6. Objectivity Algorithm

The smeared formulation described in the previous sections relates the specific energy within a Representative Volume Element (RVE) with the fracture energy of the material for each particular failure mode. Since finite elements are volume based, mesh dependency problems will arise as a result of the mesh refinement. The correction of the postfailure softening slope according to the finite element size, as reported by Bazant [52] seems an attractive solution for the problem. However the approach has some limitations. Firstly, the crack growth direction must be parallel to one edge of the finite element, which is not the case for multidirectional composite laminates where layers can have arbitrary directions. Secondly, it cannot handle nonstructured meshes required in most of the complex finite element models with geometric discontinuities. In order to overcome such limitations and to ensure the objectivity of the model for generalized situations, a methodology originally developed by Oliver [57] has been used and extended to handle composite layers [58]. The dependence of the characteristic length on the fracture energy as well as its mathematical expression, were derived based in the work proposed by Oliver [57]. The method ensures a constant energy dissipation regardless of mesh refinement, crack growth direction and element topology so that, it is still applicable to nonstructured meshes.


### 2.6.1. Crack Modelling in the Continuous Medium

Imagine a singular line in a two dimensional domain as a continuous material line, across which displacements are continuous but displacement gradients are discontinuous. The condition for a point belonging to a singular line with unit normaln at this point is that the determinant of the acoustic tensor in the n direction be zero, that is [57],
(50)det[niCijklnl]=0.
Nonpositive materials bifurcate, producing singular lines and the equation above permits their direction to be determined at each point. In the context of standard finite elements of C0 continuity, a singular line can be modelled only by the sides of the elements, these being the only points in the mesh where displacement gradient discontinuities can be obtained. However, a crack produces not only displacement gradient discontinuities but also displacement discontinuities. This latter kind of discontinuity cannot be modelled by a C0 finite element mesh for finite levels of discretization. However, a displacement discontinuity can be modeled as the limit of two parallel singular lines Γ- and Γ+ which tend to coincide with each other. The band delimited by these lines is known as singular band, and h is its width.By assuming an orthogonal curvilinear coordinate system (x′,y′) in the interior of the band, where y′ coordinates lines are parallel to the singular lines Γ- and Γ+, and x′ are the straight coordinates lines. Let u+(y′) and u-(y′) be the displacement vectors on Γ+ and Γ- the relative displacement vector can be written as,
(51)ω(y′)=u+(y′)-u-(y′)
as a vector representing the displacement “jump” between the two singular lines and
(52)δ(y′)=limh↦0ω(y′)=limh↦0(u+(y′)-u-(y′)),
If δ≠0, the singular band is modelling a discontinuous displacement field as the limit of a continuous one. This allows a crack to be idealised as a limit (with mesh refinement) of a band of finite elements where, by means of some numerical mechanisms, the condition δ≠0 is satisfied.


### 2.6.2. Displacement and Traction Vectors in the Singular Band

Consider a singular band in the solid, with a width h according to Figure1.Analysis within the singular band.
(a)(b)Along a coordinate linex′, the displacement vector u can be expanded from its value in the Γ- line using Taylor's series as
(53)u(x′,y′)=u-(y′)+[∂u∂x′]-Δx′+O(h2),
and consenquently, for a point in the Γ+ line
(54)u+(y′)=u-(y′)+[∂u∂x′]-h(y′)+O(h2).
From the (51), (53), (54) we can write
(55)u(x′,y′)=u-(y′)+[∂u∂x′]-Δx′+O(h2),(56)u(x′,y′)=u-(y′)+Δx′hω(y′)+0(h2)≅u-(y′)+ϕ(x′,y′)ω(y′),
where ϕ is a function to be determined, which approximates Δx′/h when h↦0.From (51) and (56) it can be seen that
(57)ϕ=0⇒Γ-,ϕ=1⇒Γ+.
The equilibrium across the singular band will be enforced by assuming the traction vector t acting on the plane defined by the normal n:
(58)ti=σijnj,
which is constant in the x′ direction that is
(59)t(x′,y′)=t+(y′)=t-(y′).


### 2.6.3. Energy Dissipation Within the Band

For a generic deformation process which takes place over a timeτ(0≤τ≤∞) the specific energy dissipation (energy per unit volume) within a closed domain Ω* (see Figure 1) is given by
(60)gf=∫0∞σij(x′,y′,τ)dϵij(x′,y′,τ)=∫0∞σijϵ̇ijdτ.For uniaxial deformation processgf would be, for a given point, the area under the stress-strain curve at that point. By taking the linearized geometric equations gf can be expressed as,(61)gf=∫0∞σij12(∂u̇i∂xj+∂u̇j∂xi)dτ=∫0∞σij∂u̇i∂xjdτ-∫0∞σij12(∂u̇i∂xj+∂u̇j∂xi)dτ.
The last integrand in (61) is zero, being the product of a symmetric and antisymmetric tensor, so that
(62)gf=∫0∞σij∂u̇i∂xjdτ=∫0∞∂∂xj(σiju̇i)dτ
where the Cauchy's equations for quasistatic processes and negligible body forces (∂σij/∂xi=0) have been considered. The total dissipated energy in the domain Ω* is
(63)W*=∫Ω*gfdΩ*=∫Ω*{∫0∞∂∂xj(σiju̇i)dτ}dΩ*.By applying the Gauss's theorem to (63) and using (58) we obtain,
(64)W*=∫Γ*{∫0∞σijnju̇idτ}dΓ*=∫Γ*{∫0∞tiu̇idτ}dΓ*Owing to the infitesimal width of the band, the curvilinear integral in (64) can be evaluated only on the lines Γ*+ and Γ*- (see Figure 1):
(65)W*=∫Γ*+∪Γ*-{∫0∞tiu̇idτ}dy′,
and taking into account (56) and (59) we obtain,
(66)W*=∫Γ*+∪Γ*-{∫0∞ti-(y′,τ)u̇i-(y′,τ)dτ}dy′+∫Γ*+∪Γ*-{∫0∞ti-(y′,τ)ϕ(x′,y′)ω̇(y′,τ)dτ}dy′.The first integral vanishes because the contributions onΓ*+ and Γ*- cancel each other out. Thus, the dissipated energy on Ω* is
(67)W*=∫Γ*+∪Γ*-{ϕ(x′,y′)∫0∞ti-(y′,τ)ω̇(y′,τ)dτ}dy′.Now, if the case whereΩ*=Ω is considered; that is, the whole band between points A and B in Figure 1, the total energy dissipated within the band between points A and B is
(68)W=∫ΩgfdΩ=∫yA′yB′{∫0∞ti-(y′,τ)ω̇(y′,τ)dτ}dy′.Equation (68) establishes that the energy dissipated within the idealized band can be written as a curvilinear integral along its length. The integrand of (68) represents the energy dissipated per unit of area, which in terms of fracture mechanics, is the fracture energy Gf:
(69)Gf(y′)=∫0∞ti-(y′,τ)ω̇(y′,τ)dτ.IfGf is assumed to be a material property independent of the spatial position of the point from (67) and (69) we obtain,
(70)W*=∫Γ*+∪Γ*-Gfϕ(x′,y′)dy′=Gf∫Γ*ϕ(x′,y′)dy′.By applying Green's theorem and taking into account (60), we can write the energy dissipated within the band as,
(71)W*=Gf∫Ω*∂ϕ∂x′dΩ*=∫Ω*gfdΩ*.The local form of (71) is
(72)gf=Gf∂ϕ∂x′=Gfl*,
where
(73)l*(x′,y′)=(∂ϕ∂x′)-1.The parameterl* plays the role of relating the specific energy (per unit of volume) and the fracture energy (per unit of area). l* is also identified as the characteristic length or crack band width used in existing cracking models [52]. For the unidimensional case and using the proposed Hermitian stress-strain softening law, the specific energy can be written as
(74)gf=σ0ϵf2,
where σ0 is the material strength and the degraded stress is given by
(75)σd=σ(1-d(ϵ)),
where the damage evolution law d(ϵ) is given in terms of the strain as follows:
(76)d(ϵ)=1-ϵ0ϵ[1+κ2(ϵ)(2κ(ϵ)-3)],
with
(77)κ(ϵ)=ϵ-ϵ0ϵf-ϵ0.Using  (74) the failure strain ϵf can be written in terms of the fracture energy and the characteristic length as it follows:
(78)ϵf=2Gfσ0l*.


### 2.6.4. Determination of the Functionϕ and the Characteristic Length l*

In order to apply the theory presented in the previous section to the discretized medium, consider a mesh ofC0 continuous hexahedron solid finite elements (see Figure 2).Figure 2
Determination of the characteristic length for hexahedron elements.A set of cracked elements is determined by using failure criteria for detecting the crack initiation, the crack orientation depends on the fibre direction. The cracked plane is defined here as a normal vector which is parallel to the fibres for fibre failure and normal to the fibre direction for matrix failure. The algorithm described in this section for determining the characteristic length was proposed by Oliver [57] and it has the advantages of calculating the characteristic length for arbitrary crack directions and any finite element geometry.From (71) the function ϕ has to be continuous and derivable, satisfying (57).  A simple function defined in the isoparametric coordinates ξ and η which fulfils these requirements is
(79)ϕ(ξ,η)=∑i=1ncNi(ξ,η)ϕi,
where nc is the number of corner nodes of a virtual plane located at the midplane of the element (nc=4 for our case), Ni are the standard C0 shape functions of an element of nc virtual nodes in its midplane and ϕi is the value of the ϕ at corner i. If the crack location inside the element is known, ϕi takes the value +1 if the corner node i is ahead the crack, and 0 otherwise. The function defined by (79) fulfils the required condition of continuity within elements and takes the values +1 for the nodes on the boundaries ahead the crack and 0 for the nodes on the boundaries behind the crack (see Figure 3).Figure 3
Finite element band modelling.In general, however, the exact crack location is not known, and usually only some indication of the onset of cracking and the crack directions is obtained at the integration points. The following algorithm proposed by Oliver [57] has been used for determining the characteristic length at each integration point j, as shown in Figure 4,Figure 4
Computation ofϕ values at the virtual midplane of the element.(1)  A set of local cartesian axesx′, y′ is defined at the centre of the element, this being identified by the values of the isoparametric coordinates (ξ=0, η=0 and ζ=0). The direction of the local axis x′ is defined by the normal to the fracture plane, which is the fibre angle for fibre failure (θj=θf) and (θj=θf+900) for matrix failure. 										(80){xi′yi′}=[cos(θj)sin(θj)-sin(θj)cos(θj)]{xiyi}(2)  Values ofϕ at each corner node are established according to their position with respect to the local axis x′, y′ (ϕi=1 if xi′≥0, otherwise ϕi=0).(3)  The characteristic length, at the present integration pointj with isoparametric coordinates ξj and ηj and cracking angle θj is obtained as
(81)l*(ξj,ηj)=(∂ϕ(ξj,ηj)∂x′)-1=(∑i=1nc[∂Ni(ξj,ηj)∂xcos(θj)+∂Ni(ξj,ηj)∂ysin(θj)]ϕi)-1
where
(82){∂N∂x∂N∂y}=Jxy-1{∂N∂ξ∂N∂η},
and Jxy is defined as the Jacobian matrix given by
(83)Jxy=[∂x∂ξ∂y∂ξ∂x∂η∂y∂η],
where the partial derivatives with respect to the isoparametric coordinates are written as(84)∂x∂ξ=14(1+η)x1-14(1+η)x2-14(1-η)x3+14(1-η)x4,∂x∂η=14(1+ξ)x1+14(1-ξ)x2-14(1-ξ)x3-14(1+ξ)x4,∂y∂ξ=14(1+η)y1-14(1+η)y2-14(1-η)y3+14(1-η)y4,∂y∂η=14(1+ξ)y1+14(1-ξ)y2-14(1-ξ)y3-14(1+ξ)y4,
where the pairs (xi,yi) refers to the global coordinates of the virtual nodes defining the midplane of the element.For transverse compression failure a set of cracked elements is determined by using the stress based criterion defined by (28) and the cracked plane is defined by the fracture angle θf. The function ϕ is given b(85)ϕ(ξ,ζ)=∑i=1ncNi(ξ,ζ)ϕi,
where nc is the number of corner nodes of a virtual plane located at the midplane of the element according to Figure 5, Ni are the linear shape functions defined previously and nc virtual nodes defining the virtual cracking midplane and ϕi is the value of the ϕ at corner i.Figure 5
Computation of the characteristic length for transverse compression.x′, z′ is an auxiliar coordinate system defined at the centre of the element, this being identified by the values of the isoparametric coordinates (ξ=0, η=0 and ζ=0) with the direction of the local axis x′ is defined by the normal to the fracture plane,
(86){xi′zi′}=[cos(θj)sin(θj)-sin(θj)cos(θj)]{xizi}.The values ofϕ at each corner node are established according to their position with respect to the local axis x′,z′ (ϕi=1 if xi′≥0, otherwise ϕi=0) in a similar way as the one described previously.The characteristic length associated with transverse compression at the present integration pointj with isoparametric coordinates ξj and ζj and fracture angle θf is given by
(87)l*(ξj,ζj)=(∂ϕ(ξj,ζj)∂x′)-1=(∑i=1nc[∂Ni(ξj,ζj)∂xcos(θf)+∂Ni(ξj,ζj)∂zsin(θf)]ϕi)-1
where
(88){∂N∂x∂N∂z}=Jxz-1{∂N∂ξ∂N∂ζ},
and Jxz is defined as the Jacobian matrix given by
(89)Jxz=[∂x∂ξ∂z∂ξ∂x∂ζ∂z∂ζ],
where the partial derivatives with respect to the isoparametric coordinates are written as(90)∂x∂ξ=14(1+ζ)x1-14(1+ζ)x2-14(1-ζ)x3+14(1-ζ)x4,∂x∂ζ=14(1+ξ)x1+14(1-ξ)x2-14(1-ξ)x3-14(1+ξ)x4,∂z∂ξ=14(1+ζ)z1-14(1+ζ)z2-14(1-ζ)z3+14(1-ζ)z4,∂z∂ζ=14(1+ξ)z1+14(1-ξ)z2-14(1-ξ)z3-14(1+ξ)z4,
where the pairs (xi,zi) refer to the global coordinates of the virtual nodes defining the midplane of the element. In-plane shear cracking is strongly dependent on the fibre orientarion within the element therefore, the characteristic length associated with in-plane shear failure has been assumed to be the same as the one defined for fibre failure or failure in the warp direction. For out-of-plane shear failure modes, cracks are assumed to be smeared over the thickness of the element with a crack band defined between upper and lower faces of the element which is equivalent to assume θf=900 in (87),
(91)l*(ξj,ζj)=(∂ϕ(ξj,ζj)∂x′)-1=(∑i=1nc[∂Ni(ξj,ζj)∂z]ϕi)-1


## 2.6.1. Crack Modelling in the Continuous Medium

Imagine a singular line in a two dimensional domain as a continuous material line, across which displacements are continuous but displacement gradients are discontinuous. The condition for a point belonging to a singular line with unit normaln at this point is that the determinant of the acoustic tensor in the n direction be zero, that is [57],
(50)det[niCijklnl]=0.
Nonpositive materials bifurcate, producing singular lines and the equation above permits their direction to be determined at each point. In the context of standard finite elements of C0 continuity, a singular line can be modelled only by the sides of the elements, these being the only points in the mesh where displacement gradient discontinuities can be obtained. However, a crack produces not only displacement gradient discontinuities but also displacement discontinuities. This latter kind of discontinuity cannot be modelled by a C0 finite element mesh for finite levels of discretization. However, a displacement discontinuity can be modeled as the limit of two parallel singular lines Γ- and Γ+ which tend to coincide with each other. The band delimited by these lines is known as singular band, and h is its width.By assuming an orthogonal curvilinear coordinate system (x′,y′) in the interior of the band, where y′ coordinates lines are parallel to the singular lines Γ- and Γ+, and x′ are the straight coordinates lines. Let u+(y′) and u-(y′) be the displacement vectors on Γ+ and Γ- the relative displacement vector can be written as,
(51)ω(y′)=u+(y′)-u-(y′)
as a vector representing the displacement “jump” between the two singular lines and
(52)δ(y′)=limh↦0ω(y′)=limh↦0(u+(y′)-u-(y′)),
If δ≠0, the singular band is modelling a discontinuous displacement field as the limit of a continuous one. This allows a crack to be idealised as a limit (with mesh refinement) of a band of finite elements where, by means of some numerical mechanisms, the condition δ≠0 is satisfied.


## 2.6.2. Displacement and Traction Vectors in the Singular Band

Consider a singular band in the solid, with a width h according to Figure1.Analysis within the singular band.
(a)(b)Along a coordinate linex′, the displacement vector u can be expanded from its value in the Γ- line using Taylor's series as
(53)u(x′,y′)=u-(y′)+[∂u∂x′]-Δx′+O(h2),
and consenquently, for a point in the Γ+ line
(54)u+(y′)=u-(y′)+[∂u∂x′]-h(y′)+O(h2).
From the (51), (53), (54) we can write
(55)u(x′,y′)=u-(y′)+[∂u∂x′]-Δx′+O(h2),(56)u(x′,y′)=u-(y′)+Δx′hω(y′)+0(h2)≅u-(y′)+ϕ(x′,y′)ω(y′),
where ϕ is a function to be determined, which approximates Δx′/h when h↦0.From (51) and (56) it can be seen that
(57)ϕ=0⇒Γ-,ϕ=1⇒Γ+.
The equilibrium across the singular band will be enforced by assuming the traction vector t acting on the plane defined by the normal n:
(58)ti=σijnj,
which is constant in the x′ direction that is
(59)t(x′,y′)=t+(y′)=t-(y′).


## 2.6.3. Energy Dissipation Within the Band

For a generic deformation process which takes place over a timeτ(0≤τ≤∞) the specific energy dissipation (energy per unit volume) within a closed domain Ω* (see Figure 1) is given by
(60)gf=∫0∞σij(x′,y′,τ)dϵij(x′,y′,τ)=∫0∞σijϵ̇ijdτ.For uniaxial deformation processgf would be, for a given point, the area under the stress-strain curve at that point. By taking the linearized geometric equations gf can be expressed as,(61)gf=∫0∞σij12(∂u̇i∂xj+∂u̇j∂xi)dτ=∫0∞σij∂u̇i∂xjdτ-∫0∞σij12(∂u̇i∂xj+∂u̇j∂xi)dτ.
The last integrand in (61) is zero, being the product of a symmetric and antisymmetric tensor, so that
(62)gf=∫0∞σij∂u̇i∂xjdτ=∫0∞∂∂xj(σiju̇i)dτ
where the Cauchy's equations for quasistatic processes and negligible body forces (∂σij/∂xi=0) have been considered. The total dissipated energy in the domain Ω* is
(63)W*=∫Ω*gfdΩ*=∫Ω*{∫0∞∂∂xj(σiju̇i)dτ}dΩ*.By applying the Gauss's theorem to (63) and using (58) we obtain,
(64)W*=∫Γ*{∫0∞σijnju̇idτ}dΓ*=∫Γ*{∫0∞tiu̇idτ}dΓ*Owing to the infitesimal width of the band, the curvilinear integral in (64) can be evaluated only on the lines Γ*+ and Γ*- (see Figure 1):
(65)W*=∫Γ*+∪Γ*-{∫0∞tiu̇idτ}dy′,
and taking into account (56) and (59) we obtain,
(66)W*=∫Γ*+∪Γ*-{∫0∞ti-(y′,τ)u̇i-(y′,τ)dτ}dy′+∫Γ*+∪Γ*-{∫0∞ti-(y′,τ)ϕ(x′,y′)ω̇(y′,τ)dτ}dy′.The first integral vanishes because the contributions onΓ*+ and Γ*- cancel each other out. Thus, the dissipated energy on Ω* is
(67)W*=∫Γ*+∪Γ*-{ϕ(x′,y′)∫0∞ti-(y′,τ)ω̇(y′,τ)dτ}dy′.Now, if the case whereΩ*=Ω is considered; that is, the whole band between points A and B in Figure 1, the total energy dissipated within the band between points A and B is
(68)W=∫ΩgfdΩ=∫yA′yB′{∫0∞ti-(y′,τ)ω̇(y′,τ)dτ}dy′.Equation (68) establishes that the energy dissipated within the idealized band can be written as a curvilinear integral along its length. The integrand of (68) represents the energy dissipated per unit of area, which in terms of fracture mechanics, is the fracture energy Gf:
(69)Gf(y′)=∫0∞ti-(y′,τ)ω̇(y′,τ)dτ.IfGf is assumed to be a material property independent of the spatial position of the point from (67) and (69) we obtain,
(70)W*=∫Γ*+∪Γ*-Gfϕ(x′,y′)dy′=Gf∫Γ*ϕ(x′,y′)dy′.By applying Green's theorem and taking into account (60), we can write the energy dissipated within the band as,
(71)W*=Gf∫Ω*∂ϕ∂x′dΩ*=∫Ω*gfdΩ*.The local form of (71) is
(72)gf=Gf∂ϕ∂x′=Gfl*,
where
(73)l*(x′,y′)=(∂ϕ∂x′)-1.The parameterl* plays the role of relating the specific energy (per unit of volume) and the fracture energy (per unit of area). l* is also identified as the characteristic length or crack band width used in existing cracking models [52]. For the unidimensional case and using the proposed Hermitian stress-strain softening law, the specific energy can be written as
(74)gf=σ0ϵf2,
where σ0 is the material strength and the degraded stress is given by
(75)σd=σ(1-d(ϵ)),
where the damage evolution law d(ϵ) is given in terms of the strain as follows:
(76)d(ϵ)=1-ϵ0ϵ[1+κ2(ϵ)(2κ(ϵ)-3)],
with
(77)κ(ϵ)=ϵ-ϵ0ϵf-ϵ0.Using  (74) the failure strain ϵf can be written in terms of the fracture energy and the characteristic length as it follows:
(78)ϵf=2Gfσ0l*.


## 2.6.4. Determination of the Functionϕ and the Characteristic Length l*

In order to apply the theory presented in the previous section to the discretized medium, consider a mesh ofC0 continuous hexahedron solid finite elements (see Figure 2).Figure 2
Determination of the characteristic length for hexahedron elements.A set of cracked elements is determined by using failure criteria for detecting the crack initiation, the crack orientation depends on the fibre direction. The cracked plane is defined here as a normal vector which is parallel to the fibres for fibre failure and normal to the fibre direction for matrix failure. The algorithm described in this section for determining the characteristic length was proposed by Oliver [57] and it has the advantages of calculating the characteristic length for arbitrary crack directions and any finite element geometry.From (71) the function ϕ has to be continuous and derivable, satisfying (57).  A simple function defined in the isoparametric coordinates ξ and η which fulfils these requirements is
(79)ϕ(ξ,η)=∑i=1ncNi(ξ,η)ϕi,
where nc is the number of corner nodes of a virtual plane located at the midplane of the element (nc=4 for our case), Ni are the standard C0 shape functions of an element of nc virtual nodes in its midplane and ϕi is the value of the ϕ at corner i. If the crack location inside the element is known, ϕi takes the value +1 if the corner node i is ahead the crack, and 0 otherwise. The function defined by (79) fulfils the required condition of continuity within elements and takes the values +1 for the nodes on the boundaries ahead the crack and 0 for the nodes on the boundaries behind the crack (see Figure 3).Figure 3
Finite element band modelling.In general, however, the exact crack location is not known, and usually only some indication of the onset of cracking and the crack directions is obtained at the integration points. The following algorithm proposed by Oliver [57] has been used for determining the characteristic length at each integration point j, as shown in Figure 4,Figure 4
Computation ofϕ values at the virtual midplane of the element.(1)  A set of local cartesian axesx′, y′ is defined at the centre of the element, this being identified by the values of the isoparametric coordinates (ξ=0, η=0 and ζ=0). The direction of the local axis x′ is defined by the normal to the fracture plane, which is the fibre angle for fibre failure (θj=θf) and (θj=θf+900) for matrix failure. 										(80){xi′yi′}=[cos(θj)sin(θj)-sin(θj)cos(θj)]{xiyi}(2)  Values ofϕ at each corner node are established according to their position with respect to the local axis x′, y′ (ϕi=1 if xi′≥0, otherwise ϕi=0).(3)  The characteristic length, at the present integration pointj with isoparametric coordinates ξj and ηj and cracking angle θj is obtained as
(81)l*(ξj,ηj)=(∂ϕ(ξj,ηj)∂x′)-1=(∑i=1nc[∂Ni(ξj,ηj)∂xcos(θj)+∂Ni(ξj,ηj)∂ysin(θj)]ϕi)-1
where
(82){∂N∂x∂N∂y}=Jxy-1{∂N∂ξ∂N∂η},
and Jxy is defined as the Jacobian matrix given by
(83)Jxy=[∂x∂ξ∂y∂ξ∂x∂η∂y∂η],
where the partial derivatives with respect to the isoparametric coordinates are written as(84)∂x∂ξ=14(1+η)x1-14(1+η)x2-14(1-η)x3+14(1-η)x4,∂x∂η=14(1+ξ)x1+14(1-ξ)x2-14(1-ξ)x3-14(1+ξ)x4,∂y∂ξ=14(1+η)y1-14(1+η)y2-14(1-η)y3+14(1-η)y4,∂y∂η=14(1+ξ)y1+14(1-ξ)y2-14(1-ξ)y3-14(1+ξ)y4,
where the pairs (xi,yi) refers to the global coordinates of the virtual nodes defining the midplane of the element.For transverse compression failure a set of cracked elements is determined by using the stress based criterion defined by (28) and the cracked plane is defined by the fracture angle θf. The function ϕ is given b(85)ϕ(ξ,ζ)=∑i=1ncNi(ξ,ζ)ϕi,
where nc is the number of corner nodes of a virtual plane located at the midplane of the element according to Figure 5, Ni are the linear shape functions defined previously and nc virtual nodes defining the virtual cracking midplane and ϕi is the value of the ϕ at corner i.Figure 5
Computation of the characteristic length for transverse compression.x′, z′ is an auxiliar coordinate system defined at the centre of the element, this being identified by the values of the isoparametric coordinates (ξ=0, η=0 and ζ=0) with the direction of the local axis x′ is defined by the normal to the fracture plane,
(86){xi′zi′}=[cos(θj)sin(θj)-sin(θj)cos(θj)]{xizi}.The values ofϕ at each corner node are established according to their position with respect to the local axis x′,z′ (ϕi=1 if xi′≥0, otherwise ϕi=0) in a similar way as the one described previously.The characteristic length associated with transverse compression at the present integration pointj with isoparametric coordinates ξj and ζj and fracture angle θf is given by
(87)l*(ξj,ζj)=(∂ϕ(ξj,ζj)∂x′)-1=(∑i=1nc[∂Ni(ξj,ζj)∂xcos(θf)+∂Ni(ξj,ζj)∂zsin(θf)]ϕi)-1
where
(88){∂N∂x∂N∂z}=Jxz-1{∂N∂ξ∂N∂ζ},
and Jxz is defined as the Jacobian matrix given by
(89)Jxz=[∂x∂ξ∂z∂ξ∂x∂ζ∂z∂ζ],
where the partial derivatives with respect to the isoparametric coordinates are written as(90)∂x∂ξ=14(1+ζ)x1-14(1+ζ)x2-14(1-ζ)x3+14(1-ζ)x4,∂x∂ζ=14(1+ξ)x1+14(1-ξ)x2-14(1-ξ)x3-14(1+ξ)x4,∂z∂ξ=14(1+ζ)z1-14(1+ζ)z2-14(1-ζ)z3+14(1-ζ)z4,∂z∂ζ=14(1+ξ)z1+14(1-ξ)z2-14(1-ξ)z3-14(1+ξ)z4,
where the pairs (xi,zi) refer to the global coordinates of the virtual nodes defining the midplane of the element. In-plane shear cracking is strongly dependent on the fibre orientarion within the element therefore, the characteristic length associated with in-plane shear failure has been assumed to be the same as the one defined for fibre failure or failure in the warp direction. For out-of-plane shear failure modes, cracks are assumed to be smeared over the thickness of the element with a crack band defined between upper and lower faces of the element which is equivalent to assume θf=900 in (87),
(91)l*(ξj,ζj)=(∂ϕ(ξj,ζj)∂x′)-1=(∑i=1nc[∂Ni(ξj,ζj)∂z]ϕi)-1


## 3. Numerical Implementation

This section presents details about the numerical implementation of the proposed failure model into ABAQUS FE code. The code formulation is based on the updated Lagrangian formulation which is used in conjunction with the central difference time integration scheme for integrating the resultant set of nonlinear dynamic equations. The method assumes a linear interpolation for velocities between two subsequent time steps and no stiffness matrix inversions are required during the analysis. The explicit method is conditionally stable for nonlinear dynamic problems and the stability for its explicit operator is based on a critical value of the smallest time increment for a dilatational wave to cross any element in the mesh. The numerical implementation steps are listed as follows:(1)  Stresses and strain-rates are transformed from the element axes to material axes:(92)(σld)n=T(σgd)n,ϵ̇ln+1=Tϵ̇gn+1,
where the indices l and g refer to the material (local) and element (global) coordinate systems, respectively and T is the transformation matrix defined by (2). The superscripts n and n+1 refer to the previous and current time, respectively:(2)  Compute constitutive matrixC from (3) for trial stresses.(3)  Based on the current time step compute strain increments and update the elastic stresses and strains(93)σn+1=σn+CΔϵ,ϵn+1=ϵn+Δϵ.(4)  Compute failure index for fibre failure (σ1n+1>0):(94)F1t(σ1n+1)=σ1n+1Xt≥1.
Store (ϵ1,0t)n+1 when F1t(σ1n+1)=1.(5)  Compute failure index for fibre kinking (σ1n+1<0):
(95)F1c(σ1n+1)=|σ1n+1|Xc≥1.
Store (ϵ1,0c)n+1 when F1c(σ1n+1)=1.(6)  Compute IFF index for tensile/shear matrix cracking (σ2n+1>0):
(96)F2t(σ2n+1,τ23n+1,τ12n+1)=(σ2n+1Yt)2+(τ23n+1S23)2+(τ12n+1S12)2≥1.
Store (ϵ̅m,0t)n+1 and (σ̅m,0t)n+1 when F2t(σ2n+1,τ23n+1,τ12n+1)=1.(7)  Compute IFF index for matrix crushing and intralaminar shear cracking(σ2n+1<0):(97)F2c(σntn+1,σnln+1)=(σntn+1S23A+μntσnnn+1)2+(σnln+1S12+μnlσnnn+1)2≥1.
Store (ϵ̅m,0c)n+1,(σ̅m,0c)n+1 and (θf)n+1 when F2c(σntn+1,σnln+1)=1.(8)  Compute failure index for inplane shear cracking:(98)F12(τ12n+1)=|τ12n+1|S12≥1.
Store (γ12,0)n+1,(γ12,0i)n+1 when F12(τ12n+1)=1.(9)  Compute characteristic lengths for each failure mode using objectivity algorithm:In-plane failure modes(99)l*(ξj,ηj)=(∂ϕ(ξj,ηj)∂x′)-1=(∑i=1nc[∂Ni(ξj,ηj)∂xcos(θj)+∂Ni(ξj,ηj)∂ysin(θj)]ϕi)-1θj=θfibre for fibre failure in tension and compression, θj=θfibre+900 for matrix cracking in tension/shear, θj=θfibre for in-plane shear failure.Out-of-plane failure modes(100)l*(ξj,ζj)=(∂ϕ(ξj,ζj)∂x′)-1=(∑i=1nc[∂Ni(ξj,ζj)∂z]ϕi)-1,θj=θf,where θf is the fracture plane orientation when F2c(σnt,σnl)=1 for matrix crushing/intralaminar shear failure modes.(10)  Compute final strains based on intralaminar fracture toughness and characteristic length associated with each failure mode:Fibre failure(101)ϵ1,ft=2GfibretXtl*,ϵ1,fc=2GfibrecXcl*,IFF (Matrix cracking in tension/shear)(102)ϵ̅m,ft=2σ̅m,0t[(cos2(θ)gmct)λ+(sin2(θ)gmcs)λ]-1/λ,IFF (Matrix cracking in compression/shear)(103)ϵ̅m,fc=2Gmatrixcσ̅m,0cl*=2GIIcσ̅m,0cl*,in-plane shear failure(104)γ12,f=2GshearS12l*.(11)  Update damage parametersd1, dmt, dmc, d12 when the failure criteria are met with,(105)d1={(d1)n+1if(d1)n+1>(d1)n,F1t(σ1n+1)>1orF1c(σ1n+1)>1,(d1)notherwise,dmt={(dmt)n+1if(dmt)n+1>(dmt)n,F2t(σ2n+1,τ23n+1,τ12n+1)>1,(dmt)notherwise,dmc={(dmc)n+1if(dmc)n+1>(dmc)n,F2c(σntn+1,σnln+1)>1,(dmc)notherwise,d12={(d12)n+1if(d12)n+1>(d12)n,F12(τ12n+1)>1,(d12)notherwise.(12)  Stress update(106){σ1dσ2dσ3dτ12dτ23dτ13d}n+1={(1-d1)σ1(1-dmt)σ2dmcσ3dmc(1-dmt)(1-d12)τ12dmc(1-dmt)τ23dmcτ13dmc}.(13)  Stress transformation from the material(local) to the element coordinate system (global)(107)(σgd)n+1=(T)-1(σld)n+1.(14)  Compute element nodal forces, add element hourglass control nodal forces, transform element nodal forces to global coordinates, solve for accelerations via equation of motion and update velocities and displacements.


## 4. Single Element Validation

Numerical simulations at element level were carried out to validate the proposed damage model. The numerical tests consisted of exciting each failure mode individually and verify if damage irreversibility conditions and energy concepts are always satisfied for each failure mode. The material properties used in model were taken from  [58] and they are summarised in Tables 1, 2 and 3. The material system consists of a quasi unidirectional carbon UD tape supplied by EUROCARBON, which has a plain weave pattern with T700 carbon fibres in the warp direction and a small fraction of PPG glass fibres in weft direction to hold the carbon fibres together embedded into an infusible PRIME 20 LV epoxy resin system.Table 1
Mechanical properties of each ply.
E11 (GPa)E22=E33 (GPa)G12 (GPa)G13 (GPa)G23 (GPa)ν12=ν13ν231008.114.654.655.00.30.4Table 2
Ply strengths.
Xt (MPa)Xc (MPa)Yt (MPa)Yc (MPa)S12 (MPa)20001000100160140Table 3
Intralaminar fracture toughnesses.
Gft (kJ/m2)Gfc (kJ/m2)Gmt (kJ/m2)Gmc (kJ/m2)Gs (kJ/m2)100.025.02.02.02.0The element was loaded under displacement control in each direction with prescribed displacementsU1=Up and boundary conditions shown in Figures 6, 7 and 8. The stress in the fibre direction was normalised with respect to the fibre tensile strength whilst the stress in the matrix (transverse) direction was normalised with respect to the matrix compression strength. The shear stress was normalised with respect to the shear strength. The structural responses for one element loaded in tension, compression and combined tension-compression under loading-unloading-reloading conditions in the fibre and matrix directions are shown in Figures 9 and 11, respectively. As damage commences the stresses are gradually reduced to zero. The material stiffness in each direction are also reduced as damage cummulates and during unloading the damage irreversibility condition is fully satisfied in order to avoid material self-healing as shown in Figures 10, 12 and 14. The nonlinear behaviour in shear including damage combined with plasticity effects is shown in Figure 13.Figure 6
Boundary conditions for tensile loading.Figure 7
Boundary conditions for compression loading.Figure 8
Boundary conditions for in-plane shear loading.Figure 9
Element loaded-unloaded-reloaded in tension and compression in the fibre direction.Figure 10
Damage evolution law due to combined tension and compression loading in the fibre direction.Figure 11
Element loaded-unloaded-reloaded in tension and compression in the matrix direction.Figure 12
Damage evolution law due to combined tension and compression loading in the matrix direction.Figure 13
Element loaded-unloaded-reloaded in shear.Figure 14
Damage evolution law in shear.


## 5. Coupon Tests Validation

In this section the predictions obtained using the proposed failure model are compared against experimental results at coupon level in tension compression and shear. The experimental results were taken from  [58] in which the first author tested composite laminates manufactured using the Resin Infusion under Flexible Tooling (RIFT) process. The mechanical properties of the material are listed in Tables 1, 2 and 3. The dimensions of the coupons used by the author for the tensile tests in the fibre and matrix directions were 200×20×1.35 mm3 and 200×20×2.25 mm3, respectively, with a gauge length of 100 mm for both cases. For the in-plane shear tests the specimen had dimensions of 200×20×2.7 mm3 with a gauge length of 100 mm and lay-up of [+450/-450]3. The material constants c1 and c2 required by the model to predict the inplane shear behaviour were obtained by finding the best-fit between the model and the experimental shear stress and shear strain curve. A comparison between the experimental and numerical shear stress-shear strain curves together with the material model parameters c1 and c2 are presented in Figure 15. As the experiments were carried out quasistatically, a high value in the order of 109 was assigned to c3 in order to neglect the strain rate effects in the numerical response in shear.Figure 15
Comparison between numerical and experimental curves for the in-plane shear behaviour.The lay-ups used for the specimens loaded in tension in both fibre and matrix directions were[00]3 and [00]5, respectively. The mechanical properties of each layer are listed in Tables 1, 2, 3. Each ply with a nominal thickness of 0.45 mm was individually modelled using a single solid hexahedral linear element through the thickness direction. The tabs behaviour was assumed to be orthotropic linear elastic. The virtual specimens were assumed to be fixed in one end and loaded quasistatically under displacement control using the dynamic relaxation method and the meshes used in the tensile and compression test simulations are shown in Figures 16 and 17.Figure 16
Finite element mesh used in the tensile test simulation.Figure 17
Finite element mesh used in the compression test simulation.Figures18–22 show the correlation between numerical predictions and experimental results for the tension and compression in both fibre and matrix directions and in-plane shear tests, respectively.Figure 18
Comparison between numerical predictions and experimental results for the fibre tensile tests.Figure 19
Comparison between numerical predictions and experimental results for the matrix tensile tests.Figure 20
Comparison between numerical predictions and experimental results for the fibre compression tests.Figure 21
Comparison between numerical predictions and experimental results for the matrix compression tests.Figure 22
Comparison between numerical predictions and experimental results for the in-plane shear tests.In general, the numerical predictions agree very well with the experimental results for all simulations carried out at coupon level. It can be clearly seen from Figure18 the stiffening effect in tension in the longitudinal direction. This effect has not been taken into account in the model development and it is due to the fabric arquitecture which leads to the development of high interlaminar shear stresses between fibres located in the fill and warp direction during the loading process which results in a reduction of the initial crimp angle of the fabric therefore increasing the laminate stiffness along the loading history.Comparisons between the experimental and predicted failure locations are shown in Figures23–26. The numerical predictions indicate that the use of sharp corners end tabs leads to failure near to the tab due to the high stress concentrations in these regions. This effect was also observed in the experiments. The prediction of the failure location for transverse compression, Figure 25, correlate very well with the experiments indicating that the failure in transverse direction is a shear dominated failure mode with a fracture plane orientated about 500 in respect to the through the thickness direction. In order to avoid element distortion problems, the failed elements in transverse compression have been deleted from the mesh using an erosion algorithm available in the ABAQUS Finite Element Code.Comparison between predicted and experimental failure locations for the tensile tests.
(a)
Numerical prediction(b)
ExperimentalComparison between predicted and experimental failure locations for the compression test in the fibre direction.
(a)
Numerical prediction(b)
ExperimentalComparison between predicted and experimental failure locations for the compression test in the transverse direction.
(a)
Numerical prediction(b)
ExperimentalComparison between predicted and experimental failure locations for the in-plane shear test.
(a)
Numerical prediction(b)
Experimental


## 6. Conclusions

A detailed formulation for a three dimensional failure model for composite laminates under multiaxial loading has been presented and discussed in this paper. The relevant contributions of the present formulation are the following:(i)
incorporation of strain rate effects on in-plane shear behaviour by means of a semiempirical viscoplastic constitutive law,(ii)
implementation of thermodynamically consistent degradation laws coupled with an advanced objectivity algorithm which avoids strain localization and mesh dependence problems regardless of mesh refinement and element topology,(iii)
high order damage evolution laws have been incorporated into the model in order to avoid generation of artifitial stress waves during damage growth,(iv)
the model formulation combines a quadratic stress based criterion with the matrix compression criterion based on the action plane to predict interfibre-failure without knowing a priori the orientation of the fracture plane(v)
usage of the Second Piola Kirchhoff stresses which potentially avoids material deformation problems, such as trellising,(vi)
implementation of the new failure model into ABAQUS/Explicit FE code within brick elements as an user-defined material model.Single element simulations were carried out to validade the model at element level. Comparisons between numerical predictions and experimental results at coupon level under different loading conditions were carried out to validade the proposed failure model. A very good correlation between numerical predictions and experimental results was obtained using the proposed failure model.



---

*Source: 486063-2009-09-28.xml*