A model for fast delamination analysis of laminated com- posite structures

Delamination is one of the major failure mechanisms for composites and traditionally the simulation requires high expertise in fracture mechanics and dedicated knowledge of the Finite Element Analysis (FEA) tool. Yet, the simulation cycle times are high. Geometrically nonlinear analysis approach, which is based on the Reissner-Mindlin-Von Kármán type shell facet model, has been implemented into the Elmer FE solver. Altair ESAComp software runs the Elmer Solver in the background. A post-processing capability, which enables the prediction of the delamination onset from the FEA output, has been implemented into the Altair ESAComp software. A Virtual Crack Closure Technique (VCCT) specifically developed for shell elements defining the Strain Energy Release Rate (SERR) related to the different delamination modes at the crack front is used. The onset of delamination is predicted using the relevant delamination criteria that utilize the SERR data and material allowables in the form of fracture toughness. The modeling methodology is presented for laminates including initial through-thewidth delamination. Examples include delamination in the solid laminate and debonding of the skin laminate in the sandwich structure. Rather coarse FE mesh has proved to yield good results when compared to typical approaches that utilize the standard VCCT or Cohesive Zone Elements.


Introduction
Layered composite structures exhibit various failure mechanisms namely fiber failure, matrix cracking, buckling and delamination. In this work the focus is in delamination, which is a separation between internal layers of a composite laminate caused by high through-the-thickness stresses, impacts or manufacturing defects. Delamination causes a significant structural damage, particularly in compression [1,2].
There is a clear need from industry to have easy-to-use fast delamination assessment tools for both the initial design phase and design verification. Currently, assessment of delamination requires high expertise and dedicated knowledge of particular Finite Element  [4,7].
Analysis (FEA) programs. The motivation of the work was to develop an approach that brings the delamination assessment in to the hands of a typical composites engineer. Also, the solution works as a benchmark tool to provide confidence for the users' of the general purpose FEA packages.
The solution to be developed shall be very responsive, i.e. modeling and computation time shall be relatively small so that efficient what-if-studies can be performed. On the other hand, the provided solution shall not compromise with respect to the accuracy of the results. This kind of technology can be called "quick-and-clean". The solution shall be rather general so that it serves in wide range of applications and potential future extensions shall be kept in mind.
Typical delamination shapes are circular, elliptic and peanut [3]. The simplified circular and elliptic idealized shapes of delaminations [3] were obtained by the algorithms for finding the smallest enclosing circle and smallest enclosing ellipse. The effect of delamination shape on the buckling load is presented in references [4,5,6].
Various types of delaminations exist like free edge, through-the-width and embedded, see Figure 1 [4,7]. At free edges shell elements are not capable of describing the 3D stress state. Through-the-width delaminations appear to be the right starting point due to the simplicity of the geometry and available reference data. In the scope of this work embedded delaminations seem to be reasonable target that is considered as a future extension. Virtual Crack Closure Technique VCCT [8] is the choice of this work since it is generally applicable [9]. It is based on Linear-Elastic Fracture Mechanics (LEFM). When predicting onset of delaminations growth material nonlinearities can be neglected, which is desirable. For most of the Fiber-Reinforced Polymer (FRP) structures this is applicable approach.
In the simplest form the classic application of the VCCT consists of post-processing the finite element results of a layered structure in which an initial crack is explicitly modeled. If the study is limited to the onset of delamination growth no special elements or constrained equations are needed. Some approaches involve placing of stiff spring elements at the crack front and beam results are used to determine the forces at the crack front.
In the modeling approach, single shell elements are used at the intact regions and two overlaid shell elements at the delaminated regions. Thus, single delamination in the through-the-thickness direction can be modeled. Nodes of the FE model reside at the plane of the delamination. Loads are given with respect to the geometrical mid-plane, and then internally transformed to the reference plane.
In this work, shell facet elements are used. Stress resultants are used together with in-plane strains and curvatures. The simulation is performed with a shell facet model implemented in the Altair ESAComp software [10] The shell facet model was implemented with the Elmer opensource Finite Element Method (FEM) solver [11,12,13] developed by CSC -IT Center for Science (CSC) in collaboration with Finnish universities, research laboratories, and industry.
ESAComp development started in 1992 at Helsinki University of Technology (now known as Aalto University) as a project initiated by the European Space Agency (ESA). In 2000 Componeering (currently Altair Engineering Finland), a Helsinki based company founded by the original project team, took over the development. ESAComp is a proprietary licensed software owned by Altair Engineering. ELMER technology has been integrated in ESAComp for the realization of the FE analysis capability with an agreement between Componeering and CSC.

Reissner-Mindlin-von Kármán type shell facet model
The plate bending problem is formulated for a thin or moderately thick laminated composite plate which in its undeformed configuration occupies the region Ω × (−t/2, t/2), where Ω ⊂ R 2 is the midsurface and t > 0 is the laminate thickness. The kinematical unknowns in the model are transverse deflection w, in-plane displacement u = (u x , u y ), rotation of the shell reference surface β = (β x , β y ), and drilling rotation ω. The plate is subjected to the in-plane load f = (f x , f y ) and the transverse pressure g at its reference plane.
Let us note that our model is actually a numerical "shell facet model" since we are in fact considering one element in the mesh. As far as we know, there is no mathematical analysis guaranteeing the consistency of this very classical engineering approach, but it seems to work fine in practice. We will use standard notation of tensor calculus. Dyadic and index notation with summation convention over repeated indices are used in parallel. Latin indices take their values in the set {1, 2, 3} and Greek indices in the set {1, 2}.

Constitutive relation for a single layer
Let us denote by e i andē j the cartesian basis vectors for the so called 123-coordinate system of a single ply, and for the xyz-system of material coordinates common to all plies, respectively. In the material coordinate system, i.e., the laminate coordinate system, the layer system has been rotated by a positive counter clockwise angle θ about the zaxis. Hence, we define the transformation matrix between the two coordinate systems as For linear orthotropic materials, a plane stress state is assumed and the constitutive relation for each ply has the form σ = Q : where σ = σ ij = σ ji is the second order stress tensor, = ij = ji is the strain tensor, and Q = Q ijkl = Q jikl = Q ijlk = Q klij is the fourth order tensor of elastic stiffness coefficients.
In the laminate coordinate system the constitutive equation is written as whereσ ij = T ip T jq σ pq is the laminate stress,¯ ij = T ip T jq pq is the laminate strain, and Q ijkl = T ip T jq T kr T ls Q pqrs is the tensor of stiffness coefficients in the laminate coordinate system. The six independent non-zero components of Q are computed using the orthotropic material engineering constants E 1 , E 2 , ν 12 , ν 21 = ν 12 E 2 /E 1 , G 12 , G 23 , and G 31 [14] as

Kinematic relations for a laminate
The kinematic relations for a laminate are considered in the xyz-coordinate system. For notational simplicity, laminate stresses and strains in the xyz-coordinate system are in the following denoted without bar symbol.
Let us note that the original von Kármán strains [15] do not have the quadratic in-plane displacement gradients in the membrane tensor ϕ. The quadratic in-plane displacement gradients are considered in the model for the sake of completness and they are not expected to enhance the accuracy of the original model significantly.

Constitutive relations for a laminate
In plane-stress state, the laminate membrane stress resultants N (forces per unit length) and bending moment resultants M (moments per unit length) are obtained by integration of the stress resultants of all layers z k−1 < z < z k , k = 1, . . . , n, over the thickness of the laminate as Furthermore, the resultant transvere shear forces S are obtained from Using the constitutive equation and the kinematic relations of Reissner, Mindlin, and von Kármán we get the following constitutive relations for the laminate The tensors A, B, and D are defined according to the Classical Lamination Theory (CLT) [14] as whereQ (k) defines the constitutive relation for linear orthotropic materials in plane stress state for layer k in the laminate coordinate system. The tensor for transverse shear stiffness A can be defined according to the CLT [14] as In ESAComp [10], the computation of the out-of-plane shear stress distribution and stiffness is based on the theory developed at German Aerospace Center (DLR) [16]. First Order Shear Deformation Theory (FSDT), i.e., the out-of-plane shear deformation is considered as an extension to the CLT plane stress assumption in ESAComp.

The shell facet model
The functions u, w, β, ω are determined from the condition that they minimize the potential energy of the plate. The energy is defined as +C where C > 0 is a penalty parameter for imposing the condition ω = rot(u) (see [17]), and Substituting the constitutive equations in Eq. 20, we get +C The finite element implementation of the model has been discussed by the authors in [13]. The elements used in this paper are either triangular or quadrilateral and use (bi)linear basis for all field variables. The equilibrium path is computed numerically by the arc-length method.

Principle of VCCT
In VCCT the Strain Energy Release Rate SERR is calculated at the delamination front. In VCCT-based theory it is assumed that SERR is equal to the work per unit area required to close a crack over a small distance. VCCT directly provides SERR data related to each delamination mode. Here, VCCT plate theory is used. Unlike the standard VCCT it provides the resultant SERR. When the plate theory is used SERR is the work required to change the mid-plane strains ε(u) and curvatures ε(β) in the cracked part of the plate ζ = 1, 2 at the crack front so that they are equal to the mid-plane strains and curvatures in the uncracked part of the plate at the crack front [18].
The benefit of the VCCT plate theory is that SERR can be directly derived from shell element results without introducing additional elements. These results cover resultant forces N (u, w, β) and moments M (u, w, β) that are output at the element reference plane, which is the delamination plane. Resultant forces and moments are converted to the geometrical mid-planes of each three laminate section. Then, CLT is used to calculate strains and curvatures at the geometrical mid-planes. Geometrically nonlinear analysis approach is needed to determine the nonlinear relation of the load-displacement path. However, the material is assumed to stay in the linear-elastic regime and therefore, CLT is applicable for the determination of mid-plane strains and curvatures from the FEA stress resultants. The goal is to determine the change in mid-plane strains and curvatures at the different sides of the crack. Therefore, at the virgin laminate strains and curvatures at the levels defined by the mid-planes of the cracked laminate needs to be solved. This is obtained by assuming a constant curvature in the virgin laminate. Once the strains and curvatures at both sides of the crack at the desired levels are known, the change can be determined. Laminate structures are basically the same at the different sides of the crack and knowing the change in strains and curvature leads directly to the change in forces and moments when using CLT. Finally, SERR is determined by multiplying changes in resultant forces and moments by changes in strains and curvatures and summing up both laminates. The procedure is well defined in [9,19] and described in Figure 2.
The equation for SERR is Described approach defines the total SERR (see Figure 2), but it is not capable to make distinction between the different delamination modes. Nevertheless, this is very important piece of the general solution.

Crack tip element
In [19] a three-dimensional crack tip element is presented and used to determine SERR and mode mix for different types of laminated plates containing delaminations. The approach essentially contains two parts; first so-called concentrated crack tip force and moment are determined and then, the mode I and mode II components (mode-mix) are obtained. When determining the crack tip force and moment, strains ε 22 (u), ε 22 (β), and ε 12 (β) are enforced to zero next to the crack tip and consequently, only strains ε 11 (u), ε 12 (u), and ε 11 (β) are non-zero. Thus, the problem is reduced to consider only modes I and II. In the second part the mode mix is determined using the obtained crack tip force and moment and with an experimentally determined expression for the mode mix parameter. The method is called CTE/NSF, where NSF stands for non-singular field. In standard VCCT oscillatory singularity exists for modes I and II and expertise is needed to define appropriate mesh density to obtain reliable results. The use of the nonsingular method diminishes the element size factor and it appears that CTE/NSF is more appropriate for general use for predicting delamination growth than any of the other methods that have thus far been proposed.
With CTE/NSF(for detailed description refer to [19]) SERR values related to modes I (G I ) and II (G II ) are solved. Similar to the VCCT plate theory the starting point is the resultant forces N (u, w, β) and moments M (u, w, β) extracted from the FE solution. Using the above assumption for enforced strains and curvatures and laminate reduced stiffness matrix (analogy to CLT) crack tip force and moment can be solved. The mode mix parameter is dependent on the material properties and laminate lay-ups. Once G I , G II and G (Eq. 23) are known, G III is solved from The method has been benchmarked with various loading and laminate configurations, see next chapter. For mode III delamination the method gives slightly too optimistic results, but mode III is not a governing mode for typical engineering problems.

Benchmark
Part of the software verification procedure involved setting up a benchmark model with various laminate and loading conditions. The verification cases are shortly described, and the results are presented. It should be noted that the actual implementation of the ESAComp delamination module includes different structural components than the one used in the benchmark. Figure 2. The procedure to determine SERR with the plate theory [9,19]. The FE shell model for three regions is defined so that location of the delamination is considered. The reference plane (1) is at the level of the delamination. Associated node offsets are used. From the FEA solution results N (u, w, β), M (u, w, β), ε(u), ε(β) are extracted at the geometrical mid-planes (2) of each section. Strains are defined at this level (3) assuming a constant curvature. The change of mid-plane strains ∆ε(u) and curvatures ∆ε(β) is calculated between the planes 2 and 3. The change of resultant forces and moments ∆N, ∆M , respectively, is calculated with the help of of ∆ε(u) and ∆ε(β). SERR is calculated using equation Eq. 23. An FE model very similar to [19] was constructed (see Figure 3). There were small discrepancies between the reference model and the model used in this work, however. The reference shell model was constructed using 8-noded shear deformable elements where as in the current work 4-noded shear deformable elements were used. The reference model contained two identical, but separate shell meshes that were 12 units away from each other in the thickness direction. In the reference model, the shell element reference plane was in the geometrical mid-plane of the laminate. In the uncracked region of the model, constrained equations were used to tie displacements and rotations of the node pairs. In the current approach, the uncracked region was modeled with single shell elements. The reference plane of the shell element was in the plane of the crack (see Figure 3, z = 0).
The cracked regions were modeled with two separate shell meshes. Similar to the uncracked region, the reference plane of the shell elements was in the plane of the crack. The model was clamped from the edge where x = 0. Three different loading conditions were used in [19]; an opening load, in-plane shearing load and an out-of-plane shearing load. The opening load was applied using a bending moment. For the two shearing load cases the plane of the load introduction was different in the reference work when compared to the current approach. Therefore, correction moments were needed to obtain equivalent load cases.
In this study, the total force or total moment applied to the edge of the cracked plane was 1. The study involved four different laminate configurations. In all studied cases, the location of the crack was in the same level. Studied configurations are summed up in Figure 4 and comparative results are presented for two verification cases in Figures 5  and 6. ESAComp results are presented with red curves and those should be compared to CTE/NSF-results.

Material properties
VCCT requires four material properties: fractures toughness in main directions and mixed mode interaction parameter. Test methods for the definition of the fracture toughness are illustrated in Figure 7. Mode I fracture toughness G IC is measured with Double Cantilever Figure 4. Verification Cases VC1 to VC9. In laminate code "d" indicates the location of the delamination. Component indicates the total load per cracked side. Correction moments are not presented for the shearing load cases. "k" stands for coefficient by which the load (Component) was multiplied while generating the results presented in figure 5 and 6. SERR I, II and III stand for the associated Strain Energy Release Rate components along the crack tip in the width direction for the different VCs. Comparison of mode SERR I and mode SERR II for VC1 are shown in Figures 5 and 6, respectively.
Beam (DCB). The corresponding standard is ASTM D5528. Mode II fracture toughness G IIC is measured with End-Notched Flexure (ENF) test as proposed by Martin [20], which however, appears to be obsolete. There is a brand-new Mode II ASTM standard ASTM D7905/D7905M-14 based on the 3ENF specimen (see Figure 7). The Europeans have developed ISO 15114:2014 based on the calibrated end-loaded split specimen C-ELS. For Mode III fracture toughness G IIIC test standard is not available. For G IIIC the value of G IIC is proposed, which is a conservative estimate [9].
Mixed-mode interaction parameter η is obtained from a curved fit of the BK criterion using all the results from the DCB, MMB and ENF testing. Standardized fracture mechanics tests are all for characterization of a 0-0 interface in a unidirectional laminate, whereas typical aerospace skin-stiffener interface would commonly involve plies at different angles [21]. This gives a conservative estimate for the fracture toughness in off-angle interfaces and is generally used.

Failure criteria
A widely used criterion is B-K by Benzeggah and Kenane [2,22]. The power-law shall be used when the mixed-mode interaction parameter is not available (α = β = 1). Moreover, G IIC = G IC will be usually a conservative assumption, if the value for G IIC is not available [2]. The most recent criterion is proposed by Reeder [23,24] The onset of delamination criterion is predicted using the failure index: The three above criteria are referred as the typical criteria [25]. Figure 8 shows predictions of the different mixed-mode delamination propagation criteria and corresponding experimental data for CFRP IM7/977-2 [9].

Delamination of a solid laminate
The first ESAComp example repeats the study described in [26]. A solid laminate with explicitly defined through-the-width delamination is loaded in compression (see Figure  9). The onset of delamination is typically driven by the local buckling and mode I is the most critical. In the reference, geometrically nonlinear behavior initiates at strain level of 0.18 %. Respectively, the onset of delamination begins at 0.22 %.
ESAComp geometrically nonlinear analysis was run with 20 sub-steps. ESAComp Result tracker indicates during the solution phase how the maximum resultant displacement evolves as a function of the load increments (see Figure 10). Delamination analysis is performed after the last sub-step and ESAComp prediction for the onset of delamination is well in line with the reference. B-K delamination criterion yields f = 1.1, which indicates delamination propagation. This is shown in Figure 10. Line charts illustrate the relation of G I /G Ic , G II /G IIc and f (B − K).

Debonding of a solid skin laminate
In the second ESAComp example a thick sandwich structure containing a through-thewidth debond of the skin laminate is subjected to compression (see Figure 11). The reference study is described in [27]. ESAComp delamination analysis predicts the onset of delamination with the load level of 10.3 kN. Experimental tests indicated that the debonding propagation starts at 9.3 kN.
The initial delamination amplitude corresponded to 0.2 mm and the shape of the delamination resembled a half-sine wave with the maximum deflection at the center consistent with Digital Image Correlation measurements. Identical imperfection was considered in the ESAComp analysis. The critical load for the onset of delamination was rather sensitive to the imperfection amplitude and with amplitude of 0.5 mm the critical load level decreased to 9.6 kN.
It should be noted that at the load level of 9.0 kN maximum strain level in compression was over 1.4 % at the debonded skin laminate. It can be expected that some material degradation may have been taken place before the onset of delamination, which was not considered in the simulation.
In this specific case the critical delamination mode is opening. Benchmark studies (refer to Figure 5) have been made for various laminate and load configurations. Generally, CTE/NSF method gave ∼5 % too optimistic delamination onset loads when compared to the full 3D model and standard VCCT. This partly explains slightly too optimistic prediction for the onset of delamination.

Conclusions
Conducted work presents fast and accurate approach for the prediction of onset of delamination for solid and sandwich laminates that contain explicitly defined through-the-width delamination. Mode decomposition compares well with the results obtained with 3D elements and using VCCT. Onset of delamination corresponds well against the model that utilized Reissner-Mindlin plate theory and the modified VCCT along the delamination front. Validation against experimental results shows good correspondence for the thick sandwich structure. In the current implementation the described solution is limited to single delamination in the through-the-thickness direction, but it can be extended to Figure 11. A thick sandwich structure containing a through-the-width debond of the skin laminate subjected to compression (left). FE mesh with contours (bottom right) shows the value of the Gtot at the delamination fronts. Line chart (top right) illustrate the relation of G I /G Ic at the delamination fronts. embedded delaminations. Single or multiple through-the-width and/or embedded delaminations could be defined at the same interface. For embedded delaminations contact elements shall be included to avoid inadmissible delamination modes.
ESAComp delamination module provides automated model creation for limited use cases. Solution time is at least two orders of magnitude smaller than for full 3D FE solution. All result items are easily accessible for the last load increment of the geometrically nonlinear analysis. ESAComp provides extensive material database including a collection of fracture toughness values for different fiber reinforced plastics. Together these aspects considerably lower the threshold to start simulation of delamination for laminates. For niche applications this approach provides complete solution to estimate the onset of delamination. The method does not predict the propagation of the delamination.