**Finite Element Implementation of Failure and Damage Simulation in Composite Plates**

Milan Žmindák and Martin Dudinský

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/46498

## **1. Introduction**

130 Composites and Their Properties

Tsai, S.W., & Wu, E.M. (1971) A General Theory of Strength for Anisotropic Materials.

Wilkins, M.L., & Guinan, M.W. (1973) Impact of cylinders on a rigid boundary. *Journal of* 

*Journal of Composite Materials,* Vol. 5, No. 1, pp. 58–80

*Applied Physics,* Vol. 44, No. 3, pp. 1200–1206

Composite materials are now common engineering materials used in a wide range of applications. They play an important role in the aviation, aerospace and automotive industry, and are also used in the construction of ships, submarines, nuclear and chemical facilities, etc.

The meaning of the word damage is quite broad in everyday life. In continuum mechanics the term damage is referred to as the reduction of the internal integrity of the material due to the generation, spreading and merging of small cracks, cavities and similar defects. Damage is called elastic, if the material deforms only elastically (in macroscopic level) before the occurrence of damage, as well as during its evolution. This damage model can be used if the ability of the material to deform plastically is low. Fiber-reinforced polymer matrix composites can be considered as such materials.

The use of composite materials in the design of constructions is increasing in traditional structures such as for example development of airplanes or in the automotive industry. Recently this kind of materials is used in development of special technique and rotating systems such as propellers, compressor turbine blades etc. Other applications are in electronics, electrochemical industry, environmental and biomedical engineering (Chung, 2003).

The costs for designing of composite structures is possible partially eliminate by numerical simulation of solving problem. In this case the simulation is not accepted as universal tool for analyzing of systems behaviour but it is an effective alternative to processes of experimental sciences. Simulations support development of new theories and suggestion of new experiments for testing these theories. Experiments are necessary for obtaining of input data into simulation programs and for verification of numerical programs and models.

<sup>© 2012</sup> Žmindák and Dudinský, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits © 2012 Žmindák and Dudinský, licensee InTech. This is a paper distributed under the terms of the Creative Commons

unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Laminated composites have a lot of advantages but in some cases they show different limitations that are caused by stress concentrations between layers. Discontinuous change of material properties is reason for occurrence of interlaminar stresses that often cause delamination failure (Zhang & Wang, 2009). Delamination may originate from manufacturing imperfections, cracks produced by fatigue or low velocity impact, stress concentration near geometrical/material discontinuity such as joints and free edges, or due to high interlaminar stresses (Elmarakbi, 2009)

Finite Element Implementation of Failure and Damage Simulation in Composite Plates 133

other hand the micromechanics-based approaches conduct micromechanical analysis of representative volume element (RVE) with subsequent homogenization to predict evolving material damage behavior (Mishnaevsky, 2007). Most damage models do not account for the evolution of damage or the effect of loading history (Jain & Ghosh, 2009). Significant error can consequently accrue in the solution of problems, especially those that involve nonproportional loading. Some of these homogenization studies have overcome this shortcoming through the introduction of simultaneous RVE-based microscopic and macroscopic analysis in each load step. However, such approaches can be computationally very expensive since detailed micro-mechanical analyses need to be conducted in each load

Jain & Ghosh (2009) have developed a 3D homogenization-based continuum damage mechanics (HCDM) model for fiber-reinforced composites undergoing micro-mechanical damage. Micromechanical damage in the RVE is explicitly incorporated in the form of fibermatrix interfacial debonding. The model uses the evolving principal damage coordinate system as its reference in order to represent the anisotropic coefficients, which is necessary for retaining accuracy with nonproportional loading. The HCDM model parameters are calibrated by using homogenized micromechanical solutions for the RVE for a few strain

There are many works written about damage of composite plates and many models for various types of damage, plates or loading have been developed. Shu (2006) presented generalized model for laminated composite plates with interfacial damage. This model deals with three kinds of interfacial debonding conditions: perfect bonding, weak bonding and delamination. Iannucci & Ankersen (2006) described unconventional energy based composite damage model for woven and unidirectional composite materials. This damage model has been implemented into FE codes for shell elements, with regard to tensile, compressive and shear damage failure modes. Riccio & Pietropaoli (2008) dealt with modeling damage propagation in composite plates with embedded delamination under compressive load. The influence of different failure mechanisms on the compressive behavior of delaminated composite plates was assessed, by comparing numerical results obtained with models characterized by different degrees of complexity. Tiberkak et al. (2008) studied damage prediction in composite plates subjected to low velocity impact. Fiber-reinforced composite plates subjected to low velocity impact were studied by use of finite element analysis where Mindlin's plate theory and 9-node Lagrangian element were considered. Clegg et al. (2006) worked out interesting study of hypervelocity impact damage prediction in composites. This study reports on the development of an extended orthotropic continuum material model and associated material characterization techniques for the simulation and validation of impacts onto fiber-reinforced composite materials. The model allows to predict the extent of damage and residual strength of the fiber-reinforced

Many studies of an effect of various aspects of damage process on behavior of composite plates can be found in the literature (e.g. Gayathri et al. , 2010). Many authors have been

step at every integration point in elements of the macroscopic structure.

histories.

composite material after impact.

Delaminations in layered plates and beams have been analyzed by using both cohesive damage models and fracture mechanics. Cohesive elements are widely used, in both forms of continuous interface elements and point cohesive elements (Cui, 1993), at the interface between solid finite elements to predict and to understand the damage behaviour in the interfaces of different layers in composite laminates. In the context of the fracture mechanics approach (Sládek, et al. , 2002), it allows us to predict the growth of a preexisting crack or defect. In a homogeneous and isotropic body subjected to a general loading condition, a crack tends to grow by kinking in a direction such a pure mode I condition at its tip is maintained. On the contrary, delaminations in laminated composites are constrained to propagate in its own plane because the toughness of the interface is relatively low in comparison to that of the adjoining material. Therefore a delamination crack propagates with its advancing tip in mixed mode condition and, consequently, requires a fracture criterion including all three mode components.

The theory of crack growth may be developed by using one of two approaches. First, the Griffith energetic (or global) approach introduces the concept of energy release rate (ERR) *G* as the energy available for fracture on one hand, and the critical surface energy *Gr* as the energy necessary for fracture on the other hand. Alternatively, the Irwin (local) approach is based on the stress intensity factor concept, which represents the energy stress field in the neighborhood of the crack tip. These two approaches are equivalent and, therefore, the energy criterion may be rewritten in terms of stress intensity factors.

Microcracking in a material is almost always associated with changes in mechanical behavior of the material. The problem of microcracking in fiber-reinforced composites is complicated due to the multitude of different microcracking modes which may initiate and evolve independently or simultaneously. Continuum Damage Mechanics (CDM) considers damaged materials as a continuum, in spite of heterogenity, micro-cavities, and microdefects and is based on expressing of stiffness reduction caused by damage, by establishing effective damage parameters which represent cumulative degradation of material. There are basically two categories of CDM models used for estimating the constitutive behavior of composite materials containing microcracks –phenomenological models and micromechanics models.

The phenomenological CDM models employ scalar, second order or fourth order tensors using mathematically and thermodynamically consistent formulations of damage mechanics. Damage parameters are identified through macroscopic experiments and in general, they do not explicitly account for damage mechanism in the microstructure. On the other hand the micromechanics-based approaches conduct micromechanical analysis of representative volume element (RVE) with subsequent homogenization to predict evolving material damage behavior (Mishnaevsky, 2007). Most damage models do not account for the evolution of damage or the effect of loading history (Jain & Ghosh, 2009). Significant error can consequently accrue in the solution of problems, especially those that involve nonproportional loading. Some of these homogenization studies have overcome this shortcoming through the introduction of simultaneous RVE-based microscopic and macroscopic analysis in each load step. However, such approaches can be computationally very expensive since detailed micro-mechanical analyses need to be conducted in each load step at every integration point in elements of the macroscopic structure.

132 Composites and Their Properties

to high interlaminar stresses (Elmarakbi, 2009)

criterion including all three mode components.

micromechanics models.

energy criterion may be rewritten in terms of stress intensity factors.

Laminated composites have a lot of advantages but in some cases they show different limitations that are caused by stress concentrations between layers. Discontinuous change of material properties is reason for occurrence of interlaminar stresses that often cause delamination failure (Zhang & Wang, 2009). Delamination may originate from manufacturing imperfections, cracks produced by fatigue or low velocity impact, stress concentration near geometrical/material discontinuity such as joints and free edges, or due

Delaminations in layered plates and beams have been analyzed by using both cohesive damage models and fracture mechanics. Cohesive elements are widely used, in both forms of continuous interface elements and point cohesive elements (Cui, 1993), at the interface between solid finite elements to predict and to understand the damage behaviour in the interfaces of different layers in composite laminates. In the context of the fracture mechanics approach (Sládek, et al. , 2002), it allows us to predict the growth of a preexisting crack or defect. In a homogeneous and isotropic body subjected to a general loading condition, a crack tends to grow by kinking in a direction such a pure mode I condition at its tip is maintained. On the contrary, delaminations in laminated composites are constrained to propagate in its own plane because the toughness of the interface is relatively low in comparison to that of the adjoining material. Therefore a delamination crack propagates with its advancing tip in mixed mode condition and, consequently, requires a fracture

The theory of crack growth may be developed by using one of two approaches. First, the Griffith energetic (or global) approach introduces the concept of energy release rate (ERR) *G* as the energy available for fracture on one hand, and the critical surface energy *Gr* as the energy necessary for fracture on the other hand. Alternatively, the Irwin (local) approach is based on the stress intensity factor concept, which represents the energy stress field in the neighborhood of the crack tip. These two approaches are equivalent and, therefore, the

Microcracking in a material is almost always associated with changes in mechanical behavior of the material. The problem of microcracking in fiber-reinforced composites is complicated due to the multitude of different microcracking modes which may initiate and evolve independently or simultaneously. Continuum Damage Mechanics (CDM) considers damaged materials as a continuum, in spite of heterogenity, micro-cavities, and microdefects and is based on expressing of stiffness reduction caused by damage, by establishing effective damage parameters which represent cumulative degradation of material. There are basically two categories of CDM models used for estimating the constitutive behavior of composite materials containing microcracks –phenomenological models and

The phenomenological CDM models employ scalar, second order or fourth order tensors using mathematically and thermodynamically consistent formulations of damage mechanics. Damage parameters are identified through macroscopic experiments and in general, they do not explicitly account for damage mechanism in the microstructure. On the Jain & Ghosh (2009) have developed a 3D homogenization-based continuum damage mechanics (HCDM) model for fiber-reinforced composites undergoing micro-mechanical damage. Micromechanical damage in the RVE is explicitly incorporated in the form of fibermatrix interfacial debonding. The model uses the evolving principal damage coordinate system as its reference in order to represent the anisotropic coefficients, which is necessary for retaining accuracy with nonproportional loading. The HCDM model parameters are calibrated by using homogenized micromechanical solutions for the RVE for a few strain histories.

There are many works written about damage of composite plates and many models for various types of damage, plates or loading have been developed. Shu (2006) presented generalized model for laminated composite plates with interfacial damage. This model deals with three kinds of interfacial debonding conditions: perfect bonding, weak bonding and delamination. Iannucci & Ankersen (2006) described unconventional energy based composite damage model for woven and unidirectional composite materials. This damage model has been implemented into FE codes for shell elements, with regard to tensile, compressive and shear damage failure modes. Riccio & Pietropaoli (2008) dealt with modeling damage propagation in composite plates with embedded delamination under compressive load. The influence of different failure mechanisms on the compressive behavior of delaminated composite plates was assessed, by comparing numerical results obtained with models characterized by different degrees of complexity. Tiberkak et al. (2008) studied damage prediction in composite plates subjected to low velocity impact. Fiber-reinforced composite plates subjected to low velocity impact were studied by use of finite element analysis where Mindlin's plate theory and 9-node Lagrangian element were considered. Clegg et al. (2006) worked out interesting study of hypervelocity impact damage prediction in composites. This study reports on the development of an extended orthotropic continuum material model and associated material characterization techniques for the simulation and validation of impacts onto fiber-reinforced composite materials. The model allows to predict the extent of damage and residual strength of the fiber-reinforced composite material after impact.

Many studies of an effect of various aspects of damage process on behavior of composite plates can be found in the literature (e.g. Gayathri et al. , 2010). Many authors have been

dealing with problematic of damage of composite plated under cyclic loading or problematic of impact fatigue damage (e.g. Azouaoui et al. , 2010). Composite materials are becoming more and more used for important structural elements and structures, so the problematic of fatigue damage of composites is becoming more and more actual. Numerical implementation of damage is not simple. Finite element method (FEM) is the most utilized method for modeling damage. Fast Multipole BEM method or various meshless methods are also establishing at the present time.

Finite Element Implementation of Failure and Damage Simulation in Composite Plates 135

the displacements in the *z-*th plate element, in terms of a global reference system located at

the laminate mid-surface, are expressed by (e.g. Carrera 2002; Reddy 1995)

**Figure 1.** Delaminated composite plate

**Figure 2.** Laminate subdivision in plate elements

(

i i

w (x,y,z) w (x,y)

0 i i i xi <sup>0</sup> ( <sup>i</sup> <sup>i</sup> i yi 0

u (x,y,z) u (x,y) (z z ) (x,y) v x,y,z) v x,y) (z z ) (x,y)

 

(1)

Firstly, the goal of this chapter is to present the numerical results of the delamination analysis of two laminae with different thickness with two orthotropic material properties and subjected to a pair of opposed forces. For this goal we used commercial FEM software ANSYS and the mode I, II, and III components of energy release rate (ERR) were calculated. Secondly, the goal is to present the numerical results of elastic damage of thin composite plates. The analysis was performed by user own software, created in MATLAB programming language. This software can perform numerical analysis of elastic damage using FEM layered plate finite elements based on the Kirchhoff plate theory.

This chapter is organized in three sections: Second section is focused on failure modeling in laminates by using standard shear deformable elements, whereas interface elements were used for the interface model. The delamination propagation is controlled by the critical ERR. In third section, in first, a general description of damage is provided, then damage model used is examined. Finally, some numerical results obtained for damage of plate are presented.

## **2. Theoretical background of failure modeling in laminates**

The mechanisms that lead to failure in composite materials are not yet fully understood, especially for matrix or fiber compression. Strength-based failure criteria are commonly used with the FEM to predict failure events in composite structures. Numerous continuumbased criteria have been derived to relate internal stresses and experimental measures of material strength to the onset of failure (Dávila et al., 2005). In Fig. 1 a laminate contains a single in-plane delamination crack of area Ω<sup>D</sup> with a smooth front Ω<sup>D</sup> . The laminate thickness is denoted by *h0*. The *x-y* plane is taken to be the mid-plane of the laminate, and the z-axis is taken positive downwards from the mid-plane.

### **2.1. Plate finite elements for sublaminate modeling**

Each sublaminate is represented by an assembly of first order shear deformable (FSDT) plate elements bonded by zero-thickness interfaces in the transverse direction as shown in Fig. 2. The delamination plane separates the delaminated structure into two sublaminates of thickness *h1*, *h2* and each sublaminate consist the upper *nu* plates and the lower *nl* plates. Each plate element is composed from one or few physical fiber-reinforced plies with their material axes arbitrarily oriented. Lagrangian multipliers through constraint equations (CE) are used for enforcing adhesion between the plates inside each sublaminate. Accordingly, the displacements in the *z-*th plate element, in terms of a global reference system located at the laminate mid-surface, are expressed by (e.g. Carrera 2002; Reddy 1995)

**Figure 1.** Delaminated composite plate

134 Composites and Their Properties

presented.

are also establishing at the present time.

dealing with problematic of damage of composite plated under cyclic loading or problematic of impact fatigue damage (e.g. Azouaoui et al. , 2010). Composite materials are becoming more and more used for important structural elements and structures, so the problematic of fatigue damage of composites is becoming more and more actual. Numerical implementation of damage is not simple. Finite element method (FEM) is the most utilized method for modeling damage. Fast Multipole BEM method or various meshless methods

Firstly, the goal of this chapter is to present the numerical results of the delamination analysis of two laminae with different thickness with two orthotropic material properties and subjected to a pair of opposed forces. For this goal we used commercial FEM software ANSYS and the mode I, II, and III components of energy release rate (ERR) were calculated. Secondly, the goal is to present the numerical results of elastic damage of thin composite plates. The analysis was performed by user own software, created in MATLAB programming language. This software can perform numerical analysis of elastic damage

This chapter is organized in three sections: Second section is focused on failure modeling in laminates by using standard shear deformable elements, whereas interface elements were used for the interface model. The delamination propagation is controlled by the critical ERR. In third section, in first, a general description of damage is provided, then damage model used is examined. Finally, some numerical results obtained for damage of plate are

The mechanisms that lead to failure in composite materials are not yet fully understood, especially for matrix or fiber compression. Strength-based failure criteria are commonly used with the FEM to predict failure events in composite structures. Numerous continuumbased criteria have been derived to relate internal stresses and experimental measures of material strength to the onset of failure (Dávila et al., 2005). In Fig. 1 a laminate contains a single in-plane delamination crack of area Ω<sup>D</sup> with a smooth front Ω<sup>D</sup> . The laminate thickness is denoted by *h0*. The *x-y* plane is taken to be the mid-plane of the laminate, and

Each sublaminate is represented by an assembly of first order shear deformable (FSDT) plate elements bonded by zero-thickness interfaces in the transverse direction as shown in Fig. 2. The delamination plane separates the delaminated structure into two sublaminates of thickness *h1*, *h2* and each sublaminate consist the upper *nu* plates and the lower *nl* plates. Each plate element is composed from one or few physical fiber-reinforced plies with their material axes arbitrarily oriented. Lagrangian multipliers through constraint equations (CE) are used for enforcing adhesion between the plates inside each sublaminate. Accordingly,

using FEM layered plate finite elements based on the Kirchhoff plate theory.

**2. Theoretical background of failure modeling in laminates** 

the z-axis is taken positive downwards from the mid-plane.

**2.1. Plate finite elements for sublaminate modeling** 

$$\begin{aligned} \mathbf{u}\_{i}(\mathbf{x}, \mathbf{y}, \mathbf{z}) &= \mathbf{u}\_{i}^{0}(\mathbf{x}, \mathbf{y}) + (\mathbf{z} - \mathbf{z}\_{i}) \, \boldsymbol{\nu}\_{\text{xi}}(\mathbf{x}, \mathbf{y}) \\ \mathbf{v}\_{i}(\mathbf{x}, \mathbf{y}, \mathbf{z}) &= \mathbf{v}\_{i}^{0}(\mathbf{x}, \mathbf{y}) + (\mathbf{z} - \mathbf{z}\_{i}) \, \boldsymbol{\nu}\_{\text{yi}}(\mathbf{x}, \mathbf{y}) \\ \mathbf{w}\_{i}(\mathbf{x}, \mathbf{y}, \mathbf{z}) &= \mathbf{w}\_{i}^{0}(\mathbf{x}, \mathbf{y}) \end{aligned} \tag{1}$$

**Figure 2.** Laminate subdivision in plate elements

where *ui*, *vi* refer to the in-plane displacements, and *wi* to the transverse displacements through the thickness of the *i-*th plate element, 00 0 , , *ii i uvw* , are the displacements at the midsurface of the *i-*th plate element and xi yi (x,y), x,y denote rotations of transverse normals about *y* and *x*, respectively.

Finite Element Implementation of Failure and Damage Simulation in Composite Plates 137

<sup>2</sup> (4)

σ 1 D **K**Δ (5)

1 (6)

In order to avoid interpenetration between delaminated sublaminates in the delaminated region Ω<sup>D</sup> , a unilateral frictionless contact interface can be introduced, characterized by a zero stiffness for opening relative displacements (*Δw≥0*) and a positive stiffness for closing

zz

σ 1 sign(Δw)) kzΔw

where *kz* is the penalty number imposing contact constraint and sign is the signum function. A very large value for *kz* restricts sublaminate overlapping and simulates the contact condition. Unilateral contact conditions may be implemented in ANSYS using COMBIN39. This element is a unidirectional element with nonlinear constitutive relationships with

If we introduce a scalar damage variable *D* with the value of *1* for no adhesion and the value of 0 for perfect adhesion, we get a single extended interface model with constitutive law valid both for undelaminated Ω Ω<sup>D</sup> and delaminated ΩD areas. Consequently

In this work we use the formulation via FEs, related to plate elements, interface elements and Lagrange multipliers. It is worth noting that in commercial FEA packages the Lagrange multipliers are represented by either CE or rigid links, whereas interface elements are implemented by the analyst using a combination of spring elements

In order to predict crack propagation in laminates for general loading conditions, ERR distributions along the delamination front are needed. Fracture mechanics assumes that delamination propagation is controlled by the critical ERR. Delamination grows on the region of the delamination front where the following condition is satisfied and is of the

> <sup>α</sup> <sup>β</sup> I II III cc c I II III

where α, β and γ are mixed mode fracture parameters determined by fitting experimental

The critical ERR cc c , G ,G G material properties can be I II III evaluated from experimental

GG G 

procedures. The closed-form expressions for the ERR are (Barbero, 2008)

<sup>γ</sup> Gs G s G s

1

appropriate specialization of the nonlinear constitutive law according to (4).

*2.1.2. Contact formulation for damage interface* 

constitutive law can be expressed as

(COMBIN14) and CE.

form

test results.

*2.1.3. Mixed mode analysis* 

relative displacements (*Δw≤0*), then the contact stress zz σ is

At the reference surfaces, the membrane strain vector i ε , the curvature i κ , and the transverse shear strain γ<sup>i</sup> , respectively are defined as

$$\begin{aligned} \begin{bmatrix} \varepsilon\_{\rm xxi} \\ \varepsilon\_{\rm xyi} \\ \varepsilon\_{\rm yyi} \end{bmatrix} = \begin{bmatrix} \frac{\partial \mathbf{u}\_{\rm i}^{0}}{\partial \mathbf{x}} \\ \frac{\partial \mathbf{v}\_{\rm i}^{0}}{\partial \mathbf{y}} \\ \frac{\partial \mathbf{v}\_{\rm i}^{0}}{\partial \mathbf{y}} \\ \frac{\partial \mathbf{u}\_{\rm i}^{0}}{\partial \mathbf{y}} + \frac{\partial \mathbf{v}\_{\rm i}^{0}}{\partial \mathbf{x}} \end{bmatrix}, \begin{bmatrix} \kappa\_{\rm xxi} \\ \kappa\_{\rm xyi} \\ \kappa\_{\rm yyi} \\ \kappa\_{\rm xyi} \end{bmatrix} = \begin{bmatrix} \frac{\partial \mathbf{v}\_{\rm i}^{0}}{\partial \mathbf{x}} \\ \frac{\partial \mathbf{v}\_{\rm i}^{0}}{\partial \mathbf{y}} \\ \frac{\partial \mathbf{v}\_{\rm i}^{0}}{\partial \mathbf{y}} \\ \frac{\partial \mathbf{v}\_{\rm i}^{0}}{\partial \mathbf{y}} + \frac{\partial \mathbf{v}\_{\rm i}^{0}}{\partial \mathbf{x}} \end{bmatrix}, \begin{bmatrix} \gamma\_{\rm y} \mathbf{z}i \\ \nu\_{\rm xxi} \end{bmatrix} = \begin{bmatrix} \nu\_{\rm xi}^{0} + \frac{\partial \mathbf{w}\_{\rm i}^{0}}{\partial \mathbf{y}} \\ \nu\_{\rm xi}^{0} + \frac{\partial \mathbf{w}\_{\rm i}^{0}}{\partial \mathbf{x}} \\ \nu\_{\rm xi}^{0} + \frac{\partial \mathbf{w}\_{\rm i}^{0}}{\partial \mathbf{x}} \end{bmatrix} \tag{2} \end{aligned} \tag{2}$$

The constitutive relations between stress resultants and corresponding strains are given in (Reddy & Miravete, 1995; Žmindák, 2010). In these works standard FSDT finite elements available in ANSYS software are used (ANSYS, 2007) to join these elements at the interfaces inside each sublaminate using CE or rigid links characterized by two nodes and three degrees of freedom at each node.

#### *2.1.1. Interface elements for delamination modeling*

Delamination is defined as the fracture of the plane separating two plies of a laminated composite structure (Fenske, et al., 2001). This fracture occurs within the thin resin-rich layer that forms between plies during the manufacturing process. Perfect adhesion is assumed in the undelaminated region Ω Ω<sup>D</sup> , whereas sub-laminates are free to deflect along the delaminated region ΩD but not to penetrate each other. A linear interface model, is introduced along Ω Ω<sup>D</sup> to enforce adhesion. The constitutive equation of the interface involves two stiffness parameters, *kz*,*kxy*, imposing displacement continuity in the thickness and in-plane directions, respectively, by treating them as penalty parameters. The relationship between the components of the traction vector σ acting at the lower surface of the upper sublaminate, *σzx, σzy* and *σzz*, in the out-of-plane (*z*) and in the in-plane (*x* and *y*) directions, respectively, and the corresponding components of relative interface displacement vector Δ , *Δu*, *Δv* and *Δw* is expressed as

$$\mathbf{o} = \mathbf{K} \,\,\Delta\tag{3}$$

Interface elements are implemented using COMBIN14 element. Relative opening and sliding displacements are evaluated as the difference between displacements at the interface between the lower and the upper sublaminate.

#### *2.1.2. Contact formulation for damage interface*

136 Composites and Their Properties

normals about *y* and *x*, respectively.

degrees of freedom at each node.

*2.1.1. Interface elements for delamination modeling* 

displacement vector Δ , *Δu*, *Δv* and *Δw* is expressed as

between the lower and the upper sublaminate.

where *ui*, *vi* refer to the in-plane displacements, and *wi* to the transverse displacements through the thickness of the *i-*th plate element, 00 0 , , *ii i uvw* , are the displacements at the mid-

At the reference surfaces, the membrane strain vector i ε , the curvature i κ , and the

<sup>x</sup> <sup>0</sup> <sup>w</sup> <sup>0</sup> <sup>i</sup> xxi <sup>0</sup> xxi

<sup>v</sup> <sup>γ</sup> xi <sup>i</sup> yi yzi <sup>y</sup> <sup>ε</sup> , <sup>κ</sup> , yyi yyi y y <sup>0</sup> <sup>γ</sup>xzi <sup>0</sup> wi γ κ 0 0 <sup>0</sup> xyi u v xyi <sup>0</sup> xi i i yi <sup>x</sup> xi

The constitutive relations between stress resultants and corresponding strains are given in (Reddy & Miravete, 1995; Žmindák, 2010). In these works standard FSDT finite elements available in ANSYS software are used (ANSYS, 2007) to join these elements at the interfaces inside each sublaminate using CE or rigid links characterized by two nodes and three

Delamination is defined as the fracture of the plane separating two plies of a laminated composite structure (Fenske, et al., 2001). This fracture occurs within the thin resin-rich layer that forms between plies during the manufacturing process. Perfect adhesion is assumed in the undelaminated region Ω Ω<sup>D</sup> , whereas sub-laminates are free to deflect along the delaminated region ΩD but not to penetrate each other. A linear interface model, is introduced along Ω Ω<sup>D</sup> to enforce adhesion. The constitutive equation of the interface involves two stiffness parameters, *kz*,*kxy*, imposing displacement continuity in the thickness and in-plane directions, respectively, by treating them as penalty parameters. The relationship between the components of the traction vector σ acting at the lower surface of the upper sublaminate, *σzx, σzy* and *σzz*, in the out-of-plane (*z*) and in the in-plane (*x* and *y*) directions, respectively, and the corresponding components of relative interface

Interface elements are implemented using COMBIN14 element. Relative opening and sliding displacements are evaluated as the difference between displacements at the interface

 

(x,y), x,y denote rotations of transverse

x 0

(2)

**σ K Δ** (3)

<sup>0</sup> 0u xi <sup>i</sup>

y x y x

surface of the *i-*th plate element and xi yi

transverse shear strain γ<sup>i</sup> , respectively are defined as

ε κ

In order to avoid interpenetration between delaminated sublaminates in the delaminated region Ω<sup>D</sup> , a unilateral frictionless contact interface can be introduced, characterized by a zero stiffness for opening relative displacements (*Δw≥0*) and a positive stiffness for closing relative displacements (*Δw≤0*), then the contact stress zz σ is

$$
\sigma\_{xx} = \frac{1}{2} \left( 1 - \text{sign}(\Delta \mathbf{w}) \right) \mathbf{k}\_Z \Delta \mathbf{w} \tag{4}
$$

where *kz* is the penalty number imposing contact constraint and sign is the signum function. A very large value for *kz* restricts sublaminate overlapping and simulates the contact condition. Unilateral contact conditions may be implemented in ANSYS using COMBIN39. This element is a unidirectional element with nonlinear constitutive relationships with appropriate specialization of the nonlinear constitutive law according to (4).

If we introduce a scalar damage variable *D* with the value of *1* for no adhesion and the value of 0 for perfect adhesion, we get a single extended interface model with constitutive law valid both for undelaminated Ω Ω<sup>D</sup> and delaminated ΩD areas. Consequently constitutive law can be expressed as

$$
\sigma = \left(1 - D\right) \mathbf{K} \Delta \tag{5}
$$

In this work we use the formulation via FEs, related to plate elements, interface elements and Lagrange multipliers. It is worth noting that in commercial FEA packages the Lagrange multipliers are represented by either CE or rigid links, whereas interface elements are implemented by the analyst using a combination of spring elements (COMBIN14) and CE.

#### *2.1.3. Mixed mode analysis*

In order to predict crack propagation in laminates for general loading conditions, ERR distributions along the delamination front are needed. Fracture mechanics assumes that delamination propagation is controlled by the critical ERR. Delamination grows on the region of the delamination front where the following condition is satisfied and is of the form

$$\left(\frac{\mathbf{G}\_{\mathrm{I}}\left(\mathbf{s}\right)}{\mathbf{G}\_{\mathrm{I}}^{\mathrm{c}}}\right)^{\alpha} + \left(\frac{\mathbf{G}\_{\mathrm{II}}\left(\mathbf{s}\right)}{\mathbf{G}\_{\mathrm{II}}^{\mathrm{c}}}\right)^{\beta} + \left(\frac{\mathbf{G}\_{\mathrm{III}}\left(\mathbf{s}\right)}{\mathbf{G}\_{\mathrm{III}}^{\mathrm{c}}}\right)^{\mathsf{Y}} = \mathbf{1} \tag{6}$$

where α, β and γ are mixed mode fracture parameters determined by fitting experimental test results.

The critical ERR cc c , G ,G G material properties can be I II III evaluated from experimental procedures. The closed-form expressions for the ERR are (Barbero, 2008)

$$\mathbf{G(s)} = \frac{1}{2} \lim\_{\mathbf{k}\_{\mathbf{Z}}, \mathbf{k}\_{\mathbf{XY}} \to \infty} \left[ \mathbf{k}\_{\mathbf{Z}} \Delta \mathbf{w}^{2} (\mathbf{s}) + \mathbf{k}\_{\mathbf{xy}} \Delta \mathbf{u}^{2} (\mathbf{s}) + \mathbf{k}\_{\mathbf{xy}} \Delta \mathbf{v}^{2} (\mathbf{s}) \right] \Delta \mathbf{w} (\mathbf{s}) \ge 0 \tag{7}$$

$$\mathbf{G} = \mathbf{G}\_{\mathrm{I}}(\mathbf{s}) + \mathbf{G}\_{\mathrm{II}}(\mathbf{s}) + \mathbf{G}\_{\mathrm{III}}(\mathbf{s})$$

$$\mathbf{G}\_{\mathrm{I}}(\mathbf{s}) = \begin{vmatrix} \mathbf{k}\_{\mathrm{Z}}, \lim\_{\mathbf{x} \mathbf{y} \to \infty} \frac{1}{2} \mathbf{k}\_{\mathrm{Z}} \Delta \mathbf{w}^{2}(\mathbf{s}) & \text{if} & \Delta \mathbf{w}(\mathbf{s}) \ge 0 \\ 0 & \text{if} & \Delta \mathbf{w}(\mathbf{s}) < 0 \end{vmatrix} \tag{8}$$

$$\mathbf{G}\_{\mathrm{II}}(\mathbf{s}) = \lim\_{\mathbf{k}\_{\mathrm{Z}}, \mathbf{k}\_{\mathrm{XY}} \to \infty} \frac{1}{2} \mathbf{k}\_{\mathrm{XY}} \Delta \mathbf{u}\_{\mathrm{II}}^{2}(\mathbf{s})$$

$$\mathbf{G}\_{\mathrm{III}}(\mathbf{s}) = \lim\_{\mathbf{k}\_{\mathrm{Z}}, \mathbf{k}\_{\mathrm{XY}} \to \infty} \frac{1}{2} \mathbf{k}\_{\mathrm{XY}} \Delta \mathbf{u}\_{\mathrm{I}}^{2}(\mathbf{s})$$

Finite Element Implementation of Failure and Damage Simulation in Composite Plates 139

<sup>t</sup> <sup>R</sup> <sup>Δ</sup>w R <sup>Δ</sup><sup>u</sup> <sup>R</sup> <sup>Δ</sup><sup>u</sup> ´ ´ <sup>A</sup> ´ B B nB B tB B G A GA G A II III Δ Δ nn n Δ Δ Δ Δ t t 

where <sup>z</sup> RA is the reaction in the spring element connecting node A in the z-direction, ΔwB B is the relative z-displacement between the nodes *B* and *B'*. These are located immediately ahead of the delamination front along its normal direction passing through

Similar definitions apply for reactions and relative displacement related to modes II and III. The characteristic mesh sizes in the normal and tangential directions of the delamination front are denoted by Δn and Δ<sup>t</sup> . In (9), the same element size is assumed for elements ahead of and behind the delamination front. Value of Δ<sup>t</sup> /2 must be used in (9) instead of

In order to simplify the FE modeling procedure, it is possible to introduce spring elements only along the delamination front instead of the entire undelaminated region. The perfect adhesion along the remaining portion of the undelaminated region can be imposed by CE. However, when the delamination propagation must be simulated, it is necessary to

In the next example, the delamination modeling techniques presented so far are applied to analyze typical 3D delamination problems in laminated plates. The ERR distribution along

One of the most powerful computational methods for structural analysis of composites is the FEM. The starting point would be a "validated" FE model, with a reasonably fine mesh,

the delamination front are computed for different laminates and loading conditions.

11 1 , , 22 2 (9)

z n A A

**Figure 3.** Plate assembly in the neighborhood of the delamination front

introduce interface elements in the whole undelaminated region Ω Ω<sup>D</sup> .

**2.2. Finite element modeling and numerical example** 

Δt when the node is placed at a free edge.

t

I

*A*.

They are obtained by means of the interface model, using FE code to check whether propagation occurs. Once made a global FEA of the laminate, then the calculation of *G(s)* along the delamination front reduces to a simple post-computation. The extent of the propagation of the delamination area may be established by releasing the node in which the relation (6) is first satisfied, leading to a modification of the delamination front, which in turn requires another equilibrium solution. It follows the fact that the delamination growth analysis must be accomplished iteratively. For simplicity, only the computation of ERR is described here. The study of the propagation for a 3D planar delamination requires the use of nonlinear incremental numerical computation.

The delaminated laminate is represented using two sublaminates (Fig. 2). In this case, the model is called a two-layer plate model. Multilayer plate model in each sublaminate is necessary to achieve sufficient accuracy when the mode components are needed. Sublaminates are modeled using standard shear deformable elements (SHELL181), whereas interface elements can be used for the interface model. Available interface elements (INTER204) are only compatible with solid elements, therefore interface elements are simulated here by coupling CE with spring elements (COMBIN14). Plate and interface models must be described by the same in-plane mesh.

The FE model of the plates adjacent to the delamination plane in proximity of the delamination front is illustrated in Fig. 3. Interface elements model the undelaminated region Ω Ω<sup>D</sup> up to the delamination front. The mesh of interface and plate elements must be sufficiently refined in order to capture the high interface stress gradient in the neighborhood of the delamination front, which occurs because high values for interface stiffness must be used to simulate perfect adhesion. The individual ERR at the general node A of the delamination front are calculated using the reactions obtained from spring elements and the relative displacements between the nodes already delaminated and located along the normal direction.

ERRs are computed by using (9), which is a modified version of (8) in order to avoid excessive mesh refining at the delamination front. This leads to the following expressions

Finite Element Implementation of Failure and Damage Simulation in Composite Plates 139

$$\mathbf{G}\_{\rm I}(\mathbf{A}) = \left(\frac{1}{2} \frac{\mathbf{R}\_{\rm A}^{z} \,\Delta \mathbf{w}\_{\rm B} \,\mathbf{B}}{\Delta\_{\rm II} \,\Delta\_{\rm t}}\right) \prime \,\mathbf{G}\_{\rm II}(\mathbf{A}) = \left(\frac{1}{2} \frac{\mathbf{R}\_{\rm A}^{n} \,\Delta \mathbf{u}\_{\rm B} \,\mathbf{B}}{\Delta\_{\rm II} \,\Delta\_{\rm t}}\right) \prime \,\mathbf{G}\_{\rm III}(\mathbf{A}) = \left(\frac{1}{2} \frac{\mathbf{R}\_{\rm A}^{\rm t} \,\Delta \mathbf{u}\_{\rm B} \,\mathbf{B}}{\Delta\_{\rm II} \,\Delta\_{\rm t}}\right) \tag{9}$$

where <sup>z</sup> RA is the reaction in the spring element connecting node A in the z-direction, ΔwB B is the relative z-displacement between the nodes *B* and *B'*. These are located immediately ahead of the delamination front along its normal direction passing through *A*.

**Figure 3.** Plate assembly in the neighborhood of the delamination front

138 Composites and Their Properties

 <sup>1</sup> <sup>222</sup> G s lim k <sup>Δ</sup>ws k <sup>Δ</sup>us k <sup>Δ</sup>v s,Δws 0 z xy xy <sup>2</sup> k ,k <sup>z</sup> xy

<sup>1</sup> <sup>2</sup> lim <sup>k</sup> <sup>Δ</sup>w s if Δws 0 <sup>z</sup> G s k ,k <sup>2</sup> <sup>z</sup> xy <sup>I</sup>

GGs G s G s I II III

of nonlinear incremental numerical computation.

models must be described by the same in-plane mesh.

nodes already delaminated and located along the normal direction.

<sup>1</sup> <sup>2</sup> G s lim <sup>k</sup> <sup>Δ</sup>u s II xy n k ,k <sup>2</sup> <sup>z</sup> xy

<sup>1</sup> <sup>2</sup> G s lim <sup>k</sup> <sup>Δ</sup>u s III xy <sup>t</sup> k ,k <sup>2</sup> <sup>z</sup> xy

They are obtained by means of the interface model, using FE code to check whether propagation occurs. Once made a global FEA of the laminate, then the calculation of *G(s)* along the delamination front reduces to a simple post-computation. The extent of the propagation of the delamination area may be established by releasing the node in which the relation (6) is first satisfied, leading to a modification of the delamination front, which in turn requires another equilibrium solution. It follows the fact that the delamination growth analysis must be accomplished iteratively. For simplicity, only the computation of ERR is described here. The study of the propagation for a 3D planar delamination requires the use

The delaminated laminate is represented using two sublaminates (Fig. 2). In this case, the model is called a two-layer plate model. Multilayer plate model in each sublaminate is necessary to achieve sufficient accuracy when the mode components are needed. Sublaminates are modeled using standard shear deformable elements (SHELL181), whereas interface elements can be used for the interface model. Available interface elements (INTER204) are only compatible with solid elements, therefore interface elements are simulated here by coupling CE with spring elements (COMBIN14). Plate and interface

The FE model of the plates adjacent to the delamination plane in proximity of the delamination front is illustrated in Fig. 3. Interface elements model the undelaminated region Ω Ω<sup>D</sup> up to the delamination front. The mesh of interface and plate elements must be sufficiently refined in order to capture the high interface stress gradient in the neighborhood of the delamination front, which occurs because high values for interface stiffness must be used to simulate perfect adhesion. The individual ERR at the general node A of the delamination front are calculated using the reactions obtained from spring elements and the relative displacements between the

ERRs are computed by using (9), which is a modified version of (8) in order to avoid excessive mesh refining at the delamination front. This leads to the following expressions

0 if Δws 0

(7)

(8)

Similar definitions apply for reactions and relative displacement related to modes II and III. The characteristic mesh sizes in the normal and tangential directions of the delamination front are denoted by Δn and Δ<sup>t</sup> . In (9), the same element size is assumed for elements ahead of and behind the delamination front. Value of Δ<sup>t</sup> /2 must be used in (9) instead of Δt when the node is placed at a free edge.

In order to simplify the FE modeling procedure, it is possible to introduce spring elements only along the delamination front instead of the entire undelaminated region. The perfect adhesion along the remaining portion of the undelaminated region can be imposed by CE. However, when the delamination propagation must be simulated, it is necessary to introduce interface elements in the whole undelaminated region Ω Ω<sup>D</sup> .

In the next example, the delamination modeling techniques presented so far are applied to analyze typical 3D delamination problems in laminated plates. The ERR distribution along the delamination front are computed for different laminates and loading conditions.

#### **2.2. Finite element modeling and numerical example**

One of the most powerful computational methods for structural analysis of composites is the FEM. The starting point would be a "validated" FE model, with a reasonably fine mesh, correct boundary conditions, material properties, etc. (Bathe, 1996). As a minimum requirement, the model is expected to produce stress and strains that have reasonable accuracy to those of the real structure prior to failure initiation. In spite of the great success of the FEM and BEM as effective numerical tools for the solution of boundary-value problems on complex domains, there is still a growing interest in the development of new advanced methods. Many meshless formulations are becoming popular due to their high adaptivity and a low cost to prepare input data for numerical analysis (Guiamatsia et al., 2009).

Finite Element Implementation of Failure and Damage Simulation in Composite Plates 141

of the ERR the equation (9) was used for the type I. As the biggest ERR is in the middle of the model, it is expected that the beginning of the delamination is in the middle of the model. The

**Figure 5.** Scheme of FEM Model I: a) Forces applied on the laminate model, b) ERR distribution of *GI*

(a) (b)

In this model the delamination type II (sliding type)was simulated. The applied forces are parallel to the x axis, Fig. 6 a). Two types of ERR were analyzed in this model, GII in the x direction and GIII in the y direction. ERRs were calculated separately for each direction. The reaction of the spring elements are used for the calculation of GII and the y reactions are used for the calculation of GIII . Both distributions are displaying the absolute values of ERR, at both distributions the values of ERR are smaller then the values of *GI*. The values of *GII* are in the range of (0.5, 2). 10-4 and the values of *GIII* are in the range of (0, 4). 10-5 (Fig. 7)

At this model the delamination type III (tearing type) was analyzed. The geometry of model I was used, with the mesh and its refinement around the delamination front, boundary conditions and the linking between the shell plates, but the direction of the applied forces has changed, Fig.8a). Both types of ERR are analyzed here, GII in the x direction and GIII the y direction. The values of ERR are in these ranges: value of *GII* in the range of ( 0, 14).10-3 and the value of *GIII* in the range of (0, 0.02). It is possible that better results could be achieved by increasing of the number of plate elements layers simulating the sublamina. These models can be also modeled by solid elements, but there is greater number of elements needed for accurately simulating of the stress and ERR gradients. Thereby the

number of equations and computing time increase.

for delamination type I

*Model II* 

*Model III* 

distribution of ERR through the width of the laminate is displayed on the Fig.5b.

The results of the delamination analysis of two laminae with different thickness and material are processed in this section. The laminae are fixed on one side and free on the other side. Loads are applied on the free side depending on the analyzed type of delamination.

The upper sublamina has these properties: *Ex*= 35000 MPa, *Ey* = *Ez* = 10500 MPa, *Gyz* = 10500 MPa, *Gxy* = *Gxz* = 1167 MPa, *v*xy = *v*xz = *v*yz= 0.3. The lower sublamina has these properties: *Ex=* 70000 MPa, *Ey = Ez* = 21000 MPa, *Gyz* = 2100 MPa, *Gxy = Gxz* = 2333 MPa, *v*xy = *v*xz = *v*yz= 0.3. The pair of forces applied on the laminae is *T*= 1 N/mm and the dimensions of the laminate are: a = 10mm, B = 20 mm , *L*= 20mm, *h1*= 0.5 mm, *h2*= 1mm.

**Figure 4.** a) Scheme of boundary conditions and laminate dimensions, b) scheme of FEM model

The upper sublaminate is composed by four plates *nu*= 4 and the lower by two *nl* = 2. The plates are meshed by SHELL181 elements. The zone of mesh refinement has these dimensions 5 \* 20 mm, it is centered around of the delamination front which is placed in the middle of the laminate. The interface between the sublaminates is modeled without stiffness for opening displacements and with positive stiffness for closing displacements. The interface between sublaminates is modeled by means of CE (Constrain Equation), since it is easier to apply than beam elements and delamination propagation is not solved. The delamination front is created by spring elements COMBIN 14, in each node of the delamination front by three elements. The stiffness of the spring elements binding the laminae is chosen as *kz* a *kxy* = 108 N/mm3. These elements are oriented in different directions, they were created always from a pair of nodes placed on the surface of the lower sublamina. One of the pair nodes is bounded to the upper plate by means of CE and the second one is bounded to the lower plate. ERR is calculated by using deformation along the delamination front.

#### *Model I*

At model I the ERR for delamination type I is analyzed. Model I is loaded with opening forces *T* of magnitude 1 N/mm*,* which are parallel to z axis, displayed on Fig. 5a. For the calculation of the ERR the equation (9) was used for the type I. As the biggest ERR is in the middle of the model, it is expected that the beginning of the delamination is in the middle of the model. The distribution of ERR through the width of the laminate is displayed on the Fig.5b.

**Figure 5.** Scheme of FEM Model I: a) Forces applied on the laminate model, b) ERR distribution of *GI* for delamination type I

#### *Model II*

140 Composites and Their Properties

b) scheme of FEM model

*Model I* 

correct boundary conditions, material properties, etc. (Bathe, 1996). As a minimum requirement, the model is expected to produce stress and strains that have reasonable accuracy to those of the real structure prior to failure initiation. In spite of the great success of the FEM and BEM as effective numerical tools for the solution of boundary-value problems on complex domains, there is still a growing interest in the development of new advanced methods. Many meshless formulations are becoming popular due to their high adaptivity

The results of the delamination analysis of two laminae with different thickness and material are processed in this section. The laminae are fixed on one side and free on the other side.

The upper sublamina has these properties: *Ex*= 35000 MPa, *Ey* = *Ez* = 10500 MPa, *Gyz* = 10500 MPa, *Gxy* = *Gxz* = 1167 MPa, *v*xy = *v*xz = *v*yz= 0.3. The lower sublamina has these properties: *Ex=* 70000 MPa, *Ey = Ez* = 21000 MPa, *Gyz* = 2100 MPa, *Gxy = Gxz* = 2333 MPa, *v*xy = *v*xz = *v*yz= 0.3. The pair of forces applied on the laminae is *T*= 1 N/mm and the dimensions of the laminate are: a

The upper sublaminate is composed by four plates *nu*= 4 and the lower by two *nl* = 2. The plates are meshed by SHELL181 elements. The zone of mesh refinement has these dimensions 5 \* 20 mm, it is centered around of the delamination front which is placed in the middle of the laminate. The interface between the sublaminates is modeled without stiffness for opening displacements and with positive stiffness for closing displacements. The interface between sublaminates is modeled by means of CE (Constrain Equation), since it is easier to apply than beam elements and delamination propagation is not solved. The delamination front is created by spring elements COMBIN 14, in each node of the delamination front by three elements. The stiffness of the spring elements binding the laminae is chosen as *kz* a *kxy* = 108 N/mm3. These elements are oriented in different directions, they were created always from a pair of nodes placed on the surface of the lower sublamina. One of the pair nodes is bounded to the upper plate by means of CE and the second one is bounded to the lower

At model I the ERR for delamination type I is analyzed. Model I is loaded with opening forces *T* of magnitude 1 N/mm*,* which are parallel to z axis, displayed on Fig. 5a. For the calculation

and a low cost to prepare input data for numerical analysis (Guiamatsia et al., 2009).

Loads are applied on the free side depending on the analyzed type of delamination.

= 10mm, B = 20 mm , *L*= 20mm, *h1*= 0.5 mm, *h2*= 1mm.

**Figure 4.** a) Scheme of boundary conditions and laminate dimensions,

plate. ERR is calculated by using deformation along the delamination front.

In this model the delamination type II (sliding type)was simulated. The applied forces are parallel to the x axis, Fig. 6 a). Two types of ERR were analyzed in this model, GII in the x direction and GIII in the y direction. ERRs were calculated separately for each direction. The reaction of the spring elements are used for the calculation of GII and the y reactions are used for the calculation of GIII . Both distributions are displaying the absolute values of ERR, at both distributions the values of ERR are smaller then the values of *GI*. The values of *GII* are in the range of (0.5, 2). 10-4 and the values of *GIII* are in the range of (0, 4). 10-5 (Fig. 7)

#### *Model III*

At this model the delamination type III (tearing type) was analyzed. The geometry of model I was used, with the mesh and its refinement around the delamination front, boundary conditions and the linking between the shell plates, but the direction of the applied forces has changed, Fig.8a). Both types of ERR are analyzed here, GII in the x direction and GIII the y direction. The values of ERR are in these ranges: value of *GII* in the range of ( 0, 14).10-3 and the value of *GIII* in the range of (0, 0.02). It is possible that better results could be achieved by increasing of the number of plate elements layers simulating the sublamina. These models can be also modeled by solid elements, but there is greater number of elements needed for accurately simulating of the stress and ERR gradients. Thereby the number of equations and computing time increase.

Finite Element Implementation of Failure and Damage Simulation in Composite Plates 143

S for resistance can be found. Total area of

S (11)

DS SS (10)

DD n n (12)

There are many material modeling strategies to predict damage in laminated composites subjected to static or impulsive loads. Broadly, they can be classified as (Jain & Ghosh,

We consider a volume of material free of damage if no cracks or cavities can be observed at the microscopic scale. The opposite state is the fracture of the volume element. Theory of damage describes the phenomena between the virgin state of material and the macroscopic onset of crack (Jain & Ghosh, 2009; Tumino et al., 2007). The volume element must be of sufficiently large size compared to the inhomogenities of the composite material. In Fig. 9 this volume is depicted. One section of this element is related to its normal and to its area *S*.

DS <sup>D</sup>

We note that damage *D* is a scalar assuming values between 0 and 1. For *D* = 0 the material is undamaged, for 0<*D*<1 the material is damaged, for *D* = 1 complete failure occurs. The quantitative evaluation of damage is not a trivial issue, it must be linked to a variable that is

For isotropic damage, the dependence on the normal n can be neglected, i.e.

**3. Continuum damage mechanics** 

 plasticity or yield surface approach, damage mechanics approach

failure criteria approach (Kormaníková, 2011),

Due to the presence of defects, an effective area

able to characterize the phenomenon.

The local damage related to the direction **n** is defined as:

**Figure 9.** Representative volume element for damage mechanics

fracture mechanics approach (based on energy release rates),

2009):

defects is therefore:

**Figure 6.** Scheme of FEM Model II: a) direction of loads for delamination type II, b) distribution of ERR for type II delamination of *GII*

**Figure 7.** ERR distribution of *GIII* for delamination type II

**Figure 8.** Scheme of FEM Model III: a) definitions of loads for delamination type III, b) distribution of GII for model III

#### **3. Continuum damage mechanics**

142 Composites and Their Properties

for type II delamination of *GII*

GII for model III

**Figure 7.** ERR distribution of *GIII* for delamination type II

**Figure 6.** Scheme of FEM Model II: a) direction of loads for delamination type II, b) distribution of ERR

(a) (b)

**Figure 8.** Scheme of FEM Model III: a) definitions of loads for delamination type III, b) distribution of

(a) (b)

There are many material modeling strategies to predict damage in laminated composites subjected to static or impulsive loads. Broadly, they can be classified as (Jain & Ghosh, 2009):


We consider a volume of material free of damage if no cracks or cavities can be observed at the microscopic scale. The opposite state is the fracture of the volume element. Theory of damage describes the phenomena between the virgin state of material and the macroscopic onset of crack (Jain & Ghosh, 2009; Tumino et al., 2007). The volume element must be of sufficiently large size compared to the inhomogenities of the composite material. In Fig. 9 this volume is depicted. One section of this element is related to its normal and to its area *S*. Due to the presence of defects, an effective area S for resistance can be found. Total area of defects is therefore:

$$\mathbf{S}\_{\rm D} = \mathbf{S} - \tilde{\mathbf{S}} \tag{10}$$

The local damage related to the direction **n** is defined as:

$$\mathbf{D} = \frac{\mathbf{S}\_{\text{D}}}{\mathbf{S}} \tag{11}$$

For isotropic damage, the dependence on the normal n can be neglected, i.e.

$$\mathbf{D} = \mathbf{D}\_{\mathbf{II}} \,\,\forall \mathbf{n} \tag{12}$$

We note that damage *D* is a scalar assuming values between 0 and 1. For *D* = 0 the material is undamaged, for 0<*D*<1 the material is damaged, for *D* = 1 complete failure occurs. The quantitative evaluation of damage is not a trivial issue, it must be linked to a variable that is able to characterize the phenomenon.

**Figure 9.** Representative volume element for damage mechanics

We note that several papers can be found in literature where the constitutive equations of the materials are a function of a scalar variable of damage (Barbero, 2008). For the formulation of a general multidimensional damage model it is necessary to generalize the scalar damage variables. It is therefore necessary to define corresponding tensorial damage variables that can be used in general states of deformation and damage (Tumino at al., 2007). Finite Element Implementation of Failure and Damage Simulation in Composite Plates 145

(16)

(17)

Equation (13) and (14) can be written for different simple stress states: tension and compression in fiber direction, tension in transverse direction, in-plane shear. Tensors *J* and

In the present model for damage isotropic hardening is considered and the hardening

The hardening parameters γ0, c1 and c2 are determined by approximating the experimental stress-strain curves for in-plane shear loading. If this curve is not available, we can

> <sup>σ</sup> F tanh <sup>γ</sup> <sup>F</sup>

where *F6* is the in-plane shear strength, *G12* is the in-plane initial elasticity modulus and γ6 is the in-plane shear strain (in the lamina coordinate system). This function represents

The reaching of critical damage level is dependent on stresses in points of lamina. If in a point of lamina only normal stresses in the fiber direction or transverse direction (i.e., normal stress in lamina coordinate system) occur, then simply comparing the values of damaged variables with critical values of damage variables for given material at this point is sufficient. The damage has reached critical level if at least one of the values *D1*, *D2* in the point of lamina is greater or equal to its critical value. If in a given point of lamina also shear stress occurs (in lamina coordinate system), it is additionally necessary to compare the value of the product of (1 - *D1*) (1 - *D2*) with *ks* value for given material. If the value of this product is less or equal to *ks* value, the damage has reached a critical

The Newton-Raphson method was used for solving the system of nonlinear equations. Evolution of damage has been solved using return-mapping algorithm described in (Neto, 2008). The input values are strains and strain increments in lamina coordinate system, state variables *D1*, *D2*, and δ in integration point from the start of last performed iteration, **C** matrix (gives the stress-strain relations in the effective configuration (Barbero, 2007) and

2 δ γ c exp 1 c 

> 6 12

6

G

1

6 6

*H* can be derived in terms of material strength values.

*3.2.1. Hardening parameters* 

reconstruct it using function

experimental data very well.

*3.1.3. Critical damage level* 

**3.2. Implementation of numerical method** 

level.

function was used in the form of

In this part, we focused on presenting the methodology of numerical solving of elastic damage of thin composite plates reinforced by long fibers based on continuum damage mechanics by means of the finite element method.

#### **3.1. Damage model used**

The model for fiber-reinforced lamina mentioned next was presented by Barbero and de Vivo (Barbero, 2001) and is suitable for fiber - reinforced composite material with polymer matrix. On the lamina level these composites are considered as ideal homogenous and transversely isotropic. All parameters of this model can be easily identified from available experimental data. It is assumed that damage in principal directions is identical with the principal material directions throughout the damage process. Therefore the evolution of damage is solved in the lamina coordinate system. The model predicts the evolution of damage and its effect on stiffness and subsequent redistribution of stress.

#### *3.1.1. Damage surface and damage potential*

Damage surface is similar to the Tsai-Wu damage surface and is defined by tensors **J** and **H (**Barbero, 2008 and it is commonly used for predicting failure of fiber-reinforced lamina with respect to experimental material strength values. Damage surface and damage potential have the form of (Barbero & de Vivo, 2001)

$$\log\left(\mathbf{Y},\gamma\right) = \sqrt{\mathbf{J}\_{11}\mathbf{Y}\_1^2 + \mathbf{J}\_{22}\mathbf{Y}\_2^2 + \mathbf{J}\_{33}\mathbf{Y}\_3^2} + \sqrt{\mathbf{H}\_1\mathbf{Y}\_1^2 + \mathbf{H}\_2\mathbf{Y}\_2^2 + \mathbf{H}\_3\mathbf{Y}\_3^2} - \left(\gamma + \gamma\_0\right) \tag{13}$$

$$\text{If } \mathbf{f}\left(\mathbf{Y}, \boldsymbol{\gamma}\right) = \sqrt{\mathbf{J}\_{11}\mathbf{Y}\_1^2 + \mathbf{J}\_{22}\mathbf{Y}\_2^2 + \mathbf{J}\_{33}\mathbf{Y}\_3^2} - \left(\boldsymbol{\gamma} + \boldsymbol{\gamma}\_o\right) \tag{14}$$

where the thermodynamic forces *Y1*, *Y2* and *Y3* can be calculated by means of relations

$$\begin{aligned} \mathbf{Y}\_1 &= \frac{1}{\Omega\_1^2} \left( \frac{\overline{\mathbf{S}}\_{11}}{\Omega\_1^4} \sigma\_1^2 + \frac{\overline{\mathbf{S}}\_{12}}{\Omega\_1^2} \sigma\_2 \sigma\_1 \sigma\_2 + \frac{\overline{\mathbf{S}}\_{66}}{\Omega\_1^2} \sigma\_2^2 \right) \\ \mathbf{Y}\_2 &= \frac{1}{\Omega\_2^2} \left( \frac{\overline{\mathbf{S}}\_{22}}{\Omega\_2^4} \sigma\_2^2 + \frac{\overline{\mathbf{S}}\_{12}}{\Omega\_1^2 \Omega\_2^2} \sigma\_1 \sigma\_2 + \frac{\overline{\mathbf{S}}\_{66}}{\Omega\_1^2 \Omega\_2^2} \sigma\_6^2 \right) \\ \mathbf{Y}\_3 &= 0 \end{aligned} \tag{15}$$

where stresses and components of matrix **S** are defined in the lamina coordinate system. Matrix **S** gives the strain-stress relations in the effective configuration (Barbero, 2007). Equation (13) and (14) can be written for different simple stress states: tension and compression in fiber direction, tension in transverse direction, in-plane shear. Tensors *J* and *H* can be derived in terms of material strength values.

#### *3.2.1. Hardening parameters*

144 Composites and Their Properties

**3.1. Damage model used** 

mechanics by means of the finite element method.

*3.1.1. Damage surface and damage potential* 

have the form of (Barbero & de Vivo, 2001)

We note that several papers can be found in literature where the constitutive equations of the materials are a function of a scalar variable of damage (Barbero, 2008). For the formulation of a general multidimensional damage model it is necessary to generalize the scalar damage variables. It is therefore necessary to define corresponding tensorial damage variables that can be used in general states of deformation and damage (Tumino at al., 2007).

In this part, we focused on presenting the methodology of numerical solving of elastic damage of thin composite plates reinforced by long fibers based on continuum damage

The model for fiber-reinforced lamina mentioned next was presented by Barbero and de Vivo (Barbero, 2001) and is suitable for fiber - reinforced composite material with polymer matrix. On the lamina level these composites are considered as ideal homogenous and transversely isotropic. All parameters of this model can be easily identified from available experimental data. It is assumed that damage in principal directions is identical with the principal material directions throughout the damage process. Therefore the evolution of damage is solved in the lamina coordinate system. The model predicts the evolution of

Damage surface is similar to the Tsai-Wu damage surface and is defined by tensors **J** and **H (**Barbero, 2008 and it is commonly used for predicting failure of fiber-reinforced lamina with respect to experimental material strength values. Damage surface and damage potential

222 2 2 2

 22 33 <sup>0</sup> 222

> 11 2 2 12 66 2 4 22 22 1 1 2 6 1 1 12 12

 

22 12 2 2 66

2 4 22 22 2 1 2 6 2 12 12

<sup>1</sup> S S <sup>S</sup> <sup>Y</sup> σ σσ σ <sup>Ω</sup><sup>2</sup> Ω ΩΩ ΩΩ <sup>2</sup> Y0

<sup>1</sup> S S <sup>S</sup> <sup>Y</sup> σ σσ σ <sup>1</sup> Ω Ω ΩΩ ΩΩ

3

where stresses and components of matrix **S** are defined in the lamina coordinate system. Matrix **S** gives the strain-stress relations in the effective configuration (Barbero, 2007).

where the thermodynamic forces *Y1*, *Y2* and *Y3* can be calculated by means of relations

11 1 22 2 33 3 1 1 2 2 3 3 g Y,<sup>γ</sup> J Y J Y J Y HY HY HY γ γ<sup>0</sup> (13)

11 1 2 3 f Y,γ JY JY JY γ γ (14)

(15)

damage and its effect on stiffness and subsequent redistribution of stress.

In the present model for damage isotropic hardening is considered and the hardening function was used in the form of

$$\gamma = \mathbf{c}\_1 \left[ \exp\left(\frac{\mathbf{s}}{\mathbf{c}\_2}\right) + 1 \right] \tag{16}$$

The hardening parameters γ0, c1 and c2 are determined by approximating the experimental stress-strain curves for in-plane shear loading. If this curve is not available, we can reconstruct it using function

$$\sigma\_{\boldsymbol{\theta}} = \mathbf{F}\_{\boldsymbol{\theta}} \tanh \left( \frac{\mathbf{G}\_{12}}{\mathbf{F}\_{\boldsymbol{\theta}}} \boldsymbol{\gamma}\_{\boldsymbol{\theta}} \right) \tag{17}$$

where *F6* is the in-plane shear strength, *G12* is the in-plane initial elasticity modulus and γ6 is the in-plane shear strain (in the lamina coordinate system). This function represents experimental data very well.

#### *3.1.3. Critical damage level*

The reaching of critical damage level is dependent on stresses in points of lamina. If in a point of lamina only normal stresses in the fiber direction or transverse direction (i.e., normal stress in lamina coordinate system) occur, then simply comparing the values of damaged variables with critical values of damage variables for given material at this point is sufficient. The damage has reached critical level if at least one of the values *D1*, *D2* in the point of lamina is greater or equal to its critical value. If in a given point of lamina also shear stress occurs (in lamina coordinate system), it is additionally necessary to compare the value of the product of (1 - *D1*) (1 - *D2*) with *ks* value for given material. If the value of this product is less or equal to *ks* value, the damage has reached a critical level.

#### **3.2. Implementation of numerical method**

The Newton-Raphson method was used for solving the system of nonlinear equations. Evolution of damage has been solved using return-mapping algorithm described in (Neto, 2008). The input values are strains and strain increments in lamina coordinate system, state variables *D1*, *D2*, and δ in integration point from the start of last performed iteration, **C** matrix (gives the stress-strain relations in the effective configuration (Barbero, 2007) and damage parameters related to damage model. The output variables are *D1*, *D2*, and δ, stresses and strains in lamina coordinate system in this integration point at the end of the last performed iteration. Another output is constitutive damage matrix **C**ED in lamina coordinate system, which reflects the effect of damage on the behavior of structure. Flowchart of this algorithm is described in Fig. 10.

Finite Element Implementation of Failure and Damage Simulation in Composite Plates 147

force *F* = - 4000 N in the middle of the plate. Own program created in MATLAB language was used for this analysis. Four-node layered plate finite elements based on Kirchhoff

Material properties, damage parameters and hardening parameters and critical damage values are given in Table 1 - Table 3. The parameters *J33* and *H3* are equal to zero. The plate model was divided into 8×8 elements and was analyzed in fifty load substeps. The linear static analysis shows that the largest magnitudes of stress are in parallel direction with fibers and transverse to fibers and they occur in the outer layers in the middle of the

**M30/949** 167 8.13 4.41 0.27 **M40/948** 228 7.99 4.97 0.292

**M30/949** 0.952.10-3 0.438 25.585.10-3 21.665.10-3 0.6 0.30 0.395 **M40/948** 2.208.10-3 0.214 10.503.10-3 8.130.10-3 0.12 0.10 0.395

**cr D1c**

**M30/949** 0.105 0.111 0.5 0.944 **M40/948** 0.105 0.111 0.5 0.908

The largest magnitudes of shear stress occur in the outer layers in the corner nodes. The largest magnitudes of stress in layers 2, 3, 6 and 7 occur in the center of the plate. According to the results of linear static analysis it can be expected that damage will reach the critical

Fig. 11 shows the analysis results of elastic damage of the plate made from material M30/949. Fig. 11a shows the evolution of individual stress components in dependence on strains (both in lamina coordinate system) in the midsurface of layer 1 (first layer from the bottom) in integration point (IP) 1 (in element 1, nearest to the corner). Fig. 11b shows the evolution of individual stress components in the midsurface of layer 2 in IP 872 (in element 28, nearest to the center of the plate). Fig. 12 plots described damage variables evolution in these IPs. The analysis results show that reaching the critical level is caused not by normal stresses in the lamina coordinate system, but by shear stress (in the lamina coordinate system). The analysis results of the plate made from material M30/949 show that for given load the critical level of damage was reached in layers 2 and 7 in the center of the plate and its vicinity. In IPs that are closest to the center of the plate in these layers, the critical level of damage was reached between 13th and 14th load substep. However, it is not postulated that

**E1**[GPa] **E2**[GPa] **G12[**GPa] **υ<sup>12</sup>**

**1J 2J H1 H2 <sup>0</sup> γ** <sup>1</sup> c <sup>2</sup> c

**cr D2c <sup>s</sup> <sup>k</sup>**

(classical) plate theory were used.

**Table 1.** Material properties

**Table 2.** Damage and hardening parameters

**Table 3.** Critical values of damage variables

level in some of the above points.

**cr D1t**

plate.

**Figure 10.** Flowchart of return-mapping algorithm used for solving damage evolution in particular integration points

## **3.3. Numerical example**

One problem for two different materials was simulated in order to study the damage of laminated fiber-reinforced composite plates. The composites are reinforced by carbon fibers embedded in epoxy matrix. The simply supported composite plate with laminate stacking sequence of [0, 45, -45, 90]S with dimensions of 125×125×2.5 mm was loaded by transverse force *F* = - 4000 N in the middle of the plate. Own program created in MATLAB language was used for this analysis. Four-node layered plate finite elements based on Kirchhoff (classical) plate theory were used.

Material properties, damage parameters and hardening parameters and critical damage values are given in Table 1 - Table 3. The parameters *J33* and *H3* are equal to zero. The plate model was divided into 8×8 elements and was analyzed in fifty load substeps. The linear static analysis shows that the largest magnitudes of stress are in parallel direction with fibers and transverse to fibers and they occur in the outer layers in the middle of the plate.


**Table 1.** Material properties

146 Composites and Their Properties

integration points

**3.3. Numerical example** 

Flowchart of this algorithm is described in Fig. 10.

damage parameters related to damage model. The output variables are *D1*, *D2*, and δ, stresses and strains in lamina coordinate system in this integration point at the end of the last performed iteration. Another output is constitutive damage matrix **C**ED in lamina coordinate system, which reflects the effect of damage on the behavior of structure.

**Figure 10.** Flowchart of return-mapping algorithm used for solving damage evolution in particular

One problem for two different materials was simulated in order to study the damage of laminated fiber-reinforced composite plates. The composites are reinforced by carbon fibers embedded in epoxy matrix. The simply supported composite plate with laminate stacking sequence of [0, 45, -45, 90]S with dimensions of 125×125×2.5 mm was loaded by transverse


**Table 2.** Damage and hardening parameters


**Table 3.** Critical values of damage variables

The largest magnitudes of shear stress occur in the outer layers in the corner nodes. The largest magnitudes of stress in layers 2, 3, 6 and 7 occur in the center of the plate. According to the results of linear static analysis it can be expected that damage will reach the critical level in some of the above points.

Fig. 11 shows the analysis results of elastic damage of the plate made from material M30/949. Fig. 11a shows the evolution of individual stress components in dependence on strains (both in lamina coordinate system) in the midsurface of layer 1 (first layer from the bottom) in integration point (IP) 1 (in element 1, nearest to the corner). Fig. 11b shows the evolution of individual stress components in the midsurface of layer 2 in IP 872 (in element 28, nearest to the center of the plate). Fig. 12 plots described damage variables evolution in these IPs. The analysis results show that reaching the critical level is caused not by normal stresses in the lamina coordinate system, but by shear stress (in the lamina coordinate system). The analysis results of the plate made from material M30/949 show that for given load the critical level of damage was reached in layers 2 and 7 in the center of the plate and its vicinity. In IPs that are closest to the center of the plate in these layers, the critical level of damage was reached between 13th and 14th load substep. However, it is not postulated that used damage model predicts failure: it only predicts damage evolution and its effect on stiffness and consequent stress redistribution (Barbero & de Vivo, 2001). In some cases, failure can occur before the critical level of damage is reached. For plate made from material M30/949 load F*=*-1096 N is already critical.

Finite Element Implementation of Failure and Damage Simulation in Composite Plates 149

between 12th and 13th load substep in the nearest IPs. The critical level of damage would be also reached in the center of the plate and its vicinity in layers 2 and 7 (in the nearest IPs it would be reached between 16th and 17th load substep) and also in layers 3 and 6 (in the nearest IPs it would be reached between 27th and 28th load substep). Fig. 14 shows damage variables evolution in IP 868 and IP 872. For plate made from material M40/948 load F *= -*

**Figure 13.** Stress and strain evolution for plate from material M40/948 in (a) IP 868, (b) IP 872

(a) (b)

**Figure 14.** Damage variables evolution for plate from material M40/948 in (a) IP 868, (b) IP 872

The methodology of delamination calculation in laminated plates was applied in this chapter. The analyses shows that if mixed mode conditions are involved, a double plate model is suitable to accurately capture the mode decomposition in region near the midpoint

(a) (b)

990 N is already critical.

**4. Conclusion** 

**Figure 11.** Stress and strain evolution for plate from material M30/949 in (a) IP 1, (b) IP 872

**Figure 12.** Damage variables evolution for plate from material M30/949 in (a) IP 1, (b) IP 872

Fig. 13 shows the analysis results of elastic damage of the plate made from material M40/948. Fig. 13a shows the evolution of individual stress components in dependence on strains (both in lamina coordinate system) in the midsurface of layer 1 in IP 868 (in element 28, nearest to the center of the plate). Fig. 13b shows the evolution of individual stress components in the midsurface of layer 2 in IP 872 (in element 28, nearest to the center of the plate). The results show that reaching the critical level of damage will be also caused by shear stress in lamina coordinate system. However, the critical damage level was reached in layers 1 and 7 in the center of plate and its vicinity. The critical level of damage was reached between 12th and 13th load substep in the nearest IPs. The critical level of damage would be also reached in the center of the plate and its vicinity in layers 2 and 7 (in the nearest IPs it would be reached between 16th and 17th load substep) and also in layers 3 and 6 (in the nearest IPs it would be reached between 27th and 28th load substep). Fig. 14 shows damage variables evolution in IP 868 and IP 872. For plate made from material M40/948 load F *= -* 990 N is already critical.

**Figure 13.** Stress and strain evolution for plate from material M40/948 in (a) IP 868, (b) IP 872

**Figure 14.** Damage variables evolution for plate from material M40/948 in (a) IP 868, (b) IP 872

#### **4. Conclusion**

148 Composites and Their Properties

M30/949 load F*=*-1096 N is already critical.

used damage model predicts failure: it only predicts damage evolution and its effect on stiffness and consequent stress redistribution (Barbero & de Vivo, 2001). In some cases, failure can occur before the critical level of damage is reached. For plate made from material

**Figure 11.** Stress and strain evolution for plate from material M30/949 in (a) IP 1, (b) IP 872

(a) (b)

**Figure 12.** Damage variables evolution for plate from material M30/949 in (a) IP 1, (b) IP 872

Fig. 13 shows the analysis results of elastic damage of the plate made from material M40/948. Fig. 13a shows the evolution of individual stress components in dependence on strains (both in lamina coordinate system) in the midsurface of layer 1 in IP 868 (in element 28, nearest to the center of the plate). Fig. 13b shows the evolution of individual stress components in the midsurface of layer 2 in IP 872 (in element 28, nearest to the center of the plate). The results show that reaching the critical level of damage will be also caused by shear stress in lamina coordinate system. However, the critical damage level was reached in layers 1 and 7 in the center of plate and its vicinity. The critical level of damage was reached

(a) (b)

The methodology of delamination calculation in laminated plates was applied in this chapter. The analyses shows that if mixed mode conditions are involved, a double plate model is suitable to accurately capture the mode decomposition in region near the midpoint of delamination front. The solution converges quickly because a small number of plates is needed to obtain a reasonable approximation

Finite Element Implementation of Failure and Damage Simulation in Composite Plates 151

Clegg, R. A. et al. (2006). Hypervelocity impact damage prediction in composites: Part Imaterial model and characterisation. *International Journal of Impact Engineering.* Vol. 33,

Cui W, Wisnom M. A. (1993). Combined stress-based and fracture mechanics-based model for predicting delamination in composites. *Composites* 1993; Vol. 24, pp. 467–74, ISSN

Dávila, C.G., Camanho, P.P., Rose, Ch. A. (2005), Failure criteria for FRP Laminates*. Journal* 

Elmarakbi, A., Hu, N., Fukunaga, H., (2009), Finite element simulation of delamination growth in composite materials using LS-DYNA. *Composites Science & Technology*, 69 pp.

Fenske, M.T., Vizzini, A.J. (2001), The inclusion of in-plane stresses in delamination Criteria.

Gayathri, P. et al. (2010). Effect of matrix cracking and material uncertainty on composite plates. *Reliability Engineering* & *System Safety*. Vol. 95, No. 7, pp. 716-728, ISSN 0951-8320 Guiamatsia, I. et al. (2009). Element-free Galerkin modelling of composite damage. *Composites Science and Technology.* Vol. 69, No. 15-16, pp. 2640-2648, ISSN 0266-3538 Chung Deborah, D.L. (2003). *Composite Materials: Functional Materials for Modern Technology,*

Iannucci, L. & Ankersen, J. (2006). An energy based damage model for thin laminated composites. *Composites Science and Technology*, Vol. 66, No. 7-8, pp. 934-951. ISSN 0266-

Jain, J.R. & Ghosh, S. (2009). Damage Evolution in Composites with a Homogenizationbased Continuum Damage Mechanics Model. *International Journal of Damage Mechanics*,

Kormaníková, E., Žmindák, M., Riecky, D. (2011). Strength of composites with fibres. In: *Computational Modelling and Advanced Simulations*, Murín, J., et al., pp. 167-184, ISBN

Laš, V. & Zemčík, R. (2008). Progressive Damage of Unidirectional Composite Panels.

Mishnaevsky, L., Jr (2007). Computational mesomechanics of Composites, Numerical analysis of the effect of microstructures of composites on their strength and damage

Neto, E.A. de Souza, Peric, D., Owen, D.R.J. (2008). *Computational Methods for Plasticity:* 

Reddy, J.N., Miravete, A. (1995), Practical analysis of composite laminates, ISBN 0-8493-

Riccio, A. & Pietropaoli, E. (2008). Modeling Damage Propagation in Composite Plates with Embedded Delamination under Compressive Load. *Journal of Composite Materials*, Vol.

Sládek, J., Sládek, V., Jakubovičová, L (2002). *Application of Boundary Element Methods in Fracture Mechanics,* ISBN 80-968823-0-9, University of Žilina, Faculty of Mechanical

978-94-007-0316-2 Springer Science + Business Media B.V., Dordercht

*Journal of Composite Materials*, Vol. 42, No. 1, pp. 25-44, ISSN 0021-9983

*Theory and Applications,* ISBN 978-0-470-69452-7, John Wiley & Sons, Ltd.

*Journal of Composite Materials.* Vol. 35, No.15, pp. 1325-1342. ISSN 0021-9983

*of Composite Materials.* Vol. 39, No.4, pp. 323-345. ISSN 0021-9983

No. 1-12, pp. 190-200, ISSN 0734-743X

2383-2391, ISSN 0266-3538

ISBN 1-85233-665-X, Springer.

Vol. 18, No. 6, pp. 533-568, ISSN 1056-7895

42, No. 13, pp. 1309-1335, ISSN 0021-9983

1359-835X.

3538

resistance.

9401-5, CRC Press

Engineering, Žilina

Damage model presented in this chapter has been utilized in this solution. This damage model is suitable for elastic damage of fiber-reinforced composite materials with polymer matrix. The postulated damage surface reduces to the Tsai-Wu surface in stress space. The problem of elastic damage is considered as material nonlinearity, so we get system of nonlinear equations. The Newton-Raphson method has been used for solving this system of nonlinear equations. Evolution of damage has been solved using the return-mapping algorithm. Flowchart of this algorithm was also presented. Numerical example of one problem for two different materials was presented next. Own program created in MATLAB language was used for this analysis. Four-node layered plate finite elements based on Kirchhoff (classical) plate theory were used. The analysis results show that change of material as well as the presence and values of shear stress have significant influence on the evolution of damage as well as on location of critical damage and load at which the critical level of damage will be reached. Critical damage level has not necessary to be reached in places with maximum magnitude of equivalent stress, but can be reached in other places.

## **Author details**

Milan Žmindák and Martin Dudinský *University of Žilina, Slovakia* 

## **Acknowledgement**

The authors gratefully acknowledge the support by the Slovak Grant Agency VEGA 1/1226/12 and Slovak Science and Technology Assistance Agency registered under number APVV-0169-07.

## **5. References**

Ansys v.11*, Theory manual* (2007), ANSYS, Inc., Southpointe, PA.

Azouaoui, K. et al. 2010. Evaluation of impact fatigue damage in glass/epoxy composite laminate. *International Journal of Fatigue*, Vol. 32, No. 2, pp. 443-452, ISSN 0142-1123

Bathe, K. J. (1996). *Finite Element Procedures*, Prentice Hall, ISBN 0-13-301458-4, New Jersey


Clegg, R. A. et al. (2006). Hypervelocity impact damage prediction in composites: Part Imaterial model and characterisation. *International Journal of Impact Engineering.* Vol. 33, No. 1-12, pp. 190-200, ISSN 0734-743X

150 Composites and Their Properties

**Author details** 

*University of Žilina, Slovakia* 

5433-3, Boca Raton

73-93, ISSN 1056-7895

**Acknowledgement** 

APVV-0169-07.

**5. References** 

3060.

Milan Žmindák and Martin Dudinský

needed to obtain a reasonable approximation

of delamination front. The solution converges quickly because a small number of plates is

Damage model presented in this chapter has been utilized in this solution. This damage model is suitable for elastic damage of fiber-reinforced composite materials with polymer matrix. The postulated damage surface reduces to the Tsai-Wu surface in stress space. The problem of elastic damage is considered as material nonlinearity, so we get system of nonlinear equations. The Newton-Raphson method has been used for solving this system of nonlinear equations. Evolution of damage has been solved using the return-mapping algorithm. Flowchart of this algorithm was also presented. Numerical example of one problem for two different materials was presented next. Own program created in MATLAB language was used for this analysis. Four-node layered plate finite elements based on Kirchhoff (classical) plate theory were used. The analysis results show that change of material as well as the presence and values of shear stress have significant influence on the evolution of damage as well as on location of critical damage and load at which the critical level of damage will be reached. Critical damage level has not necessary to be reached in places with maximum magnitude of equivalent stress, but can be reached in other places.

The authors gratefully acknowledge the support by the Slovak Grant Agency VEGA 1/1226/12 and Slovak Science and Technology Assistance Agency registered under number

Azouaoui, K. et al. 2010. Evaluation of impact fatigue damage in glass/epoxy composite laminate. *International Journal of Fatigue*, Vol. 32, No. 2, pp. 443-452, ISSN 0142-1123 Bathe, K. J. (1996). *Finite Element Procedures*, Prentice Hall, ISBN 0-13-301458-4, New Jersey Barbero, E. J. (2007). *Finite Element Analysis of Composite Materials*, CRC Press, ISBN 1-4200-

Barbero, E. J. & De Vivo, L. (2001). A Constitutive Model for Elastic Damage in Fiber-Reinforced PMC Laminae*. International Journal of Damage Mechanics.* Vol. 10, No. 1, pp.

Carrera, E. (2002). Theories and Finite Elements for Multilayered, Anisotropic, Composite Plates and Shells. *Arch. Comput. Meth. Engng*. Vol. 9, Nr. 2, 2002, pp. 87-140, ISSN 1134-

Ansys v.11*, Theory manual* (2007), ANSYS, Inc., Southpointe, PA.


Shu, X. (2006). A generalised model of laminated composite plates with interfacial damage. *Composite Structures*, Vol. 74, No. 2, pp. 237-246, ISSN 0263-8223

**Chapter 8** 

© 2012 Wang and Zhang, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

© 2012 Wang and Zhang, licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

**Numerical Modelling of Damage Evolution** 

Demands on high performance materials for use in a spectrum of structural and nonstructural applications are increasing to the point that monolithic materials cannot fully satisfy multifunctional requirements. One approach that has emerged is to develop advanced fiber composites whose properties are tailored for wide applications in engineering, such as turbine generators and ultrasonic aircrafts (Miracle, 2005). Composite materials are composed of reinforcements and matrices, in which mechanical behaviors can be seriously affected by their micro-structures. Therefore, it becomes necessary to deal with these materials from a micromechanical view point. Most importantly, an urgent task in the field of fiber-reinforced composites is to develop a relationship between the material

Damage and failure in fiber-reinforced composites evolves at several different length scales. At the smallest scale, pre-existing defects in fibers grow. Due to statistical distribution of such defect, the fiber strength exhibits large variability. Therefore, as increasing load is applied to the composite, the weakest fiber will first break. The loads carried by the broken fiber are redistributed among the remaining unbroken fibers and matrix as determined by the constitutive response of the fibers, matrix and interface. And this then causes other fibers near the failure site to fail and thus shed further load to intact fibers. Consequently, the failure of fibrous composites goes through a very complicated damage evolution, which is a combination of fiber fracture, matrix deformation and interfacial debonding and slipping around the fiber breaks, before it reaches ultimate failure. It is obvious that the connection between the microstructural scale and the macroscopic scale is nontrivial and involves

structure and its performance under more severe loads and environments.

mechanics, stochastics, and volume scaling (Curtin, 1999).

**and Failure Behavior of Continuous** 

**Fiber Reinforced Composites** 

Additional information is available at the end of the chapter

F. Wang and J. Q. Zhang

http://dx.doi.org/10.5772/48481

**1. Introduction** 

