Modelling of fracture processes in laminate composite plates with embedded delamination

This work is devoted to research and development of models for the processes of nucleation and growth of delamination in laminated composite materials. Possibilities of fracture mechanics approaches using virtual crack closure technique as well as progressive failure analysis based on the failure criteria are demonstrated. Influence of the initial delamination size on the rate of growth of the defect is studied. A comparative analysis of the applicability of fracture criteria to modeling of samples with delamination were performed, including implementation of element deletion algorithm. The numerical results obtained with the studied methods were compared.


INTRODUCTION
omposite materials are widely used in demanding and high-performance applications, largely due to their relatively high strength and stiffness and, at the same time, low weight and density.However, these materials are affected by development of internal defects caused by the peculiarities of their microstructure.To predict the behavior of composites under operational loads, an important and urgent challenge is to develop numerical approaches for description of the processes of accumulation of damage and subsequent failure of the materials.Modelling of failure processes is challenging, especially for the composite materials, due to the specifics of their microstructure.One of the most common types of composite materials are laminated composites.Their numerical simulation usually considers several scale levels: • micro-scale, where the constituents that composing a representative volume of the material can be distinguished (fiber and matrix); • ply scale, which consists of fibers, stacked in a certain way, and linking matrix.Usually a separate ply has orthotropic or transversal-isotropic properties.• laminate scale, which is formed by superimposing of several plies with different fiber orientations.Depending on the stacking of plies, the laminate may possess anisotropic, orthotropic or quasi-isotropic properties.For simplicity, the effective properties of the laminate describing homogeneous behavior are often used in calculations.

C
In most cases, only prediction of initiation of fracture is not sufficient because the local failure does not necessarily lead to the loss of the bearing capacity of structure due to stress redistribution between the components of the material [1].It is necessary to understand underpinning mechanisms of the process at each scale level to be able to create reliable models of failure of composites.Constituents of composites are described by phenomenological approaches.Their strength characteristics can be described in terms of traditional theories or determined experimentally.During the processes of manufacturing and exploitation of composite materials, their microstructure on the ply scale may be subject to various structural changes, contributing to development of the failure processes.These include matrix cracking, fiber breakage or bending, potentially accompanied by delamination between matrix and fibers, as well as between the plies [1][2][3][4].Delamination is one of the most commonly occurring types of defects in composites, it is defined as separation of plies with forming of voids and growing of cracks, oriented by normal to the ply [5].It may be caused by local geometric nonuniformities that appeared in the manufacturing process and results in out-of-plane loads.Also, debonding can be caused by differences between the properties of the constituents of the composite, which leads to emergence of interlaminar stresses.Cracking of the matrix, which occurred during the deformation process, gives the same effect.Another cause of interlaminar stresses and, consequently, delamination, may be the influence of temperature and humidity.Finally, the different types of operating loads, such as prolonged static and cyclic loading as well as impact effects, can directly or indirectly cause the appearance of defects in the form of deboning and delamination [6][7][8].Whatever the original cause of the formation of the delamination is, its development can lead to premature decrease of structures' resource, and in some cases, even cause catastrophic loss of load-bearing capacity.Understanding influence of the defects on the properties of composite structures is necessary for the prediction of their durability.In this regard, urgent problem is development of models describing the delamination process as well as development of the initiation of damage in composites.This work is devoted to investigation of numerical models for processes of delamination nucleation and development in laminate composite materials, numerical modeling of these processes depending on the initial parameters of the defect, as well as studying of the influence of the defect on the strength characteristics of the materials.Currently, many approaches had been developed for creation of mathematical models of the mechanical behavior of composites considering the start of fracture process and its further development.There are two main groups: approaches based on the theory of strength, and approaches based on fracture mechanics [9].Approaches related to the first group use strength criteria to determine the moment of the first ply failure in the laminate.Such criteria may include various parameters, such as the value of the stress tensor components, material characteristics, strain, force, displacement and others.The foundations for the development of strength criteria for the composites were laid in the works of Hill [10], Tsai [11] and others, who formulated the criteria on the basis of the relations between the components of the stress tensor.More sophisticated criteria based on a set of different possible failure modes, including separating of the matrix and fibers failure, were proposed by Hashin [12,13], Puck [14], Chang [15] and were developed in the works of other authors [9].Theories that are able to describe development of existing damage had been developed within the framework of fracture mechanics.At the same time, usually the classical theories do not say anything about the origin of the defect.One of the most widely used models are based on monitoring of the rate of deformation energy released during the propagation of the defect, and its comparison with a threshold value of the strain energy release rate for the particular material [9].To determine the energy release rate, a number of techniques, implemented with the finite element methods (FEM), have been developed [16][17][18].Over the past decades many models and analytical tools for modeling of debonding processes were created using modifications and combinations of the approaches outlined above.Detailed reviews of works in this area are given, for example, in [1] and [6].In this paper a practical analysis of models of delamination growth will be carried out, involving the virtual crack closure technique as well as the progressive failure analysis using a range of failure criteria.

Virtual crack closure technique (VCCT)
ccording to the concept of this approach, the energy released during the crack expansion onto a certain distance is equal to the energy required to close the crack at the same distance [17,19,20].The energy expended to the crack closure is calculated from the coupling of values of the strength at the crack tip and displacement during movement of the crack tip: where I X and I Z are shear and opening forces at the node I , u  and w  are displacements, corresponding to shear and opening at the node L (Fig. 1).The energy release rate can be calculated then as: where ∆ is formed crack surface.This hypothesis can be used if the crack profile does not change significantly during its growth, and when the value of the distance, for which the rate of strain energy release is being calculated, is relatively small compared with the overall dimensions of the crack.
In the real materials, the crack growth during delamination occurs simultaneously in all three modes of deformation (opening, in-plane shear and out-of-plane shear).In order to take account of this fact, the energy release rate is calculated for each mode ( ) , then they are summed up: and are compared with the critical value C G .The opening of the two nodes and the crack growth occurs with satisfaction of condition: The critical value C G depends on all three modes of deformation and is determined by the mixed criterion.The three dimensional Benzeggah and Kenane (B-K) criterion is often used [21]: where values Ic G , IIc G , and exponential parameter  are determined experimentally.Such experiments are in details described in [22,23].Besides, guidelines for experiments determining the critical strain energy release rate are presented in the ASTM standards [24].
In this paper, VCCT method is used to study the influence of the size of the initial delamination to its propagation speed.
The laminated unidirectional carbon-fiber thermoset polymer with plies orientation   7 0 / 90 and fiber volume fraction 60% has been considered [25].The sample dimensions were 50 x 250 mm (Fig. 2a).The defect in the form of circleshaped delamination with diameter d 10, 20, 30  and thickness 0.3 mm was placed between 6 th and 7 th layers.To implement example calculations, the material physical-mechanical properties were taken from [25,26].They are shown in Tab.The finite element model of the laminate consists of two bonded plates, with 6 and 8 plies, respectively.The plates were modelled using shell elements with a single point of integration for each ply.Computations were performed in SIMULIA Abaqus package.The finite element model takes into account the presence of a defect by defining contact conditions between nodes (Fig. 2b).Thus, the contact conditions are defined between all nodes of the two plates except the nodes within a radius of the defect (Fig. 2c).In the top and bottom of the sample nodes are bonded rigidly, and in the central part their slippage is permitted during deformation.Distributed compression load was set at the upper end of the sample with a fixed displacement value u mm 3.5  . Boundary conditions provide fixed lower end of the plate and restrict movement of side faces of the plate along axis 2. Values (3) and ( 5) are being calculated at each step and for each node during analysis to test condition (4).When the criterion is met, nodes are being disconnected and start to move, simulating the defect growth.Fig. 3 represents the sequential defect growth through the state of bonded nodes at the different values of applied displacements load.The dark color corresponds to debonded nodes, which were initially disconnected or for which the condition (4) was satisfied later.Fig. 4 shows the condition of the samples after the loss of stability with superimposed field of displacement along the axis 3 (directed along the normal to the sample), for sample with a defect size of 20 mm.Fig. 5 shows the force-displacement curve for samples with different initial size of defects.For each plot, there is point marking start of delamination growth from the initial defect, as well as the point corresponding to the time when delamination reaches the edge of the sample.At this stage, loss of strength of the finite elements is not modeled, so strength of the whole sample can be regarded as buckling of the plate, which is characterized by a drop-down jump on the charts.The obtained results confirm that the initial size of the delamination has a significant impact on the ability of the sample to resist applied load.So, for delamination with a diameter of 10 mm, its spread to the edges of the sample occurs at a load of 142 kN and defect growth begins at a value of the compression force of 92.2 kN.Buckling occurs simultaneously with delamination achievement of the edge of the sample.With increasing the diameter of the initial defect up to 20mm, almost half less of load is sufficient for the loss of strength of the sample: delamination reaches the edges at 90 kN.The force needed to start the growth of the defect differs not so significantly: 72.2 kN.Buckling of the plate happens with the force of 112kN applied.In the third case, when the initial delamination zone has a diameter of 30mm, the defect spreads to the border almost instantly (obviously, not without the influence of the proximity to the border).The values of the applied force for the initiation of growth and achievement of the borders are 68.1 kN and 70 kN, respectively.However, the sample continues to resist buckling till 110kN load, which corresponds to the sample with the defect diameter of 20mm.It can be concluded that in some cases the sample continues to maintain the bearing capacity even at the critical growth of delamination.It is also worth noting that the force corresponding to the beginning of the defect, are almost matched for the samples with defect diameter 20mm and 30mm.The next part of the work is devoted to studying the effect of the delamination presence on the strength properties of the individual plies of the laminate.

Progressive failure analysis
This part of the work is devoted to studying the influence of the defect on the models of progressive failure of composites.Progressive failure analysis allows to evaluate criteria of the composite strength at each step of deformation that enables to modify properties of the elements during the analysis, depending on whether the failure criterion for them is fulfilled or not.The studied sample in this case was modelled using three-dimensional solid finite elements.In addition, for a more complete account of the defect, each layer was modeled as a separate set of finite elements, with the orientation set corresponding to the ply orientation.A defect with a diameter of 20mm was placed in the center of the sample by adding embedded ply (film) in the selected group of finite elements (Fig. 6).The properties of this layer, simulating a delamination, were taken to be the properties of the laminate, multiplied by 10 -6 .Boundary conditions were specified, as discussed in the previous case, in the form of the upper end displacements of 3.5 mm.The properties of the material correspond to the previously studied model and are presented in Tab. 1.
The following criteria of failure of the finite elements were used criteria: Maximal Stress, Tsai-Hill [27], MCT [28] and Hashin [13].The progressive failure was implemented in Autodesk Helius PFA.Below is a brief description of each used criterion.
In the maximal stress criteria, values of the stress tensor components are compared with the effective homogenized longitudinal and transverse tensile, compressive and shear loading strength properties:  Another criterion, that takes into account the same deformation modes, is the Tsai-Hill criterion [27], according to which failure begins when the following conditions are met: Two of the above criteria operate effective characteristics of the laminate, and does not imply a separate analysis of the behavior of the internal constituents of the composite.The criteria presented below, allow to separate deformation modes of matrix and fibers individually.The required input data, such as the properties of internal constituents of the materials, are presented in Tab. 1.After fulfillment of the fiber breaking conditions in an element, its stiffness constants are reduced to 1% of their initial value.After the matrix failure, the values of an element stiffness constants are reduced to 10% of their original.The Hashin criterion distinguishes four modes of failure: failure of fiber and matrix in tension or compression [13].The fiber breaks when tension or compression where cr 11   and cr 11   is the laminate longitudinal tension and compression strength (see Tab. 2).
The matrix failure occurs when the following inequalities are satisfied:    are transverse and longitudinal shear strength constants.The multicontinuum theory (MCT) criterion of composite materials, proposed in [28], also allows to split degradation of matrix and inclusions.It is assumed that fiber and matrix have transversely isotropic properties.Within multicontinuum theory for composite materials, the following relationships between stress and strain fields in the material with effective properties and the individual constituents of the composite are introduced [28]: where С σ and С ε are stress and strain fields for homogenized material; f σ and f ε are stress and strain in fibers; m σ , m ε are stress and strain in matrix; f p and m p is volume fraction of fiber and matrix; structural elasticity moduli of homogenized composites, fibers and matrix, respectively.Stresses and strains fields of in the matrix m and the fiber f can be expressed as follows: where Thus, the strain field in inclusions and matrix may be expressed from the effective strains using relations ( 7), ( 9) and ( 10).This representation is used for finite element modeling for verification of fulfillment of failure criteria of composites.Fiber failure criterion is set with the following equation: where Matrix failure occurs if the following identities are true: where Fig. 7 compares von Mises stress tensor fields for the plies of different orientations, closest to the zone with a defect at the time just before reaching the breaking load.The difference is caused by the fact that during the deformation mechanical properties of the parts of the elements have already been reduced (according to the above stiffness reduction schemes) due to implementation of criteria.For Hashin and MCT methods, additional algorithms for removing elements have been implemented.Thus, elements were subject to removal if the conditions of fiber failure criterion were met.Fig. 8 shows a moment of fracture for the modelled sample using the four selected approaches.The figures presented for the different orientation of the plies, closest to the defective area.Red color allocates elements for which the failure criteria were met and, therefore, the strength properties have been decreased.For the Hashin criteria and MCT cases, red color designates elements for which the conditions of matrix failure were implemented, since the elements with fractured fiber have been removed, as was presumed by the algorithm.

RESULTS AND DISCUSSION
he results shown in Fig. 7 and 8 indicate that according to all models, except the Tsai-Hill model, failure of the plies occurs in the defect area.Tsai-Hill model (Fig. 9b) had appeared to be not sensitive to the defect area and shows damage initiation in the zone closest to the end with applied load.Since the Hashin and MCT criteria distinguish modes of fiber and matrix deformation, models based on them most reliably describe failure of plies of different orientation.Thus, Fig. 8c and 8d demonstrate that for the transversely oriented plies the number of failed finite elements is much bigger as compared to longitudinal oriented plies.Criterion of maximal stress (Fig. 8a) is implemented in the area around the defect, but does not give an opportunity to explore the influence of orthotropic properties of the layers.Fig. 9 shows a comparison of the deformation diagrams for the progressive failure analysis models based on the criteria and models of fracture mechanics using the virtual crack closure technique for the sample with an initial diameter of delamination 20 mm.The diagram of the sample without defect, modelled using MCT criterion, is also presented for comparison.Based on the results, the difference in the load, required to break the sample with a defect, for a set of methods used is about 7kN.The presence of the defect reduces the load down to 20 kN.Model based on the criterion of maximal stress produces results for breaking load comparable to a model based on the Hashin criterion, since both criteria use similar instruments for strength assessment of the laminate (in the case of the maximal stress criterion) or fiber (in the case of the Hashin criterion), consisting in direct comparison of stress with tensile strength.According to the Tsai-Hill criterion, the sample will collapse at a load of 120kN, which indicates that although this model fails to localize the defect, its presence has indirect impact.The energy approach using virtual crack closure technique simulates the lowest breaking load (about 118 kN) even though, as mentioned above, in this case buckling of the sample is modeled without loss of elements' stiffness.
T CONCLUSIONS everal methods of simulation of the delamination in composite materials samples have been investigated.The effect of defect size on the parameters of deformation and failure of the sample was established.Computational experiments were performed to explore the possibility of the delamination simulation using progressive failure models, applying various criteria.It can be stated that a defect in the form of delamination has a significant influence on the strength of the samples and its presence can be modeled using both methods of fracture mechanics and methods based on failure criteria.The energy approaches for modeling the growth of the delamination seem more promising, however, because the methods using the values of stress fields in the failure criteria, in general, impose additional conditions on the model's geometrical properties.Thus, the appearance of a singularity due to distortion of the geometry or features of the finite-element mesh, can prevent achieving reliable results.At the same time, disadvantage of the energy techniques is that they require a predetermine position of the delamination in the structure.Thus, currently the combined methods are being developed, which utilize failure criteria for locating the start of delamination and the fracture mechanics methods to modeling its subsequent growth [29].The initial information about the presence of the delamination can be also obtained by experimental approaches.Among them are the non-destructive testing methods such as acoustic emission, ultrasonic examination and others [30][31][32][33].Another promising direction is embedding of the sensor elements, such as piezoelectric sensors and sensors, based on optical fibers [34][35][36].Indications of such sensors may be also used for verification of failure models as well as predictions of durability of structures based on them.Another approach to improve the models of delamination in composites is account of a different nature of defects.Thus, in the case of manufacturing-induced defects, the occurrence of delamination is often accompanied by the presence of socalled "resin pockets".In this case, a more accurate model is required, that considers resin as a separate component with isotropic properties.Also, the temperature factors can affect the rate of growth, the effect of which may be incorporated in the models with the appropriate constitutive relations.Thus, based on the basic principles, including those reflected in this paper, for each specific case of the composite structure the model that most closely simulates the formation and development of the defects can be created.

Figure 2 :
Figure 2: Geometrical properties of the sample.

Figure 4 :
Figure 4: State of sample with defect of 20 mm in diameter after buckling with superimposed displacements along axis 3.

22   or cr 22  23  and cr 12
 is the laminate transversal tension and compression strength, cr

, GPa 22 E , GPa 12 ν 12 G , GPa 13 G , GPa 23 G , GPa
1. Critical strain energy release rate values for the material are IC 11 E

Table 1 :
Properties of the ply and constituents.
The strength characteristics for the sample are shown in Tab. 2. Failure begins when corresponding stress tensor components in any finite element exceed the effective strength properties.Triggering criterion indicates simultaneously fracture of fibers and matrix, while the finite element effective strength properties are being decreased by 10%.