Computational multi-scale analysis of simultaneous processes of delamination and damage accumulation in laminated composites

When studying mechanical behavior of structures made of laminated polymer composite materials, issues related to their relatively high susceptibility to damage, such as delamination, fiber breakage and matrix failure, play significant role. Key aspects are also related to accurate description of the internal structures and heterogeneity of the material. This work is aimed at a numerical study of the processes of damage accumulation in plies as well as development of delamination in laminated polymer composite L-shaped specimen taking into account microstructural parameters. The influence of contact conditions between the plies on the model’s behavior has been studied. It was established that the delamination processes inside the specimen can initiate at the early stages of deformation, which is reflected in the difference of the results obtained using models with various numbers of separated plies. Comparison of the results of modelled elastic behavior of the specimen and models with the developed delamination leads to conclusions about influence of a number of model parameters on the results obtained before and after propagation of delamination. In particular, such factors as the geometry and mechanical properties of models, the number and conditions for specifying contact pairs in the model were analyzed.


INTRODUCTION
he difficulties in modeling of failure of laminated composite structures are related to the inhomogeneities of microstructure of composite materials. A characteristic feature of laminated polymer composites is that the origin of the fracture processes in the material does not always lead to complete failure of the structure. Thus, until the structural integrity is lost, damage can be accumulated at various scale levels. T When analyzing the mechanical behavior of polymer composite structures, an important task is to create numerical models that allow to reproduce the deformation behavior of the structure not only in the elastic zone, but also after the onset of the development of internal defects and subsequent failure process. To obtain more accurate numerical results it is necessary to take into account both the structural features of the laminated composites (distribution of epoxy between plies, presence of epoxy "pockets" and internal defects) and a number of mechanisms that cause the occurrence and progression of the defect, implemented in the form of appropriate criteria and fracture models. In addition, interaction of microstructural elements plays a key role in the propagation of fracture processes in composites. Accordingly, besides the development of defects arising at the scale level of the composite ply, it is necessary to model the effect of damage accumulation at the micro-scale level. The most common types of fracture in laminated composites are fiber rupture, matrix cracking and delamination [1][2][3]. The traditional approach when creating failure models of composites is to use the failure criteria on a scale of a ply. These criteria compare the stress state in a ply with the strength characteristics or allowable deformation of the material and are defined in the local coordinate system. At present, a large number of fracture criteria have been developed for composite materials that take into account these mechanisms to some extent and which depend on a different sets of critical constants [4][5][6][7][8]. The foundations for the development of strength criteria for the composites were laid in the works of Hill [9], Tsai [10] 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, separately, matrix and fibers failure, were initially proposed, by Hashin [11,12], Puck [13] and Chang [14]. In order to take into account the structural defects, the approaches of continuum damage mechanics (CDM) are used, according to which after the failure criterion in a ply is met, the material with defects is replaced by a fictious continuous material without defects, but with lowered elastic properties [15]. To simulate the growth of delamination, approaches of the fracture mechanics are used. Since laminated composites demonstrate brittle behavior, linear elastic fracture mechanics (LEFM) are used to model delamination, one of the common approaches of which is the virtual crack closure technique (VCCT) ) [16][17][18][19][20]. The essence of the method in application to the finite element approach lies in the separation of the nodes of the finite element mesh located on the initially connected surfaces of the plies when specified fracture criterion is fulfilled. For the same model of composite structure, the use of different fracture theories can lead to significant differences in the results of evaluation of its mechanical response. When creating models of structures made of laminated composites in order to take into account all the basic mechanisms of failure it is necessary to (i) to create a laminated geometric ply-byply model of the composite material; (ii) analyze the state of contact between the plies on the basis of approaches of fracture mechanics that allow to numerically implement the process of separation of plies; (iii) evaluate the strength characteristics and damage accumulation in the ply not only by using its effective properties, but also by establishing its relation to the analysis of local fields in microstructural components. This work is aimed at studying the influence of damage accumulation processes and changes in the stiffness properties in the woven ply of a polymer composite material on the development of delamination in laminated polymer composites structures. An L-shaped specimen subject to load leading to delamination and, at the same time, accumulation of defects and damage in the plies, is investigated. The method of estimation of microstructural parameters as well as physical and mechanical properties of the composite material is based on a multi-scale mathematical model taking into account the geometry of the material, anisotropy of physical and mechanical properties of the plies and mechanical properties of the matrix. The virtual crack closure technique (VCCT) with the application of the Benzeggagh-Kenane (BK) criterion is used for modelling the propagating delamination. The assessment of plies damage is based on two criteria: a multicomponent criterion, which uses the simplest fracture indicators for various components of the stress tensor, and more complex Hashin criterion [12], taking into account the analysis of local microstructural fields. The numerical results of modeling of the mechanical behavior of the specimen using different criteria are compared and the relationship between the processes of delamination and damage accumulation in plies is studied.

MULTI-SCALE MODEL OF FRACTURE IN LAMINATED COMPOSITES
Virtual crack closure technique ccording to the concept of this technique, the energy released by the expansion of the crack at a certain distance is equal to the energy required to close the crack at the same distance. In the framework of the finite element simulation, this energy can be calculated from the nodal forces at the delamination front and the corresponding A nodal displacements behind the delamination front [18]: where I X and I Z are shear and opening forces in the node I , u  and w  are displacements corresponding to the shear and opening in the node L ( Fig. 1). The energy release rate is then calculated as: where ∆ is a crack surface. In real materials, the debonding crack usually grows simultaneously in all three modes of deformation (opening, in-plane shear and out-of-plane shear). In order to take this into account, the energy release rates are calculated for each mode ( , , ) I II III G G G and the subsequently summed: The latter is compared with the critical value C G . In this case, the beginning of the opening of two nodes and growth of the crack occurs when the following condition is met: The critical value c G depends on all three modes of deformation and is defined using the mixed criterion. One of the most used criterion in three-dimension space is the Benzeggah and Kenane (BK) criterion [21]:   1   I  II  III   II  III  IC  IIC  IC  I  I I  I I where the constants Ic G , IIc G are determined experimentally for each laminated composite material.

Progressive failure model
The progressive failure model allows to evaluate the performance of the strength criteria in every ply of composite material at each step of deformation, which makes it possible to change the properties of finite elements during the analysis, depending on whether the criterion is fulfilled for them or not [22][23][24][25][26][27][28][29][30][31].
Considering the material with initially linear elastic behavior, the ratio of stresses of the initial and damaged material is expressed, respectively, as follows: where Ĉ denotes an effective (undamaged) stiffness tensor,   C D is the stiffness tensor depending on the damage tensor D. The relationship between the effective stress  and conditional stress  is written using the damage tensor D: Using the principle of compatibility (equivalence) of strains, the tensor   To calculate the damage variables ij D , it is necessary to specify the criterion of failure and the damage evolution law. In this work, the damage accumulation in plies is modeled using the multicomponent fracture criterion and the Hashin criterion.
The multicomponent fracture criterion includes five fracture indicators corresponding to mechanical response of the composite ply when the specimen is loaded in different directions: 11 11 , 0; where t X is tensile strength in direction 1; c X is compressive strength in direction 1; t Y is tensile strength in direction 2; c Y is compressive strength in direction 2; S is shear strength in the plane (1,2). Hereinafter, it is assumed that the plane 1-2 coincides with the surface of the woven ply.
where 1  changes according to the power damage evolution law: Here 1  , according to [32].

Homogenization procedure
Since failure of a composite structure is a direct result of the processes occurring inside the material, the accuracy of the numerical methods can be improved by analyzing behavior of a ply at the microstructural level, including studying of stress and strain fields in the individual microstructural components. Modern approaches of micromechanics, in addition to calculating the effective properties of a composite material based on microstructural data, allow to simulate the mechanical behavior of local fields. The goal of the homogenization procedure is to give an approximate estimate of the average values of stress and strain fields, both at the macro level and in each phase. In the present work the mean field method of homogenization is employed to study the individual components' behavior. This analytical method is based on the assumption of interrelationship of the average stresses and strains in each phase of the representative volume of the material [9]. The sum of the phases volume fractions of two-phase composites is 1 where M is for matrix, F is for inclusions (fiber in our case). Then, according to the mean field method, the strain concentration tensors can be determined from the following equations: Average volumetric values for strain field F  in the inclusions are associated with average volumetric values for strain field in the matrix M  via the concentration tensor  B , and with macro-strains  via the second concentration tensor  A .
For woven composite materials, the mean field method can be implemented on the basis of a modified Mori-Tanaka homogenization approach [33]. The initial method was proposed by Mori and Tanaka in 1973 and is based on the approximate use of the Eshelby's solution. The application of the Mori-Tanaka approach is obvious when the reinforcing particles can be effectively approximated as ellipsoids. For the textile composites with an organized structure consisting of yarns, the homogenization procedure is performed in the several stages [34]. First, a mechanically equivalent material's microstructure is created, with a simplified geometry for inclusions, and a segmentation of the geometry of the textile ply is carried out. Then, properties of each textile segment are calculated using homogenization formulas for a unidirectional fiber array using the local volume fraction of fibers in the segment, fiber properties and elastic properties of the matrix. The result is a stiffness matrix, expressed in the local coordinate system. For the Mori-Tanaka homogenization, the spatial arrangement of the inclusions does not make difference, the only important geometric factors are orientation and size of the inclusions. After the equivalent set of inclusions has been created, the effective stiffness tensor is obtained as follows: the Eshelby tensors i S are calculated for the inclusions in local coordinates, and the result is converted for the global coordinate system. Then the strain concentration tensors are calculated for all the inclusions: where j p is relative volume fraction the inclusions phase j , M p is relative volume fraction of the matrix, 1

M i
A is calculated by the formula: where i S is the Eshelby tensor for inclusions phase i , i C is inclusions' stiffness tensors, M C is stiffness tensor of the polymer matrix. Further, the effective stiffness matrix of the composite is expressed as follows:

NUMERICAL EXAMPLE
he key parameters for modeling of the fracture processes in the laminated composites are the critical constants included in the fracture criteria, the contact properties between the plies, the number of contact pairs in which delamination can emerge. The algorithm of delamination analysis of laminated composites consists of three main phases -analysis of the elastic zone, the moment of initiation of delamination and the subsequent development of delamination. The laminated L-shaped fiberglass specimen consisting of 12 plies was chosen for the study (Fig. 2). The unit cell of the material's woven microstructure is shown in Fig. 3. The thickness of the ply was set to 0.26 mm. The boundary conditions simulated the real experiment loading -pre-tension of the fastening bolts and displacements of 8 mm applied to the lower part of the structure while the upper bolts remained fixed (Fig. 2). The main features of the model operating in the elastic zone are its geometric structure, elastic properties of structural elements, and the type of boundary conditions. In the studied model, the possibility of simultaneous growth of the debonding zone between (i) five plies (each second pair) and (ii) each pair of plies of the material is implemented. The initial elastic properties of the structural elements of the composite specimen model are presented in Tab. 1. The values of the critical constants used in the BK criterion (5) are given in Tab. 2. The latter properties were obtained from the in-house set of the experiments. The ply's critical constants to be used in multicomponent (9)

RESULTS AND DISCUSSION
he Figs. 4-7 show the results for the finite element model of the specimen with delamination between all pairs of plies and degradation of the plies' properties predicted by the Hashin criterion (12)- (16). In particular, the state of contact between the surfaces of the second and third, as well as the fifth and sixth plies, at different steps of loading are presented: before delamination (u = 2.65mm), at the beginning of delamination (u = 2.7mm) and in the final stage of delamination (u = 8mm). Red color shows an existing contact. For the same loading steps, the field of the damage tensor variable D11 in one of the plies of the corresponding pair is presented. These values change from 0 (blue) to 1 (red). Figs. 8-11 show the same set of results for the model with multicomponent criterion. The results show that in both cases the propagation of delamination and the accumulation of damage processes are not the same in different plies. The degree of delamination and damage is higher in the plies located in the middle of the specimen. It was found that using different criteria of damage accumulation in the ply, different rates of propagation of delamination between various pairs of plies are observed. However, in some cases, the activation of fracture indicators according to the selected criterion and the weakening of the strength characteristics of the ply is due to the already occurring process of delamination and the associated redistribution of stresses. In particular, this can be observed in Fig. 4 for the third ply (Hashin criterion) and in Fig. 10 for the sixth ply (multicomponent criterion). At the same time, there is a reverse process -the accumulation of damage in the matrix of the plies enforces development of delamination. This effect was observed for the sixth ply in the model with Hashin criterion (Fig. 5), as well as for the third ply in the model with multicomponent criterion (Fig. 8). Thus, the delamination process can begin regardless of the degree of damage accumulation in the ply, and the subsequent increase in damage contributes to the growth of the delamination rate. Figs. 6 and 7 show the stress distribution in the matrix obtained at different loading steps using the reverse application of the Mori-Tanaka homogenization approach. As can be seen from these figures, the maximum stresses in the matrix decrease after initiation of the delamination process of the adjacent pair of plies. From the Fig. 9 it follows that the stress values in the matrix, in addition to the growth of damage in the ply, are also affected by the state of contact of the neighboring plies -a stress concentrator is present in the zone of the delamination initiation. However, after separation of the plies, the magnitude of the stresses in the ply decreases sharply. In the framework of the analysis of the initial parameters of the model, the influence of the number of the specified contact conditions between the plies of the model that allows modeling the development of delamination is studied. Two cases were investigated for this purpose. In the first case (VCCT-1), the contact conditions with the virtual crack closure technique were set only between each second contact pair (5 in total), while tie connectivity conditions were set between the others. In the second case (VCCT-2), the delamination conditions were used for each contact pair of plies (11 in total). The Fig. 12 demonstrates how these two model setups with varying number of the contact pairs may influence the numerically simulated mechanical response of the macroscopic sample. The 22  strains values at the control point depending on the displacements of the narrow part of the specimen were compared. The control point was taken in the middle of back surface of the wide part of the specimen. It was established that the processes of delamination inside the specimen can begin at the early stages of deformation, which is reflected in the results obtained by models with different number of separated plies. Due to the fact that delamination propagates between different pairs of plies unevenly, there are variations in the control values of strains on the surface of the specimen. As the VCCT-1 model consisted of more contact pairs with allowable delamination, the realization of the applied load leaded to redistribution of stresses between all these plies and triggered the delamination onset criterion for every contact pair. Since there were only five contact pairs with allowable delamination in VCCT-2 model, it had become more rigid, what consequently resulted in later initiation of delamination between contact pairs. Performing such numerical test may be useful for the tasks where the direct experimental comparison is performed, since in many works the VCCT delamination models are simplified to lower number (sometimes just one) of contact pairs, which may be the reason for inaccurate comparison results.

CONCLUSIONS
n this work, the behavior of an L-shaped specimen of a laminated polymer composite material under the influence of a load, providing development of delamination and accumulation of damage in individual plies, was studied numerically. To simulate delamination, a ply-by-ply models with various numbers of contact pairs were investigated, for which the virtual crack closure technique using the BK criterion was implemented. For the degradation of the properties of the plies, a model of progressive failure was applied to the plies using two criteria: the Hashin criterion and the multicomponent criterion. The mean field method with Mori-Tanaka approach was used in order to analyze the behavior of microstructural components in the numerical implementation of delamination and damage accumulation models. The obtained numerical results allow to trace the influence of model parameters on the mechanical behavior and fracture of the specimen both at the ply level and at the micro-scale level, while simultaneously taking into account delamination and damage accumulation in the ply. It has been established that the degradation of the plies' properties, geometric parameters and the parameters of the progressive failure model may lead to uneven numerically simulated growth of delamination. This indicates the importance of selecting criteria and conducting preliminary experiments to evaluate the strength constants in order to reliably reflect the fracture processes in models of composite structures using the interconnection of scale levels.