Overcoming the cohesive zone limit in composites delamination: modeling with slender structural elements and higher‐order adaptive integration

Cohesive element (CE) is a well‐established finite element for fracture, widely used for the modeling of delamination in composites. However, an extremely fine mesh is usually needed to resolve the cohesive zone, making CE‐based delamination analysis computationally prohibitive for applications beyond the scale of lab coupons. In this work, a new CE‐based method of modeling delamination in composites is proposed to overcome this cohesive zone limit on the mesh density. The proposed method makes use of slender structural elements for the plies, a compatible formulation with adaptive higher‐order integration for the CEs, and the corotational formulation for geometrically nonlinear analysis. The proposed method is verified and validated on the classical benchmark problems of Mode I, II, mixed‐mode delamination, a buckling‐induced delamination problem and a double‐delamination problem. The results show that elements much larger than the cohesive zone length can be used while retaining accuracy.

as pronounced as in the case of composites delamination. This is because the cohesive zone length decreases with smaller substrate thickness, higher substrate stiffness, and lower fracture toughness. 4 CFRP composite laminates possess all three, making the cohesive zone limit particularly stringent for the modeling of composite delamination. Several authors have given suggestions on the minimum number of linear CEs (ie, CEs conforming to first-order continuum/continuum shell ply elements) required to span the cohesive zone in order to ensure the prediction accuracy. Harper and Hallett suggested a number of 3, 6 while Moes and Belytschko proposed a number of 10. 3 Camanho and Davila used approximately three elements within the cohesive zone. 14 Similarly, Yang and Cox suggested to use from three to five elements. 4 The cohesive zone length for Mode I delamination in a composite laminate could be smaller than 1 mm, which means that the elements would be at most around 0.3 mm in size. This limitation often results in prohibitively high computational cost to perform simulations of delamination in component-scale composite structures. The above-mentioned problem of cohesive zone limit has drawn the attention of many researchers in the past. Turon et al 5 proposed an engineering solution to the problem by lowering the strengths of the CE, which enlarged the numerical cohesive zone length, thereby allowing larger CEs to be used. A possible consequence of this procedure, particularly with excessively reduced interfacial strengths, is that the simulation would predict a numerical cohesive zone much larger than the physical one, leading to overcompliant predictions of structural responses and underpredictions of the peak loads. 6,15 Dávila et al 16 developed a shell CE which was compatible with Mindlin-Reissner shell elements for plies. However, their method did not resolve the issue of cohesive zone limit. Yang et al 7 proposed to use Gauss quadrature coupled with level-set method and subdomain integration on linear CEs, which allowed the integration scheme to adapt to the partially damaged elements. Do et al 9 extended the previous study onto more types of integration improvements. Their approaches managed to enlarge the linear CE to be comparable to the cohesive zone length. However, the first-order shape functions were not well suited for modeling the curved geometry of the interface within the cohesive zone, thereby compromising the effectiveness of the method in further enlarging the CE beyond the scale of the cohesive zone. Samimi et al developed an adaptive CE where the linear CE was enriched with piecewise linear shape functions with adaptive peaks. The adaptive peaks were updated with respect to crack growth so as to better model the evolving cohesive zone within the CE domain. 17,18 The enrichment helped obtain more accurate and smoother load-displacement responses than those obtained from using standard linear CEs. However, the issue of the overprediction of the peak load with large CE was not addressed. van der Meer et al 19 employed the level-set method together the Hansbo's method (or otherwise known as the phantom node method) 20 to represent the delamination front within elements. They adopted a fracture mechanics-based approach and calculated the energy release rate for delamination based on the stress discontinuities across the interface under the assumption of classical lamination theory. Their approach effectively and efficiently overcame the cohesive zone limit, allowing elements much larger than the cohesive zone to be used. However, the approach was limited to the modeling of single delamination in the laminate. Nevertheless, their work showed that even coarse-mesh models could predict delamination accurately when equipped the right driving force. Guiamatsia et al [21][22][23] employed the analytical displacement solutions of a beam on elastic foundation as enrichment functions to simulate two-dimensional (2D) delamination, which enabled the use of CEs much larger than the cohesive zone.
However, the computational cost incurred by using enrichment could offset the benefit of using larger elements. 23 In addition, the extension of their method to three-dimensional (3D) were not evident and, to the authors' knowledge, are not yet available. Nevertheless, their work showed the importance of correctly reproducing beam theory solutions in overcoming the cohesive zone limit in delamination. Alvarez et al 11 used quadratic shape functions for the substrates and CEs, and increased the number of integration points in the CEs to a maximum of 30. They showed that quadratic shape functions modeled the kinematics of delamination much more efficiently than linear ones, and that increasing the integration points had a clear positive influence on the accuracy of the predicted load-displacement curves. The authors claimed that using quadratic CE achieved mesh-independent modeling of Mode I cohesive fracture. However, the node-to-node distance in the coarsest mesh was still comparable to the cohesive zone length. Latifi et al 24 developed a new interface element formulation based on the thick level-set method. 25 The cohesive zone length was essentially replaced with a user-defined length parameter, the so-called damage zone length, which defined the thickness of the level-set boundary. This gave more room for the user to relax the mesh constraint. However, the full potential of this approach was not investigated, as the eventual element size used for the delamination simulation was still less than 1 mm, comparable to the cohesive zone length. Cohesive cracks have also been modeled within the framework of Isogeometric Analysis to benefit from its higher solution smoothness. [26][27][28][29][30][31] However, the specific issue of overcoming the cohesive zone limit was not resolved as the coarsest mesh used was still conforming to the cohesive zone limit. Främby et al 32 developed an equivalent single-layer shell approach which adaptively refined itself through the thickness to model delamination, thereby reducing the DoFs needed in the thickness direction. However, the mesh in the plane of delamination still needed to satisfy the cohesive zone limit. Lu et al 33 employed the floating node method (FNM) 34 to design adaptive substrate and CEs to model delamination. Extra internal DoFs were used to adaptively refine and coarsen the elements containing the cohesive zone. Since refined elements were only used where cohesive zones were present, their method allowed the use of an overall coarse mesh and significantly reduced the CPU time needed to simulate delamination with CEs. However, their approach would become much more complex in 3D, especially if both matrix cracks and delamination were to be modeled as in the context of FNM. 35,36 From the above review, one can see that there is still a need for a solution to the problem of the cohesive zone limit in modeling delamination using CEs. We hereby propose a new method which overcomes the cohesive zone limit. In this article, we demonstrate the method in 2D for modeling delamination between slender beams. However, the method can be extended to model delamination between thin plates and shells. The article is organized as the following. Section 2 presents the idea and the method in detail. Section 3 demonstrates the effectiveness of the method on a series of applications. It investigates the influence of geometrical nonlinearity, makes a comparison of the computational costs against those of the standard CE models, and shows some parametric studies on the order and type of integration schemes. In the end, Section 4 concludes with the main findings of this work.

Overall idea
Several of the reviewed methods in Section 1 share some recurring features which shed light on the potential enabling techniques for overcoming the cohesive zone limit: • The method of Guiamatsia et al [21][22][23] enriched the CE's shape functions with Euler-Bernouilli (EB) beam solutions, while the method of van der Meer et al 19 calculated the delamination driving force based on stresses obtained from classical lamination theory. Although the two methods appear to be completely different, they share the same fundamental kinematic assumption, that is, the Kirchhoff-Love hypothesis. This kinematic assumption is not enforced in delamination models made of continuum and CEs or Mindlin-Reissner shell and CEs.
• In terms of numerical integration of CEs, recent work have moved beyond standard Gauss or Newton-Cotes scheme.
To reinforce the above two points, the problem of cohesive zone limit is revisited by simulating the DCB test featured in Turon et al. 5 The choice of DCB is motivated by the fact that mode I delamination exhibits a smaller cohesive zone length than the mode II and mode III cases. 6 Therefore, it is the most relevant case to demonstrate the cohesive zone limit. The geometry and boundary conditions of the DCB test are shown in Figure 2A,B, respectively. First-order plane-strain quadrilateral elements are used to model the two arms, connected by a layer of zero-thickness linear CEs for delamination. Reference results of displacements and stresses are obtained using a fine mesh of ∼0.01 mm along the expected delamination direction and five elements across the thickness of each arm. The CEs are assigned with the classical bilinear cohesive law and the standard two-point Newton-Cotes integration scheme. The material properties are listed in Table 1.
A zoomed view of the mesh near the crack tip region in the DCB model is shown in Figure 3A. The normal stress ( Y ) distribution around the crack tip is shown in Figure 3B. In Figure 3C, the normal opening and damage variable distributions along the CEs at the interface are presented. Two observations are highlighted from these results: • Figure 3B indicates that the normal stress distribution is not entirely smooth, as the peak is a kinking point which is  0 continuous only. If a CE spans across the cohesive zone, then the peak lies within the CE and the stress field will be nonsmooth at the peak. When integrating the internal energy of this CE, the integrand, being the product of the stress with the separation, will also be nonsmooth at the peak. Hence, the classical Gauss or Newton-Cotes schemes for polynomials will not be sufficient as they would assume the integrand to be smooth up to a corresponding order, which is not entirely the case here. Therefore, a more refined integration scheme should be used for CEs larger than the cohesive zone.
• While the normal stress is not smooth at the peak, the corresponding displacement field, that is, the normal opening, has no kink at the location of peak stress. It is in fact approximately the deflection solution of an slender (EB) beam under distributed loading, which would be  1 -continuous along the neutral axis. In 3D, this would correspond to the Kirchhoff-Love plate/shell deflection solution. The classical Lagrangian continuum finite elements do not enforce the Kirchhoff-Love hypothesis associated with slender structures and they are  0 continuous across element boundaries. Hence they are not well suited to represent accurately the deflection of a slender composite ply near the cohesive zone.
The second highlight above contradicts the finding from the work of Alvarez et al, 11 where using quadratic CE with higher-order fixed integration achieved mesh-independent simulation of Mode I cohesive fracture. However, the cases studied in Alvarez et al 11 were not representative of delamination between slender structures and that the same technique would still suffer from cohesive zone limit if it were to be applied on the DCB simulation in Figure 3. This argument is supported with results in the appendix section A1. Based on the above observations and analysis, we propose a new method for modeling composite delamination which is composed of the following main ingredients: • slender structural elements formulated based on the Kirchhoff-Love hypothesis are used to model the plies; • a compatible CE formulation is developed to model delamination between the ply elements; and • a higher-order adaptive integration scheme is developed for integrating the CEs.
In the following sections, the above ingredients will be addressed in more details.

Slender ply element
Fiber-reinforced composite plies have a typical thickness of around 0.125 mm. They are thin compared with the sizes of the structures that they are applied on, such as automobile and aerospace components. Therefore, the theories of slender structures can be safely applied to the modeling of composite plies. Such slender structural elements already exist, such as the EB beam element in 2D and the Kirchhoff-Love plate/shell element in 3D. They provide  1 -continuous approximation to the deflection solution of the slender composite plies along the neutral axis. In 2D, the EB beam element has been well established and widely implemented in commercial FE packages. Hence, the full details of its formulation will not be revisited here. Instead, only the DoFs and interpolation functions are presented for reference in the later sections. A 2D EB beam element of initial length l has two nodes. Each node has three DoFs, that is, u,v, and . They are the horizontal displacement, vertical displacement and rotation, respectively. The horizontal displacement u is interpolated through linear polynomials: where and The vertical displacement v is interpolated through cubic Hermite polynomials: where , and The rotation is interpolated through the derivatives of the above shape functions: where The above displacement interpolations will be used in the following section to compute the opening vector of the CE between two EB beam elements. Although the EB beam theory itself is formulated for small-strain applications, the EB beam element as a whole can undergo large rigid-body displacement and rotation while still having small strain in the element's local frame of reference, that is, the corotated frame of reference. Therefore, a mesh of EB beam elements could model the class of geometrically nonlinear problems with large displacement but small-strain. This is usually the class that composite delamination problems belong to.

2.3
Structural CE between slender plies

Opening vector
To model delamination between the plies, the CEs must be kinematically compatible with the ply elements, that is, they need to share the same displacement interpolations along the interface (c.f., Figure 5). In order to achieve this, a 2D CE F I G U R E 4 Nodes and DoFs of the CE in its local frame of reference

F I G U R E 5 Compatible CE between two beam elements
is constructed as shown in Figure 4. Four nodes are used, where nodes 3 and 4 are shared with the upper ply element and nodes 1 and 2 are shared with the lower ply element. Similarly as in the ply element, each node has three DoFs, u,v, and . The element has initially zero thickness. However, to ensure a better visualization of the CE geometry, it has been depicted as having a finite thickness. Note that the displacement interpolations in Equations (1) and (4) are for the neutral axis of the EB beam element. They need to be projected onto the outer surfaces to represent the displacements of the CE surfaces ( Figure 5). The Mode I opening of the CE can be expressed as: where v topCE and v botCE are the vertical displacements of the top and bottom CE surfaces, respectively. They are not the same as those of the neutral axis of the top and bottom EB beam elements. In fact, an additional term is required, which includes the additional surface displacements of the CE induced by the rotations of the beams' neutral axis and the offsets of the beam surfaces from the neutral axis. Hence, using EB beam theory and Equations (4) and (7), they can be expressed as: where h is the beam thickness. The Mode II opening of the CE can be expressed as: where u top and u bot are the horizontal displacements of the top and bottom CE surfaces, respectively. Considering again the rotations of beams' neutral axis and the offset of beam surfaces from neutral axis, they can be written as: The DoF vector of the CE can be defined as: Combining Equations (9) to (13), we have: whereΔ I andΔ II are the additional openings induced by the rotations of the neutral axis of the EB beams and the offsets of the beam surfaces from the neutral axis. The opening vector of the CE, that is, , can be expressed as: where and̃= In the geometrically linear formulation, the opening vector in Equation (16) is described in the local coordinate system of the initial interface. In a geometrically nonlinear formulation, the opening vector can be described in the corotated coordinate system of the current interface. Figure 6 illustrates the differences between the two formulations in the fracture mode decomposition of the opening vector. In this work, the corotational formulation is adopted to accurately evaluate the fracture modes of the current interface. While in the corotational formulation, if the rotations are assumed small, we could approximate the sin top∕bot and cos top∕bot terms by their first-order Taylor expansions and linearize the expression of̃(and hence also ) with respect to q CE . As the small-rotation assumption should be applied consistently across beams and their interfaces, the case of moderate rotations with nonlinear (eg, von Kármán) strains for the beams can no longer be considered. To keep the method generic, we do not perform this linearization at this stage.
The rest of the section details the steps to express the opening vector in the corotated coordinate system. First, as in the case of the geometrically linear CE, the local coordinate system ( , ) of the initial CE is used to express the vector quantities in Section 2.3.1. Γ 0 denotes the mid-surface of the initial CE domain and Γ denotes the mid-surface of the current CE domain. It is a parametric curve of the coordinate . Before deformation, the points on Γ 0 , represented by the position vector 0 (in the local coordinates ( , )), can be expressed as: where l is the initial length of the CE, which coincides with the ply beam element. After deformation, the displacement vector of the mid-surface is the middle interpolation of the displacement vectors of the top and bottom surfaces of the CE. Hence, after deformation, the points on Γ, denoted by the vector , can be expressed as:

F I G U R E 6
The components of the opening vector in the geometrically linear formulation (blue) and the corotational formulation (red) The rotation of Γ with respect to its initial configuration is the angle that Γ forms with the axis (ie, Γ in Figure 6): Here, the above expression is obtained by inspecting the geometry of the deformation. In a general case, the angle of rotation of an interface can be obtained by using the Nanson's formula, as shown in the appendix in Section B1. From Equation (20), we have: Using Equations (10) and (12), we obtain: Substituting the above equations into Equations (22) and (23), we get:

F I G U R E 7 A cohesive interface under equilibrium [Colour figure can be viewed at wileyonlinelibrary.com]
The differential dΓ (which will be used later for integration) can be evaluated as: The angle Γ is used to perform the transformation of the opening vector to the corotated coordinate system (̂,̂) to obtain̂. The transformation matrix, here referred to as , is the following: The corotational opening vector in the current configuration can finally be obtained as: Figure 7 shows the equilibrium of a cohesive interface in a loaded body. The opening is exaggerated for plotting. It is assumed that before the cohesive interface fails completely, the top and bottom surfaces practically coincide with Γ, that is, the opening at final failure is small compared with the geometrical dimensions of the body. The strong form of the equilibrium equation is:̂=p for all points on Γ,

Stiffness matrix
wherêis the internal Cauchy traction on the top CE surface, andp is the external Cauchy traction on the top CE surface applied by the neighboring material of the body. Both quantities are expressed in the corotated coordinate system: where and p are the internal and external Cauchy traction vectors on the top CE surface, respectively, expressed in the initial (reference) local coordinate system. Based on the virtual work principle, the following weak form of the equilibrium equation can be written: where the two terms on the LHS are the internal virtual work on the top and bottom CE surfaces, respectively; and those on the RHS are the external virtual work, respectively.û top andû bot are variations of u top and u bot in the corotational coordinate system, respectively. The above weak form of the equilibrium equation of the CE surfaces can be simplified as: Based on Equation (32), we have: where Plugging Equation (37) into Equation (36) and using the arbitrariness of q CE , Equation (36) can be simplified into: where f CE int and f CE ext are the internal and external force vectors of the CE, respectively. In this work, we assume that the cohesive law, which relates the cohesive traction with the opening, is defined on the corotational quantitieŝand̂: The above relationship is often expressed through the constitutive matrix D CE : where K and d are the penalty stiffness and the damage variable of the CE, respectively. The role of K is to approximate the bonding of the interface before any damage onset and to prevent excessive interpenetration of the two CE surfaces under compression. Note that the term |Δ I | ensures that d is not affecting the normal stress when the CE is under compression as the crack would be closed. d is a function of̂through the cohesive law. In this work, the classical bilinear cohesive law is adopted. The mixed-mode formula implemented is the one proposed by Turon et al. 37 Plugging Equation (41) into Equation (39), we have: Using Equation (30) to change the integration domain from Γ to Γ 0 , we obtain: The above equation is a nonlinear equation of the DoF vector q CE through B, D CE ,̂and J. It needs to be linearized to allow Newton-Raphson solvers to find the solution through iterations. The consistent tangent stiffness matrix K would be: In the above equation, the second and third terms on the RHS form the material tangent stiffness matrix. With cohesive laws having negative traction-separation slope (eg, the bilinear cohesive law), the true material tangent loses positive-definiteness as damage grows. To aid numerical convergence with Newton-Raphson methods, only the secant material tangent (ie, the third term on the RHS) is used in the element stiffness matrix. The geometric terms (ie, the first and last terms on the RHS) are omitted. Including the geometric terms should in theory improve the convergence rate. However, for the problems we have encountered so far, using the secant material stiffness in combination with the Quasi-Newton solver is already providing satisfactorily fast convergence which will be demonstrated in Section 3. Hence, the implemented stiffness matrix of the CE, K CE , is: The evaluation of f CE int and K CE requires the calculation of B and J. While J has been adequately expressed in Equation (30), B is, however, much more complicated to express. The following section will derive in detail the expression of B used for the implementation.

B matrix
Now let us derive the more detailed expressions of B and B needed for the implementation of B in Equation (38). With the expression of B and Equation (31), we can use the chain rule to obtain: where ⊗ indicates the tensor product. Denoting = ∕ and using the formula that arctan ′ (x) = 1∕(1 + x 2 ), we can write Γ ∕ q CE as: In the above expression, the only terms needing further derivations are 2 ∕ ∕ q CE , the rest are straightforward to calculate using earlier equations. With the expressions of ∕ ∕ in Equations (28) and (29), we can calculate their second derivatives with respect to q CE : In the above equations, the only terms needing further derivations are bot ∕ q CE and top ∕ q CE , the rest are straightforward to calculate. Using Equation (7), we can write: Hence, we arrive at the following expressions for bot ∕ q CE and top ∕ q CE : The final form of Equations (49) and (50) would then be: With the above equations, the expression of B in Equation (47) is fully defined. With the expression of in Equations (16) and (18), B in Equation (38) can be expanded as: All terms in the above equation are already derived in earlier equations. With both B and B well defined, the matrix B and the final matrix B in Equation (38) can be readily implemented. f CE int and K CE can then be integrated according to Equations (44) and (46), respectively. Because of the dependency of D CE on d which evolves according to a cohesive law, special care must be taken when integrating the CE, especially when d at different locations of the CE domain are on different phases of the cohesive law. Section 2.4 will explain the integration scheme proposed in this work to accurately and efficiently integrate for f CE int and K CE .

Linearized versions for the cohesive opening
For problems with small rotations, two linearized versions of the formulation can be derived: a corotational but smallversion where linearization is performed with respect to top and bot (ie, Small-version), and a classical geometrically linear version where the cohesive mid-surface is always the initial mid-surface without corotation. These two versions are presented in this section. The different versions are summarized in Table 2.
In the small-formulation, the displacements of the top and bottom surfaces become:

Nonlinear
Corotational formulation with full expression of̃.
The above equations return the following cohesive opening vector: The small-assumption implies that Γ is also small, as it is an angle in between top and bot . A consistent linearization would reduce the transformation matrix to be the following: The matrix B then becomes: Note that here we have retained the full expression of when taking its derivative with respect to Γ . Had we followed along the previous linearization and used the linearized expression of in Equation (61) for the derivation of B , we would have obtained zeros on the diagonal blocks of B , which could lead to numerical instability based on our experience. The calculation of Γ still follows Equation (21) and that of Γ ∕ q CE follows Equation (48). They require the following derivatives of the interface with respect to the parametric coordinate: The B matrix B is simply: In the Geo-Linear formulation, on top of the simplifications of the small-version, no corotation is considered for the cohesive mid-surface. Therefore, the integration is always performed in the initial domain Γ 0 . Consequently, the angle Γ is no longer used, hence Γ = 0 and B = 0. The Jacobian is a constant as there is no change of integration domain during analysis: When taken with respect to the natural coordinate between -1 and +1 for numerical integration, the Jacobian is simply:

Material frame-indifference
The principle of material frame-indifference states that the constitutive law must be invariant under a change of frame. 38 This principle is relevant to the observed space, that is, the physical space in which the body moves and deforms. 39 Under a change of frame, a spatial point x (ie, a point in the observed space) is mapped to another spatial point x * through the following relationship: where Q and c are the rotation and translation of the observed space with respect to the new frame, respectively. The symbol * marks quantities after the change of frame. A spatial vector g is objective (frame-indifferent) if it simply rotates along: g * = Q g. A second-order spatial tensor G is objective if G * = Q G Q T . Note that the change of frame does not affect the reference configuration as it is a mathematical object (defined in some reference space) associated to the body (ie, the so-called "body manifold" 38,40 ) for the purpose of labeling material points only. 38,39 Hence, the reference configuration does not need to exist physically, though it is often set to be the initial configuration of the body. Therefore, quantities defined entirely in the reference space such as its frames, material points, and vectors, the Green-Lagrange strain tensor and the second Piola-Kirchhoff stress tensor are invariant under a change of frame in the observed space. 39 To ensure that the constitutive law is also invariant, a convenient choice is to define the constitutive law directly on quantities in the reference configuration. 41,42 Corotational quantities are similar in nature to reference quantities, in the sense that they are both described in a frame of reference that is associated to the body rather than to an independent observer. Hence, intuitively, corotational quantities should also be invariant under a change of frame. In the rest of this section, we will derive that the corotational quantities in the constitutive law of this work Equation (40), that is,̂and̂, are invariant. In this work, the reference configuration is set to be the initial configuration of the body. The reference space effectively coincides with the observed space at the initial instant of time, but they are still two different spaces in the sense that the former is only mathematical and the latter is physical. Both of them are, however, in an underlying Euclidean space and their global frames coincide when there is no change of frame in the observed space. The two global frames will be referred to as "reference global frame" and "current global frame," respectively. Looking back at Figure 6, the local frame ( , ) in the reference space may not coincide with the current global frame. We can define a matrix 0 which transforms a vector from the observed space to the reference space and describes it under the coordinate system of the reference local frame (hereafter referred to as the reference local coordinate system). Since the reference local frame is defined in the reference configuration, it is invariant. If a vector is described in the reference local coordinate system, it is also invariant as its components will not change under a change of frame in the observed space. Hence, for an arbitrary objective spacial vector g (eg, a fiber in the current configuration), we can write: In the case where the reference local frame coincides with the reference global frame (ie, beam initially lying along x-axis), 0 would be an identity matrix which transforms to Q T when there is a change of frame in the observed space. The transformation of a vector from the current global frame to the corotational frame would be the compound transformation 0 . is entirely defined by the current normal vector of the interface described in the reference local coordinate system (see section B1 of the appendix). As argued earlier, the components of the normal vector will be invariant, hence is invariant. For an objective spacial vector g, we can derive: The above equation shows that the corotational transformation of an objective spacial vector is invariant under a change of frame in the observed space. With this result, we can show that̂and̂are invariant by showing that and are objective spacial vectors. The following existing result will be useful for the derivation: • the Cauchy stress tensor is objective, that is, * = Q Q T . 39 Let n be the unit normal vector at a point on the cohesive mid-surface in the current configuration, with a starting point x a and ending point x b . Under a change of frame, making use of Equation (66), we have: Therefore, n is an objective spacial vector. Based on the assumption that the upper, lower and middle cohesive surfaces practically coincide with each other, we can equate the cohesive traction with the surface traction obtained from the Cauchy's theorem: Since both and n are objective, it is straightforward to verify that is objective. The opening vector is a spacial vector between two material points in the current configuration. Using Equation (66), it is easy to show that it is objective.
Finally, based on the above results, we arrive at the conclusion that botĥand̂are invariant. Hence, the constitutive law of the form in Equation (40) is invariant, thereby satisfying the frame-indifference principle.

Higher-order adaptive integration scheme
As the stress field near the cohesive zone is not smooth, the standard Gauss or Newton-Cotes integration schemes would not suffice to provide accurate evaluation of the integral and more refined schemes are needed. Compared with the adaptive subdomain integration schemes, 7,9,17 a fixed higher-order scheme has the advantage of being easier to implement, without having to worry about the update of integration point locations and the transfer of history variables all the time. 18 However, as the number of integration points increases, so does the computational cost. In this work, we propose to establish an adaptive scheme which employs a higher-order quadrature rule only in the elements that have overlaps with the cohesive zone. The CE will essentially switch between two Gauss quadrature rules: a lower-order one and a higher-order one. The lower-order scheme is used for intact and completely failed elements. Note that with the beam shape functions being cubic, the integrand within the expression of the internal force vector in Equation (39) could be of sixth order even before damage onset. The 4-point Gauss scheme, accurate up to the seventh order, is the minimum to sufficiently integrates the internal force vector when the element is intact. In practice, the 3-point Gauss scheme, obviously more economic, seems to work well in many cases (c.f., Section 3.1.5). The 4-point scheme is, however, set as the default lower-order scheme as its seventh-order accuracy ensures a wider range of numerical stability in adaptive situations (as the change of integration schemes in neighbor elements could induce local variations of nodal forces) across all the cases studied in this work. The choice of 30 Gauss points as the default for the high-order integration ensures sufficient coverage of the stress variation within the CE domain when a cohesive zone is passing through, particularly in large elements of 5 mm in length. The new scheme, although being adaptive, does not have to transfer history variables when changing the quadrature rule. This will be explained in the following paragraphs. The flowchart in Figure 8 illustrates the whole integration process. Every CE is attributed with two status flags: an iterating status which represents the CE status in the last iteration (which may or may not converge) and a converged status which represents its last converged status, that is, its status at the end of the last converged iteration. They can be either "Intact," "Damaged," or "Failed." The status is decided based on the index of the damage initiation criterion and the values of the damage variable d at the integration points of the CE. They will decide the order of quadrature rule that the CE will use.
Initially, the default converged and iterating status of all CEs are "Intact." During the analysis, the damage initiation criterion is probed at 30 integration points. Note that only the initiation is probed at this stage, not the value of d. At this step, two scenarios could occur: 1. If damage does not initiate in any of the 30 integration points, the CE remains intact and its subsequent integration will be switched to a 4-point quadrature rule. Note that no transfer of history variable across the integration points is needed here as the CE is still intact and the values of d would all be 0. 2. If damage initiates in any of the 30 integration points, the iterating CE status is updated to be "Damaged" and the subsequent numerical integration will be switched to the 30-point rule where d at all the 30 integration points will be calculated for the current iteration.

F I G U R E 8 Flowchart of the adaptive integration scheme
Once the values of d at all the 30 integration points are obtained, a check is performed to see if all of them have reached the value of 1. At this step, the two possible outcomes are: 1. If not all of the 30 values of d have reached 1, this indicates that some of the points are still going through the cohesive damage process and are not yet fully failed. In this case, the iterating status of the CE remains to be "Damaged." 2. If all of the 30 values of d have reached 1, this means that the entire CE has fully failed. In this case, the iterating status of the CE will be updated to be "Failed." After the numerical integration of the element stiffness matrices and internal force vectors, the solver assembles them into the system stiffness matrix and internal force vector, and decide if the current iteration converges. In case of convergence, a new loading increment will be applied if the total loading is not yet finished; otherwise the analysis would be completed. In case of no convergence, the solver estimates the new trial solution based on the residual and the stiffness matrix and start the next iteration. Therefore, the integration in the current iteration of the CE will be affected by whether the last iteration converges or not.

F I G U R E 9 Integration points (marked with red dots) in the adaptive integration scheme
Assuming that in the last iteration, the iterating CE status is updated to "Damaged," that is, damage onset has been detected within the 30 Gauss points. In the current iteration, the integration of this CE will then have multiple scenarios: • if the last iteration is converged, then the converged CE status will be updated to be the iterating CE status, that is, converged CE status = "Damaged." The integration will directly go to the 30-point scheme; • if the last iteration is not converged, then the converged CE status remains intact. But since the iterating CE status is damaged, the integration will also go to the 30-point scheme according to the flowchart in Figure 8. The reason for this construct is that even if the last iteration did not converge, the fact that the last probing action detected damage onset means that damage onset is happening or close to happen in this CE, hence it should switch to the 30-point scheme for higher accuracy.
Therefore, as soon as damage onset is detected in one iteration, the integration of the CE in the next iteration will go directly to the 30-point rule and d at all the 30 integration points will be updated.
Once the converged status of the CE becomes "Failed," the integration will be switched back to a 4-point rule, as the cohesive zone has left the domain of the CE and the stress profile within the CE would be relatively simple. The only nonzero stress would be the normal stress due to contact between the CE surfaces. Note that at this step no transfer of history variable is needed in practice, as the values of d at the 4 quadrature points will obviously be 1.
The above adaptive integration scheme ensures that only a limited number of CEs need to be integrated on a high-order quadrature rule. The intent here is to use a high number of quadrature points only in the regions near the cohesive zone. This is shown in Figure 9, where a zoomed view of a crack-tip region in a DCB simulation is presented. With the adaptive integration scheme, the higher-order quadrature points essentially follow the propagation of the cohesive zone.

RESULTS
The method in Section 2 is implemented in 2D within the commercial finite element software package, Abaqus. The EB beam elements of Abaqus are directly used to model the plies. The structural CE in Section 2.3 is implemented as an Abaqus user-defined element. The solution-dependent variables of the user-defined CE contain the status information of the CE, which are used to determine the corresponding integration scheme according to the algorithm presented in Section 2.4. The validity of the method will first be demonstrated on three classical test cases, namely, the DCB, the ENF and the 50%-MMB cases. Aspects such as load-displacement curves, geometrical nonlinearity, CPU time gains and parametric studies on integration will be covered. Second, applications on a buckling-induced delamination problem and a multiple-delamination problem will be presented. Finally, cohesive traction profiles and a technique of getting out-of-plane stresses for the EB beam plies will be demonstrated. All the simulations were performed with the Quasi-Newton solver and no form of damping or viscosity was used.

The standard set: DCB, ENF, and MMB
The details of the DCB case has already been introduced in Section 2.1. The ENF and the 50%-MMB cases are drawn from Turon et al. 43 These two test specimens have the same geometry (c.f., Figure 10C) and are made of the same material (c.f., Table 3). The boundary conditions for ENF are shown in Figure 10A and those for MMB in Figure 10B.  The measured cohesive zone length in the benchmark tests are 0.84, 4.8, and 1.8 mm for DCB, ENF, and MMB, respectively. The nonlinear formulation of the CE has been used to model the interface delamination. To investigate the effect of bending-stretching coupling in the beams, both the linear beam elements and the nonlinear beam elements have been used to model the plies.

Load-displacement curves
The load-displacement curves of the simulations with the linear beam elements are shown in Figure 11, along with the analytical solutions based on modified beam theory 43  In the DCB case, the standard model requires an element size of 0.2 mm in order to reach a good agreement with the analytical curve. With 1 mm and larger elements, the results of the standard model are far off from the theoretical solution. The predicted peak load at which delamination initiates is ∼185% of the analytical or experimental value. With the structural CE method proposed here, the prediction by 2.5 mm elements is significantly improved, with peak load predicted at ∼104% of the analytical or experimental value. The slight overprediction of the initial stiffness is due to the lack of transverse shear deformation in the EB beam elements, which is, however, accounted for in the analytical DCB solution based on the modified beam theory. 6 In the case of even larger CEs of 5 mm in length, the standard model predicts no delamination to occur within the duration of the simulation, while the structural CE method still delivers close agreement with the experimental curve, with peak load predicted at 108% of the experimental value. After the peak load, the proposed method's predictions with 2.5 and 5 mm elements still oscillate closely around the analytical and experimental curves. In the ENF case, the standard model with 2.5 mm linear CEs would predict accurately up to the peak load, but would have a slight overprediction of the postpeak load-displacement curve of stable crack propagation. The standard ENF simulation using 5 mm linear CEs would overpredict the initial stiffness and fail to capture the load-drop completely. These results are consistent with the measured numerical cohesive zone length of 4.8 mm from the fine-mesh benchmark simulation. The predictions of the structural CE method, with both 2.5 and 5 mm meshes, are in excellent agreement with the analytical curve. The standard MMB simulation using 2.5 mm linear CEs does not capture the load-drop. With the structural CE model, the prediction of the 2.5 mm mesh follows closely the analytical curve. The 5 mm mesh still predicts well the peak load, but suffers some overprediction on the softening part of the curve. Overall, the results demonstrate the capability of the proposed structural CE method in predicting accurately the load-displacement responses of the different delamination tests with large elements of up to 5 mm in length.

Influence of using geometrically nonlinear beam formulation for the plies
In this section, the beam elements for the plies are switched between geometrically linear and nonlinear formulations to investigate its effect on the standard set of DCB, ENF, and MMB results, while keeping the CE formulation fixed to the nonlinear formulation. We found that the DCB and ENF results are practically unaffected by the switch. The MMB results, however, show a small but clearly visible difference. The softening curves with nonlinear beam formulation shift towards the bottom-left with respect to their linear counterparts, as shown in Figure 12A. We suspect that this shift, which represents a lower mixed-mode fracture toughness than analytically predicted by linear elastic fracture mechanics, is due to the bending-stretching coupling in the top beam. In MMB, the top beam bends more than the bottom beam, particularly after the damage onset of the interface. If we consider the effect of bending-stretching coupling, then the top beam would contract with respect to the bottom beam due to its relatively larger bending, and this contributes to lower the percentage of Mode II opening with respect to the mid-surface of the CE, as shown in Figure 12B. With the material parameters here, this effect is probably small, but MMB simulation is also known to be much more sensitive to parameter variations than DCB and ENF. 43 Further investigation is, however, needed to confirm this hypothesis.

Performance of linearized formulations
The performances of the two linearized CE formulations (Small-and Geo-Linear formulations in Table 2) are also tested on this set of problems with a 2.5 mm mesh and a 30-point Gauss integration scheme. As the two linearized versions assume small rotations, they are used in combination with geometrically linear beam elements for the plies. It can be seen in Figure 13 that the linearized formulations predict the correct load-displacement curves in all cases. Hence, for these problems, the nonlinear formulation is not necessary and the two linearized versions are good enough. In a general case where nonlinear beam elements are used for nonnegligible rotations, the small-or geometrically linear simplifications cannot be made and the nonlinear formulation should be used for the CEs.

Computational time gain
Although being able to enlarge the size of CE beyond the cohesive zone limit is important, the ultimate goal of the new method is to reduce the computational time of delamination simulations using CEs. In order to test the time gain, the  results of the proposed method on 2.5 mm meshes are compared against those of the standard model on 0.2 mm meshes, on the basis that both sets of results are in close agreement with the analytical curves. The same solver settings 2 have been used across the standard models and structural CE models. The comparisons, summarized in Table 4, show that the structural CE method is faster in all three cases by at least one order of magnitude. If the curves of the structural CE models on 5-mm meshes can be accepted, then the structural CE 5-mm models would achieve an additional 33% speed up over their 2.5-mm counterparts, thereby reducing the simulation time of this set of problems to between 9 to 11 CPU seconds.

Parametric studies on integration
In this section, the differences between different integration schemes with different numbers of integration points are investigated. The 5 mm CE mesh is chosen in this investigation because the differences on a coarse mesh are more visible. In Figure 14, the load-displacement curves are reported for 3, 4, and 10 integration points. Here, the number of integration points is not made adaptive during analysis to keep the comparisons clear and focused. Overall speaking, using more integration points lead to better agreement with the analytical solutions. With 3 and 4 integration points, Gauss scheme gives much better predictions than the Newton-Cotes scheme. As the number of integration points increase to 10, the two

F I G U R E 13
Load-displacement predictions of the standard set by the linearized formulations listed in Table 2 [Colour figure can be viewed at wileyonlinelibrary.com] schemes converge to similar curves. The ENF simulations are, however, practically unaffected by the change of integration schemes, as shown in Figure 14C,D.
In order to obtain a more quantitative comparison between the different integration schemes, the following normalized L2 error is used to quantify the differences between numerical load-displacement curve f * (x) and analytical load-displacement solution f (x) in Figure 14: Since the differences between the analytical and numerical curves are mostly happening after the load-drop point, the location of the peak force in the analytical curve (the red dot on the x-axis in Figure 15) is set as the left limit of the integral domain (a in Equation (72)). The right limits of integration are 12, 8, and 8 mm for the DCB, ENF, and MMB cases, respectively. As the integrand in Equation (72) is quite complicated, the trapezoidal integration scheme is used to evaluate Equation (72). The integration domain is equally subdivided into segments and on each segment a trapezoid is formed to approximate the integral. The number of trapezoidal subdivisions is increased until the summed integral converges to a fixed value. The results have been reported in Figure 16. In DCB and MMB, clear differences can be observed between the two integration schemes and also between the different numbers of integration points. In ENF, with the used shear strength of 58.9 MPa, very low error is obtained as compared with DCB and MMB in all cases. The choice of integration scheme or number of integration points used does not seem to have an effect (c.f., Figure 16B). However, the value of   Figure 16D). We believe that this is because the initial shear strength used (ie, 58.9 MPa) is in the lower range of typical shear strength values of composites, resulting in a relatively larger numerical cohesive zone which helps smooth the stress field there. Overall speaking, the results verify the difference between Newton-Cotes and Gauss schemes at lower orders and their similarity at higher orders. Gauss scheme shows less error than Newton-Cotes in general. Using odd numbers of Gauss integration points seems to give better performance than using even numbers of them in the ENF and MMB cases.

Buckling-induced delamination
The geometrically nonlinear formulation of the higher-order structural CE is assessed in a buckling-induced delamination test originally proposed by Allix and Corigliano. 44 This same test has been used by different authors as benchmark for their numerical models, such as de Borst 45 and Zhi et al. 46 The geometry and boundary conditions of the test are shown in Figure 17 and the material properties are reported in Table 5. To induce buckling, a small (6 × 10 −4 mm) vertical perturbation was initially applied, and subsequently an horizontal displacement drove the buckling-induced delamination. They are indicated in Figure 17 with P and P, respectively. In their original article, Allix and Corigliano used four different values of fracture toughness for the CE's cohesive law, ranging from 0.2 to 1.6 kJ/m 2 . The four tests have been simulated with the proposed method. The test with the smallest value of fracture toughness (0.2 kJ/m 2 ) has the smallest cohesive zone length, which was found to be 0.25 mm through fine-mesh simulations (0.01 mm elements).
The main purpose of this test is to show that the proposed method is suitable for geometrically nonlinear analysis. In addition, the attainable mesh size and CPU time are compared with the values reported in the literature. de Borst 45 reported the use of 100 elements (0.2 mm) along the length for the G c I = 0.8 kJ/m 2 case. Zhi et al 46 reported that they used a mesh of 0.05 mm elements, and that 0.3 CPU-hours were required to complete the simulation of the G c I = 1.6 kJ/m 2 case. With the proposed method, using 0.75 mm CE for the G c I = 0.2 kJ/m 2 case still gives a solution very close to the asymptotic one, as shown in Figure 18. The G c I = 1.6 kJ/m 2 case, performed with 0.25 mm CEs, requires 51 CPU-seconds to be completed.

Double delamination
To demonstrate the capability of the proposed method in handling multiple delaminations, the problem with two initial cracks developed by Robinson et al 47 is modeled. The geometry and boundary conditions of the problem are depicted in Figure 19. The thicknesses t 1 , t 2 , and t 3 are 1.59, 0.265, and 1.325 mm, respectively. The problem has been previously modeled by Alfano and Crisfield using CEs under geometrically linear assumptions, 48 hence initially the same material properties are used here (HTA913), as shown in Table 6, and the geometrically linear formulation is used for the beam elements. The result of the geometrically linear analysis ("Linear beam" curve in Figure 20B) is consistent with the result of Alfano and Crisfield. 48 In Figure 20A, the evolution of the deformation is illustrated. In Figure 20B, a comparison is made between the simulated load-displacement curve and the experimental load-displacement data extracted from the work of Robinson et al. 47 For the simulation, elements of 1 mm in length are used with a 10-point Gauss  integration scheme. It can be seen that the predicted curve from the geometrically linear analysis follows accurately the experimental data. As the interface uses the fully nonlinear formulation of the structure CE, in principle the correct result should be reproducible when the geometrically nonlinear beam elements are used for the plies. However, it is observed that once the geometrically nonlinear formulation is activated for the beam elements, the experimental data could no longer be reproduced: the delamination always continued to propagate along the upper interface. The obtained curve is labeled "Nonlinear beam" in Figure 20B. Davies and Ankerson 49 reported in their study that this apparently incorrect curve was caused by using excessively reduced strengths for the interface for the sake of artificially increasing the numerical cohesive zone length, which is the approach taken in Alfano and Crisfield. 48 We found no report of the strength data in the original article of Robinson et al, 47 but the strength values of 3.3 and 7.0 MPa for the interfacial normal and shear strengths of HTA913 composites, respectively, are obviously artificially reduced by roughly one order of magnitude. Hence, we performed a new set of simulations using 33 and 70 MPa for the normal and shear strengths, respectively, and the results are shown in Figure 21. Now both curves (with or without geometrically nonlinear option for the beams) predict correctly the experimental data on two different meshes. This shows the danger of using excessively reduced strengths-the downside of this practice could be masked when used in combination with other modeling simplifications, such as the case here with the geometrically linear beam formulation which ignores the bending-stretching coupling. As the arms are loaded and the precrack on the top interface opens up, the very thin middle layer undergoes nonnegligible rotation (c.f., Figure 20A). This middle layer is responsible for the propagation of the top delamination: as it is pulled, its top-right-end peels off from the top arm to propagate the top delamination. If bending-stretching coupling is considered, the end loadings on the middle layer would induce combined bending and axial tension, which gives a much stiffer reaction than in the geometrically linear case. This stiffer middle layer behavior, coupled with excessively reduced interfacial strengths, helps the top-right-end of the middle beam layer peel off from the top arm more easily, thereby favoring the propagation of the delamination on the top interface and altering the whole solution.

F I G U R E 23
Longitudinal and transverse shear stress distributions in the plies near the crack tip, extracted from the reference DCB model. The X coordinate measures the distance from the crack tip, while Y measures the distance from the crack plane as in Figure 3 [Colour figure can be viewed at wileyonlinelibrary.com]

3.4
Cohesive traction and ply stresses

Cohesive traction
The traction values computed at the integration points of the CEs are stored for postprocessing. The normal traction distributions near the cohesive zone are plotted for the DCB simulations. The results are shown in Figure 22 for 1, 2.5, and 5 mm CEs. The predicted traction distribution in the cohesive zone seems to be unaffected by the change of mesh. The length of the cohesive zone is predicted to be around 2.7 mm, larger than that predicted by the reference linear CEs. The traction profile in the intact region adjacent to the cohesive zone shows that the normal stress first quickly descends into compression, then rises up again into tension. Compared with the reference traction solution obtained with the finely meshed standard model, the compression magnitude is much higher. This compression-tension oscillation also grows larger with coarser meshes. This could be due to the lack of transverse shear deformation in the EB beam elements. Figure 23 shows the longitudinal and shear stresses in the plies near the crack tip. The stresses are extracted from the integration points of the reference solid-element DCB model. It can be seen that although the transverse shear stress is small compared with the longitudinal stress, it is in the same order of magnitude as that of the normal stress in Figure 3B. Hence, neglecting the transverse shear strain in the plies (as done in the EB beam elements) will lead to overstiff behavior of the plies and this could be the reason of this high compression-tension oscillation ahead of the cohesive zone.

Ply stresses
For the analysis which requires stresses in the plies (such as intralaminar damage analysis), the stress output of the EB beam element is limited to the axial stress. However, the CEs on the interfaces have out-of-plane normal and shear tractions which make up the other components of the stress tensor in 2D. A typical laminate would have multiple plies and interfaces through the thickness, hence it is possible to interpolate the out-of-plane normal and shear stresses between the CEs to approximate the intralaminar out-of-plane normal and shear stresses. To illustrate this idea, the DCB and ENF simulations are modeled with three EB beam elements per arm with two CEs in between the EB elements. The out-of-plane stresses are then linearly interpolated from the CEs across the thickness of each ply (for the surface plies the stresses are interpolated to zero at the outer surface). In Figure 24, the interpolated stresses through the thickness at the crack tip location before damage onset are compared with the integration-point stresses from the reference solid-element models (meshed with 5 plane-strain elements across the thickness of each arm). It can be seen that the two curves differ the most near the delamination plane: the reference curve shows a sharp drop of stress in the solid elements adjacent to the delamination plane. This could be due to the limitation of the linear solid elements in capturing the stress concentration near the crack plane, resulting in an underprediction of stresses. The two solutions actually match reasonably closely on the rest part of the curve in both cases. Hence, with the proposed method, the out-of-plane stresses in plies can be well estimated through interpolation between CEs, making the method viable for intralaminar analysis. However, special care must be taken when interpolating the stress ahead of the cohesive zone due to the presence of high compression-tension oscillations as discussed in the previous paragraph. Further study is needed to come up with ways to regulate this high compression-tension oscillation to allow for reliable stress interpolations from CEs to the plies in the region ahead of the cohesive zone. In addition, it is difficult to judge if this technique would work well in a general case with an arbitrary layup. This will be investigated in a follow-up study where the plate version of the method will be implemented. Alternatively, methods on the recovery of out-of-plane stresses from beams and shells have been well studied in the literature. A stress-recovery method based on the a-posteriori integration of the equilibrium equations could be employed to equip the proposed method with the full 3D stresses for intralaminar analysis. 50-52

SUMMARY AND CONCLUSIONS
A new method of modeling delamination in composites is proposed in this work. It makes use of slender structural elements for the plies and compatible structural CEs for the interfaces. A corotational formulation is developed to allow for geometrically nonlinear analysis. It is shown that a cohesive law defined on corotational traction and separation vectors are invariant under a change of frame, hence satisfying the frame-indifference principle. A higher-order adaptive integration scheme is proposed which integrates accurately the internal force vector with the presence of the cohesive zone, and at the same time preserves the computational efficiency of lower-order integration schemes when there is no cohesive zone present. The proposed method has been implemented in 2D and validated on DCB, ENF, MMB, buckling-driven delamination and double-delamination problems. Linearized versions of the method have also been developed and validated on the DCB, ENF, and MMB problems. They also show excellent agreement with the reference solutions. Including the geometrical nonlinear effect (ie, bending-stretching coupling) in the plies is found to influence slightly the result of MMB, where it appears to lower the mixed-mode fracture toughness possibly due to its influence on the mode decomposition of the opening vector. The double-delamination problem is found to be sensitive to the cohesive strength values when bending-stretching coupling is included in the simulation. Physically correct strength values should be used. Compared with the reference method where the plies are modeled with continuum elements and the interfaces with linear CEs, it has been shown that the proposed method does not suffer from the cohesive zone limit. It gives accurate strength and postdamage predictions even with elements ten times larger than the limit size for standard CEs. As a result, the new method has achieved a computational speed-up of at least one order of magnitude over the reference method.
The Gauss and Newton-Cotes integration schemes are compared on different orders of integration. It is observed that Gauss scheme shows superior accuracy on lower orders, while both schemes show similar accuracy on higher orders. Furthermore, the simulations of the proposed method do not require the use of any numerical damping or viscosity for convergence. This showcases the superior numerical stability of this method in addition to its accuracy and efficiency.
While the EB beams only provide axial stress for the plies, interpolation can be used to obtain out-of-plane stresses from the CEs. The interpolated stresses are comparable to stresses at the integration points of the reference models. This opens up the possibility of the method to be applied on intralaminar problems. However, there is some mesh-dependent stress oscillations ahead of the cohesive zone which is probably due to the lack of transverse shear deformation in the EB beam elements. These stress oscillations need to be considered or regulated when interpolating for out-of-plane stresses in the plies. Moreover, the effectiveness of the interpolation technique in a general laminate with arbitrary layups still need to be assessed. Apart from using this interpolation technique, it is also possible to combine the proposed method with an a posteriori stress recovery method to obtain the out-of-plane stresses. This could be investigated in a future work.
The future work will be naturally focused on the extension of the proposed method to the modeling of delamination between thin plates and shells using Kirchhoff-Love plate and shell elements and the corresponding structural CEs. Once established, it is expected to overcome the cohesive zone limit of element size for the modeling of delamination in general thin-walled structures. This would bring along a step reduction in the computational cost associated with the current CE-based models, thereby widening the possibility of performing direct delamination simulation on models of large-scale thin-walled structures.

APPENDIX B. DERIVING THE ROTATION ANGLE USING THE NANSON'S FORMULA
The Nanson's formula relates a differential surface in the reference configuration, n 0 dA 0 , to itself in the current configuration, n dA, through the following equation: where n 0 and n are the normal vectors of the surface in the reference and current configurations, respectively. F is the deformation gradient. n contains the directional sins and cosines that can be used to extract the rotation angle and form the rotation matrix. When applied to the mid-surface of the CE, we have: Combining the above two equations, we obtain: Since we also have: therefore, we obtain:

APPENDIX C. ENF WITHOUT SNAP-THROUGH
In this section, an ENF test without snap-back is simulated with the proposed method to see if it captures the full load-displacement curve. The boundary conditions are the same as the ones depicted in Figure 10A. However, the geometry and the material properties are slightly different. Compared with the ENF geometry in Figure 10C, this ENF specimen is shorter overall (100.8 mm in length) and has a relatively longer precrack (30 mm in length) compared with the length of the specimen. This geometry gives a stable crack propagation in ENF without snap-back. The material properties are listed in Tables C1 and C2. The elements used in the simulations are 1 mm long with 30 Gauss integration points for the CEs. The load-displacement curve of the simulation agrees closely with the analytical curve, as shown in Figure C1.