Open access

Finite Element Implementation of Failure and Damage Simulation in Composite Plates

Written By

Milan Žmindák and Martin Dudinský

Submitted: 16 November 2011 Published: 22 August 2012

DOI: 10.5772/48248

From the Edited Volume

Composites and Their Properties

Edited by Ning Hu

Chapter metrics overview

3,849 Chapter Downloads

View Full Metrics

1. Introduction

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.

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)

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 micro-defects 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.

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 fiber-matrix 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.

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 continuum-based 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 ΩDwith a smooth frontΩD. 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


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, ui0,vi0,wi0, are the displacements at the mid-surface of the i-th plate element and ψxi(x,y),ψyi(x,y) denote rotations of transverse normals about y and x, respectively.

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

{εxxiεyyiγxyi}={ui0xvi0yui0y+vi0x}, {κxxiκyyiκxyi}={ψxi0xψyi0yψxi0y+ψyi0x}, {γyziγxzi}={ψxi0+wi0yψxi0+wi0x}

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ΩΩD, 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ΩΩD 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


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

In order to avoid interpenetration between delaminated sublaminates in the delaminated regionΩD, 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


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 ΩΩDand delaminated ΩD areas. Consequently constitutive law can be expressed as


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


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

The critical ERR (GIc,GIIc,GIIIc) material properties can be evaluated from experimental procedures. The closed-form expressions for the ERR are (Barbero, 2008)

G=GI(s)+GII(s)+GIII(s)GI(s)={limkz,kxy 12kzΔw2(s)    if    Δw(s)00                 if          Δw(s)<0 }GII(s)=limkz,kxy    12kxyΔun2(s)GIII(s)= limkz,kxy   12kxyΔut2(s)E8

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 ΩΩDup 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

GI(A)=(12RAzΔwBB´ΔnΔt) , GII(A)=(12RAnΔunBB´ΔnΔt) , GIII(A)=(12RAtΔutBB´ΔnΔt)E9

where RAz is the reaction in the spring element connecting node A in the z-direction, ΔwBB 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

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Δt. In (9), the same element size is assumed for elements ahead of and behind the delamination front. Value of Δt/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ΩΩD.

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).

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, vxy = vxz = vyz= 0.3. The lower sublamina has these properties: Ex= 70000 MPa, Ey = Ez = 21000 MPa, Gyz = 2100 MPa, Gxy = Gxz = 2333 MPa, vxy = vxz = vyz= 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 GIfor delamination type I

Model II

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, GIIin the x direction and GIIIin 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 ofGIII. 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, GIIin 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.

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 ofGII for model III


3. Continuum damage mechanics

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):

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

  • fracture mechanics approach (based on energy release rates),

  • plasticity or yield surface approach,

  • damage mechanics approach

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:


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


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


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).

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)


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

Y1=1Ω12(S¯11Ω14σ12+S¯12Ω12Ω22σ1σ2+S¯66Ω12Ω22σ62)Y2=1Ω22(S¯22Ω24σ22+S¯12Ω12Ω22σ1σ2+S¯66Ω12Ω22σ62) Y3=0E15

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

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


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


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 CED in lamina coordinate system, which reflects the effect of damage on the behavior of structure. Flowchart of this algorithm is described in Fig. 10.

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


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.

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

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

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.


  1. 1. Ansys v.11, Theory manual2007ANSYS, Inc., Southpointe, PA.
  2. 2. al.2010Evaluation of impact fatigue damage in glass/epoxy composite laminate. International Journal of Fatigue, 3224434520142-1123
  3. 3. BatheK. J.1996Finite Element Procedures, Prentice Hall, 0-13301-458-4Jersey
  4. 4. BarberoE. J.2007Finite Element Analysis of Composite Materials, CRC Press, 1-42005-433-3Raton
  5. 5. BarberoE. J.De VivoL.2001A Constitutive Model for Elastic Damage in Fiber-Reinforced PMC Laminae. International Journal of Damage Mechanics. 10173931056-7895
  6. 6. CarreraE.2002Theories and Finite Elements for Multilayered, Anisotropic, Composite Plates and Shells. Arch. Comput. Meth. Engng. 9Nr. 2, 2002, 871401134-3060
  7. 7. CleggR. al.2006Hypervelocity impact damage prediction in composites: Part I-material model and characterisation. International Journal of Impact Engineering. 331-121902000073-4743X
  8. 8. CuiW.WisnomM. A.1993Combined stress-based and fracture mechanics-based model for predicting delamination in composites. Composites 1993; 244674740135-9835X.
  9. 9. DávilaC. G.CamanhoP. P.RoseCh. A.2005Failure criteria for FRP Laminates. Journal of Composite Materials. 3943233450021-9983
  10. 10. ElmarakbiA.HuN.FukunagaH.2009Finite element simulation of delamination growth in composite materials using LS-DYNA. Composites Science & Technology, 69 238323910266-3538
  11. 11. FenskeM. T.VizziniA. J.2001The inclusion of in-plane stresses in delamination Criteria. Journal of Composite Materials. 3515132513420021-9983
  12. 12. al.2010Effect of matrix cracking and material uncertainty on composite plates. Reliability Engineering & System Safety. 9577167280951-8320
  13. 13. al.2009Element-free Galerkin modelling of composite damage. Composites Science and Technology. 6915-16264026480266-3538
  14. 14. ChungDeborah. D. L.2003Composite Materials: Functional Materials for Modern Technology, 185233665Springer.
  15. 15. IannucciL.AnkersenJ.2006An energy based damage model for thin laminated composites. Composites Science and Technology, 667-89349510266-3538
  16. 16. JainJ. R.GhoshS.2009Damage Evolution in Composites with a Homogenization-based Continuum Damage Mechanics Model. International Journal of Damage Mechanics, 1865335681056-7895
  17. 17. KormaníkováE.ŽmindákM.RieckyD.2011Strength of composites with fibres. In: Computational Modelling and Advanced Simulations, Murín, J., et al., 167184978-9-40070-316-2Springer Science + Business Media B.V., Dordercht
  18. 18. LašV.ZemčíkR.2008Progressive Damage of Unidirectional Composite Panels. Journal of Composite Materials, 42125440021-9983
  19. 19. MishnaevskyL.Jr 2007Computational mesomechanics of Composites, Numerical analysis of the effect of microstructures of composites on their strength and damage resistance.
  20. 20. NetoE. SouzaPeric. D.OwenD. R. J.2008Computational Methods for Plasticity: Theory and Applications, 978-0-47069-452-7John Wiley & Sons, Ltd.
  21. 21. ReddyJ. N.MiraveteA.1995Practical analysis of composite laminates, 0-84939-401-5Press
  22. 22. RiccioA.PietropaoliE.2008Modeling Damage Propagation in Composite Plates with Embedded Delamination under Compressive Load. Journal of Composite Materials, 4213130913350021-9983
  23. 23. SládekJ.SládekV.JakubovičováL.2002Application of Boundary Element Methods in Fracture Mechanics, 8-09688-230-9of Žilina, Faculty of Mechanical Engineering, Žilina
  24. 24. ShuX.2006A generalised model of laminated composite plates with interfacial damage. Composite Structures, 7422372460263-8223
  25. 25. al.2008Damage prediction in composite plates subjected to low velocity impact. Composite Structures, 83173820263-8223
  26. 26. TuminoD.CampelloF.CatalanottiG.2007A continuum damage model to simulate failure in composite plates under uniaxial compression. Express Polymer Letters 1, 115230178-8618X
  27. 27. ZhangZ.WangS.2009Buckling, post-buckling and delamination propagation in debonded composite laminates. Part 1: Theoretical development, Composite Structures, 881211300263-8223
  28. 28. ŽmindákM.RieckyD.DanišovičS.2010Finite element implementation of failure and damage models for composite structures, Proceedings of the XV international conference “ Machine Modelling and Simulation, Warszaw university of technology, Krasiczyn, August, 2010.
  29. 29. ŽmindákM.RieckyD.SoukupJ.2010Failure of Composites with Short Fibers. Communications. 4/2010, 3339ISSN-13354205.

Written By

Milan Žmindák and Martin Dudinský

Submitted: 16 November 2011 Published: 22 August 2012