1
A COMPUTATIONAL MODEL FOR PREDICTION OF PROGRESSIVE
DAMAGE IN LAMINATED COMPOSITES
A. Ahmed1,2*, L. J. Sluys1 1
Department of Civil Engineering and Geoscience, Delft University of Technology, The Netherlands
2
Department of Civil Engineering, N-W.F.P University of Engineering & Technology, Pakistan *A.Ahmed@tudelft.nl
Keywords: Progressive damage, Discontinuous solid-like shell, Mesh independent matrix
cracking, Delamination
Abstract
A finite element model based on solid-like shell elements is presented for the simulation of progressive damage in laminated composite structures. To model mesh independent matrix cracking, a discontinuous solid-like shell element (DSLS) is utilized. The shell element has only displacement degrees of freedom, thus avoids the need for a complicated update of rotation degrees of freedom in nonlinear applications. To model delamination phenomena, a shell interface model is presented. The model allows computationally efficient simulation of delamination. To model the coupled response of matrix cracking and delamination under large deformations, a computational framework is developed. The combined modeling of matrix cracking and delamination is achieved without incorporation of additional degrees of freedom. In addition to physical nonlinearities, the numerical model is also able to simulate geometrical nonlinearities. Numerical examples are presented to simulate failure resulting in cracking and delamination in laminated composites.
1 Introduction
The increasing use of fiber-reinforced composites in modern day world has also increased the significance of performing failure analysis for reliable and safe design of fiber reinforced laminated composite structures. This necessitates a need for developing robust and efficient numerical tools, able to properly take into account the progressive damage in these structures. For instance, impact on composite structures causes significant damage in terms of matrix cracking and delamination. It has been experimentally observed, e.g. [1,2,3], that these two damage mechanisms appear concurrently and there is a strong interaction between them. However, complexities exist in developing efficient numerical tools due to the presence of different failure mechanisms and their interactions. In several studies, e.g. [2,4,5,6,7], failure based or continuum damage models have been proposed for the prediction of impact damage in laminated composites. However, such models do not always perform well in localization problems, which may lead to mesh dependent results and a wrong failure mode.
In order to better represent localization phenomena at the interface, cohesive zone models were also explored to model delamination phenomena. [8] used interface elements with cohesive zones to model delamination phenomena in fiber reinforced laminated composites subjected to in-plane loading. [9] modeled each ply of the laminate with a single layer of solid
2
elements. The process of delamination cracking was modeled by doubling the nodes at the interface, while matrix cracking damage was modeled with a continuum damage model. Although, the use of solid elements helps in obtaining the three dimensional stress field, which is crucial for delamination onset and propagation, these elements tend to lock (Poisson thickness locking) and create numerical difficulties when used in thin shell applications. [10,11] presented a failure model based on solid-like shell elements and used a plasticity based approach to model matrix cracking and interface elements for delamination cracking. The process of matrix cracking also results in strain localization and therefore continuum damage models or plasticity models create numerical difficulties in finite element computations. [12] and [13] used interface elements to model both matrix cracking and delamination to simulate the in-plane and out-of-plane damage in laminated composites, respectively. However, the use of interface elements requires the finite element mesh to be aligned with the crack geometry and the cracks can only grow along predefined locations. A different approach to model cracking in a material is to use the partition of unity approach [14], which allows modeling of arbitrary propagating cracks through the finite element mesh. Such class of methods has been explored for modeling the in-plane response of composite laminates e.g. [15,16]. For modeling mesh independent matrix cracking in laminated plates and shells, [17] presented a discontinuous shell model based on the phantom node method [18].
2 Progressive failure model
In this contribution a meso-scopic failure model for laminated composite plates and shells is presented. Two key damage mechanisms, i.e. matrix cracking/splitting and delamination damage are considered in this work. A computational framework is developed to take into account the coupled response of matrix cracking and delamination during events of impact damage, see figure 1. Owing to the three dimensional nature of damage, especially delamination damage, solid-like shell elements are used to model each ply with orthotropic material properties. This not only provides a three dimensional state of stress for failure analysis but also removes Poisson thickness locking, commonly found in solid elements, and thus avoids the need of a high degree of mesh refinement in thickness direction to model thin plies of laminates. The progressive failure model allows for arbitrary propagation of matrix cracks through the finite element mesh. Delamination cracking is modeled using a shell interface model. Details of mathematical models to simulate each failure mechanism are given below.
2.1 Matrix cracking/splitting damage
To simulate matrix cracking, a discontinuous solid-like shell element (DSLS) [19] is used. The discontinuity in the shell mid-surface, shell director and more importantly in the internal stretching field is incorporated by exploiting the phantom node method [18]. This enables the element to model arbitrarily propagating cracks through a finite element mesh. In [17], the discontinuous solid-like shell element (DSLS) has been successfully used to model matrix cracking in laminated fiber reinforced shell structures. The crack growth direction is taken equal to the fiber direction.
2.2 Delamination damage
To model delamination cracking, a finite element method for delamination proposed by [20] is utilized, which is an alternative approach for modeling interfacial phenomena compared to traditional interface elements. This method allows for complete kinematic description of
3
interfaces as opposed to interface elements. As a result, it becomes possible to obtain a fully consistently linearized tangent for the interface contribution, which is important for quadratic convergence of the Newton-Raphson scheme. Particularly, for geometrically nonlinear problems, where large changes in normal vector and area of the cohesive crack surfaces, require computation of cohesive geometric stiffness [21,26]. In case of traditional interface elements, this cannot be achieved due inadequate kinematic description of the interfaces.
Figure 1. Progressive failure model
(a) Adapted interface (b) Un-adapted interface
Figure 2. Kinematics of adapted and un-adapted interface 2.3 Progressive damage - Numerical aspects
2.3.1 Coupling between matrix cracking and delamination
The presence of a matrix crack in one or both plies connected to a common interface results in a discontinuous interface (see figure 1). In order to properly take the interaction between matrix cracking and delamination into account, one or both planes of the interface have to be updated. In case of bending dominated problems, if the cracked interface is not updated (see figure 2), it will result in incorrect computation of normals to the interface and therefore will result in incorrect computation of damage. As a consequence the load capacity will be over-predicted, as was demonstrated in [23].
4
In order to model an adapted interface, a partition of unity approach can be used. However, due to the fact that interfaces in our shell interface model are defined as an integral part of the continuum elements, advantage can be taken from the database which has already been generated for matrix cracking. This will automatically update the connectivity of the elements connected to an interface. As a consequence, no additional model and degrees of freedom are added to incorporate the discontinuity in the two planes of the interface.
2.3.2 Crack growth under large deformation
In cohesive zone models the tractions are transferred through a unique crack surface, but when a body undergoes large deformation, no unique crack surface can be defined. Usually the non-uniqueness in defining a normal to an interface is avoided by defining an average crack surface, [21]. However, this assumption leads to incorrect kinematic description of the interface. This can be demonstrated through a simple example. Consider an interface (figure 3), which is given a rigid rotation.
Figure 3. Cracked surface under large deformation/rotation
The deformation gradients of the interfaces CD, C′D and a fictitious average interface E′′ D are given as:
− = = − = ′ ′′ 5 . 0 0 5 . 0 0 1 0 5 . 0 0 5 . 0 , 1 0 0 0 1 0 0 0 1 , 0 0 1 0 1 0 1 0 0 D E D C CD F F F
Accordingly, the area ratios for the three surfaces are
7071 . 0 , 0 . 1 , 0 . 1 = = = ′ ′′ dA da dA da dA da D E D C CD Where D E D C CD da da da ′′ ′ and
, are the deformed areas of the three surfaces and dA is the reference crack surface area. Above computation shows that there is a reduction of the deformed area, which is un-realistic for the problem under consideration. The reasons for such a behavior are obvious and are also schematically illustrated in figure 3. It can be inferred from the figure that an average kinematic assumption for the interface, inherently assumes that the interface will not follow the curved path. The correct average crack surface should have been E′ , instead crack surface D E′′ is predicted and the interface front is located at D
5
the chord. This results in a decreased interface length lc = 0.3536. However, it is worth to note
that the magnitude of the error induced by using averaging kinematics is also dependent upon the geometry, the boundary conditions of the body and the magnitude of interface rotation. 2.4 Constitutive laws
In meso-scopic analysis, where each ply of the laminate is modeled with a single layer of solid-like shell elements in thickness direction, the ply is considered to be of homogeneous material with orthotropic material properties. A cohesive constitutive law based on the cohesive law of Xu and Needleman [22] is used, both for matrix and delamination cracking. The mode-mixity is taken into account by the Benzeggagh-Kenane [24] mode-mixity criterion. More details on the cohesive model can be found in [23].
3 Numerical Results
3.1 Peel test
A peel test is performed to demonstrate the rapid convergence and efficiency of the shell interface model. A similar type of problem was also analyzed by [21]. The geometry of the model is shown in figure 4.
Figure 4. Peel test model geometry
A displacement controlled analysis is performed. A total displacement of 4mm is applied in 20 equal increments. The bulk material is considered to be isotropic, homogeneous with Young’s modulus, E = 100Mpa and Poisson ratio equals to zero. A bi-linear isotropic cohesive law is assumed for the cohesive crack. The cohesive strength is assumed to be 3MPa and fracture toughness is 0.74 N/mm. The performance of the present model is judged by comparing the relative residual norms ( R/Ro) at a particular displacement step with the results
obtained without the contribution of cohesive geometric tangent stiffness. Ro is the residual
norm of the initial iteration. The solution is considered to converge when the relative residual norm becomes less 1.0e-10. Table 1 compares the residual norms at various load steps. Figure 5 shows the deformed shapes at different displacement increments. Note that the deformation is not scaled. Iter; No.
u
=
1
u
=
2
u
=
3
u
=
4
W it h o u t G eo -s ti ff1 6.4360e-01 5.4281e-01 4.8312e-01 7.2589E-01 2 1.5532E-02 2.8160E-02 3.4767E-02 2.5194E-02 3 2.8335E-04 2.8160E-02 3.7535E-03 3.9396E-03 4 1.6041E-05 3.9804E-03 1.6415E-04 5.7263E-04 5 8.4909E-07 3.5507E-04 1.2344E-05 6.3328E-05 6 4.4167E-08 4.1491E-05 8.7764E-07 8.9142E-06 7 2.2913E-09 5.4756E-06 6.1688E-08 1.2855E-06 8 1.1883E-10 7.1784E-07 4.3281E-09 1.8603E-07 9 6.1623E-12 9.3925E-08 3.0360E-10 2.6934E-08 10 --- 1.6069E-09 2.1297e-11 3.8998e-09 11 --- 2.1018E-10 --- 5.6466e-10 12 --- 2.7491E-11 --- 8.1759e-11
6
Ro 7.6992e-01 5.5356e-01 4.5706e-01 5.0098e-01
W it h G eo -s ti
ff 1 2 6.44E-01 1.51E-02 5.65E-01 2.57E-02 5.06E-01 3.26E-02 8.04E-01 3.38E-02
3 3.95E-05 5.54E-03 1.75E-03 1.07E-02 4 1.71E-09 1.29E-04 1.12E-06 1.03E-03 5 3.59E-14 6.51E-09 3.16E-12 1.57E-06
6 1.2313e-13 2.0949e-11
Ro 7.6948e-01 5.5957e-01 4.6187e-01 5.7969e-01
Table 1. Relative residual norms after every iteration in different displacement increments
It can be observed that the shell interface model converges rapidly and a quadratic convergence is achieved. The number of iterations required to obtain a converged solution is almost half compared to the analysis without the cohesive geometric stiffness contribution, which is usually the case with traditional interface elements.
Figure 5. Deformed shapes at displacement increments, u=1,2,3and4
3.2 Two-ply laminated plate,[0/90]
A two-ply [0/90], square laminated plate with length 20mm, loaded quasi-statically with a center point load, is analyzed (see also [23]). The thickness of each ply is 0.2mm. The material is considered to be carbon-epoxy. Due to symmetry, only one-half of the plate is modeled. Figure 6 shows the damage progress at different load steps. It can be observed that progressive damage in terms of initiation and growth of matrix cracks and simultaneous development of delamination is modeled properly. Note that the two processes i.e. progress of a traction free crack and delamination, are growing side by side, signifying strong interaction between the two mechanisms. A typical peanut-shape delamination area, as observed during experiments [25], is predicted.
4 Conclusions
A progressive failure model for damage in laminated composite plates is presented. The model uses solid-like shell elements which, on one hand, are able to model thin plies of the
7
laminate and, on the other hand, give a complete three dimensional state of stress. This is crucial for delamination damage. The model is capable of simulating mesh independent matrix cracking. In addition to this, a strong interaction between the matrix crack and delamination damage is captured.
Figure 6. Damage progression in a two-ply laminated plate subjected to quasi-static transverse load
Acknowledgements
Financial support from N-W.F.P University of Engineering & Technology Peshawar, Pakistan, under the HEC approved project titled ``Strengthening of existing earthquake engineering center'' and Delft University of Technology, Netherlands is gratefully acknowledged.
References
[1] Joshi S. P., Sun C. T., Impact induced fracture in a laminated composite. J. Compos. Mater.
19, p 51-66 (1985).
[2] Choi H. Y., Chang FK., A model for predicting damage in a graphite/epoxy laminated composites resulting from low-velocity point impact. J. Compos. Mater. 26, pp. 2134-2169 (1992).
[3] Lammerant L., Verpoest I., The interaction between matrix cracks and delaminations during quasi static impact of composites. Compos. Sci. Technol. 51, pp. 505-516 (1994).
[4] Hou J.P., petrinic N.,Ruiz C., Hallet S.R., Prediction of impact damage in composite plates. Compos. Sci. Technol. 60, pp. 273-281 (2000).
8
[5] Zhao G., Cho C., On impact damage of composite shells by a low velocity projectile. J. Compos. Mater. 38, pp. 1231-1254 (2004).
[6] Iannucci L., Ankersen J., An energy based damage model for thin laminated composites. J. Compos. Mater. 66, pp. 934--951 (2006).
[7] Donadon M.V., Iannucci L., Falzon B.G., Hodgkinson J.M., de Almeida S.F.M., A progressive failure model for composite laminates subjected to low velocity impact damage. Compos. Struct. 86, pp. 1232--1252 (2008).
[8] Yang Q., Cox B., Cohesive models for damage evolution in laminated composites. Int. J. Fact. 133, pp. 107-137 (2005).
[9] Collombet F., Bonni J., Lataillade J.L., A three dimensional modelling of low velocity impact damage in composite laminates. Int. J. Numer. Methods Engng. 39, pp. 1491-1516 (1996). [10] Hashagen F., Schellekens J.C.J, de Borst R., Parisch H., Finite element procedure for
modelling fibre metal laminates. Compos. Struct. 32, pp. 255-264 (1995).
[11] Hashagen F., de Borst R., de Vries T., Delamination behavior of spliced fiber metal laminates. Part 2. Numerical investigation. Compos. Struct. 46, pp. 141--162 (1999).
[12] Wisnom M.R., Chang F., Modelling of splitting and delamination in notched cross-ply laminates. Compos. Sci. Technol. 60, pp. 2849-2856 (2000).
[13] Bouvet C., Castanie B., Bizeul M., Barrau J.J, Low velocity impact modelling in laminate composite panels with discrete interface elements. Int. J. Solids Struct. 46, pp. 2809--2821 (2009).
[14] Melenk J.M., Babuska I., The partition of unity finite element method: Basic theory and application. Comput. Methods Appl. Mech. Engrg. 139, pp. 289--314 (1996).
[15] Iarve E.V., Mesh independent modelling of cracks by using higher order shape functions. Int. J. Numer. Methods. Engrg. 56, pp. 869-882 (2003).
[16] van der Meer F.P., Sluys L.J., A phantom node formulation with mixed mode cohesive law for splitting in laminates. Int. J. Fract. 158, pp. 107-124 (2009).
[17] Ahmed A., van der Meer F. P., Sluys L. J., A geometrically exact, discontinuous shell model for transverse matrix cracking in composite laminates in “Proceeding of ECCOMAS thematic conference on the mechanical response of composites”, Hannover, Germany, pp. 371-378, (2011).
[18] Hansbo A., Hansbo P., A finite element method for the simulation of strong and weak discontinuities in solid mechanics. Comput. Methods Appl. Mech. Engrg. 195, pp. 3523-3540 (2004).
[19] Ahmed A., van der Meer F. P., Sluys L. J., A geometrically nonlinear discontinuous solid-like shell element (DSLS) for thin shell structures. Comput. Methods Appl. Mech. Engrg. 201-204, pp. 191-207 (2012).
[20] Larsson R., A discontinuous shell-interface element for delamination analysis of laminated composite structures. Comput. Methods Appl. mech. Engrg. 193, pp. 3173-3194 (2004). [21] Wells G.N., de Borst R., Sluys L.J., A consistently geometrically nonlinear approach for
delamination, Int. J. Numer. Methods Engng.,54, pp. 1333-1355 ( 2002).
[22] Xu X.P., Needleman A., Numerical simulations of fast crack growth in brittle solids. J. Mech. Phys. Solids 42, pp. 1397-1434 (1994).
[23] Ahmed A., Sluys L.J., A three dimensional progressive failure model for impact damage in laminated composites. (In preparation).
[24] Benzeggah M.L., Kenane M., Measurement of mixed-mode delamination fracture toughness of unidirection glass/epoxy composites with mixed mode bending apparatus. J. Compos. Sci. Technol. 56, pp. 439-449 (1996).
[25] Liu D., Impact induced delamination-A view of bending stiffness mismatching. J. Compos. Mater., 22, pp. 674-692 (1988).
[26] Gasser T.C., Holzapfel G.A., Modeling 3D crack propagation in unreinforced concrete using PUFEM. Comput. Methods Appl. Mech. Engrg., 194, pp. 2859-2896 (2005).
