## 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 *G*_{r} 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 *h*_{0}. 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 *h*_{1}, *h*_{2} and each sublaminate consist the upper *n*_{u} plates and the lower *n*_{l} 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)

where *u*_{i}, *v*_{i} refer to the in-plane displacements, and *w*_{i} to the transverse displacements through the thickness of the *i-*th plate element, *i-*th plate element and *y* and *x*, respectively.

At the reference surfaces, the membrane strain vector

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*k*_{z},*k*_{xy}, 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 *σ*_{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*Δw≥0*) and a positive stiffness for closing relative displacements (*Δw≤0*), then the contact stress

where *k*_{z} is the penalty number imposing contact constraint and sign is the signum function. A very large value for *k*_{z} 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

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

(8) |

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

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

where *B* and *B'*. These are located immediately ahead of the delamination front along its normal direction passing through *A*.

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

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

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: *E*_{x}= 35000 MPa, *E*_{y} = *E*_{z} = 10500 MPa, *G*_{yz} = 10500 MPa, *G*_{xy} = *G*_{xz} = 1167 MPa, *v*_{xy} = *v*_{xz} = *v*_{yz}= 0.3. The lower sublamina has these properties: *E*_{x}*=* 70000 MPa, *E*_{y}
*= E*_{z} = 21000 MPa, *G*_{yz} = 2100 MPa, *G*_{xy}
*= G*_{xz} = 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, *h*_{1}= 0.5 mm, *h*_{2}= 1mm.

The upper sublaminate is composed by four plates *n*_{u}= 4 and the lower by two *n*_{l} = 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 *k*_{z} a *k*_{xy} = 10^{8} N/mm^{3}. 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.

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, _{I}. The values of G_{II} 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, _{II} in the range of ( 0, 14).10^{-3} and the value of G_{III} 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.

## 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

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.

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 Y_{1}, Y_{2} and Y_{3} can be calculated by means of relations

(15) |

where stresses and components of matrix

#### 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}, c_{1} and c_{2} 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 F_{6} is the in-plane shear strength, G_{12} 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 D_{1}, D_{2} 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 - D_{1}) (1 - D_{2}) with k_{s} value for given material. If the value of this product is less or equal to k_{s} 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 D_{1}, D_{2}, and δ in integration point from the start of last performed iteration, _{1}, D_{2}, 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.

### 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 J_{33} and H_{3} 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.

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 |

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.

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.

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