Efficient modelling of delamination growth using adaptive isogeometric continuum shell elements

The computational efficiency of CAE tools for analysing failure progression in large layered composites is key. In particular, efficient approximation and solution methods for delamination modelling are crucial to meet today’s requirements on virtual development lead times. For that purpose, we present here an adaptive continuum shell element based on the isogeometric analysis framework, suitable for the modelling of arbitrary delamination growth. To achieve an efficient procedure, we utilise that, in isogeometric analysis, the continuity of the approximation field easily can be adapted via so-called knot insertion. As a result, the current continuum shell provides a basis for an accurate but also computationally efficient prediction of delamination growth in laminated composites. Results show that the adaptive modelling framework works well and that, in comparison to a fully resolved model, the adaptive approach gives significant time savings even for simple analyses where major parts of the domain exhibit delamination growth.


Introduction
To accurately predict damage growth in large, thin-walled laminated composite structures, it is required to have models that are able to capture relevant deformation mechanisms in a computationally efficient manner. Thus, the deformation kinematics as well as the potential material degradation mechanisms should be well represented in areas where damage occurs, at the same time as unnecessary detail should be avoided. When non-linear material phenomena such as intraply 1 damage development or delamination growth need to be considered in the analysis, an accurate prediction of the three-dimensional stress field in critical areas is key. Both because intraply damage is generally driven by the interaction of three-dimensional stresses in the ply, and because delamination is physically driven by interlaminar normal and shear stresses.
Traditional finite element (FE) based shell elements fall short in these cases, since the availability of only planar stress components is insufficient. Therefore, FE continuum shell elements, cf. e.g. Parisch [1], which possess three-dimensional kinematic descriptions and stress fields, are a natural alternative. However, such shell theories are commonly based on first shear deformation theory, meaning that the laminate must be represented by at least one continuum shell element per ply to obtain a good prediction of the three-dimensional stress state. Moreover, to capture delamination growth, all ply interfaces must be modelled using interface elements. In the end, this leads to a computationally, very expensive approach (both in terms of simulation time and memory requirements).
A strategy to alleviate the computational cost associated with such detailed modeling of each ply, which is currently increasing in popularity in the literature, is to make the through-thickness kinematics adaptive, e.g. by exploiting the partition-of-unity property of Lagrange polynomials as inspired by [2]. Several alternative, but similar, methods have been developed that support laminate failure analyses requiring initially only one shell element through the thickness with the common feature that (arbitrary) delamination propagation is considered only in areas where it is needed [3][4][5][6][7]. Common challenges associated with many of these methods are to correctly predict the onset of delamination growth since transverse stress predictions in general are of low accuracy in traditional (inexpensive) shell elements, and how to enable the kinematic enrichment in an efficient manner.
In this respect, continuum shell elements based on the concept of isogeometric analysis (IGA) provide an interesting option. In IGA [8,9], rather than using Lagrange polynomials as in traditional FE, higher-order splines (NURBS or T-splines) are used as basis functions in the approximations. A benefit of the use of such splines is that the displacement approximation is higher-order continuous, thereby providing a basis for well-performing shell [10][11][12] and continuum shell formulations [13][14][15], where in the latter case an accurate geometric description of the shell mid-surface (via in-plane NURBS) is combined with the three-dimensional stress representation of conventional continuum shell elements, realised via a higher-order through-thickness B-spline approximation.
Furthermore, as shown by Hosseini et al. [15], by adopting an isogeometric continuum shell modelling framework it is rather straightforward to, by so-called knot-insertion, modify the through-thickness kinematics to incorporate weak and strong discontinuities, see also [16]. By introducing weak discontinuities at ply interfaces, the through-thickness strain discontinuities at these locations are explicitly accounted for. This enables a much better 3D strain and stress prediction, which is key for a good estimation of the amount of intraply damage. By introducing strong discontinuities, the element is also capable to represent initiation and growth of one or several delamination cracks.
In the current paper, we extend the shell formulation from [15] into an adaptive continuum shell that allows for an update of the through-thickness kinematics at any required time instant during the simulation. The adaptivity is facilitated by that the knot insertion can be fully automated due to the hierarchical nature of the isogeometric approximation functions.
In the current paper we also demonstrate that the higher order in-plane continuity obtained with an IGA approach allows for an element-local recovery procedure for accurate prediction of out-of-plane (transverse) stresses, similar to [17,18]. Thereby, the possibility for delamination initiation can be assessed element-wise, something which is not generally possible in traditional FE based approaches if a single shell element (through the thickness) is used to represent the laminate. 2 As a result, the current shell provides a good basis for an accurate but also computationally efficient prediction of the progressive failure in laminates, without a-priory knowledge of where damage will occur. Results show that the adaptive modelling framework works well, both to predict the full 3D stress states in multiaxial laminates, but also to capture growth of delaminations. Furthermore, in comparison to a fully resolved model, the adaptive approach gives significant time savings even for simple analyses where significant parts of the domain exhibit delamination growth. This implies that computational efforts (time and memory) can be reduced considerably when using this adaptive concept in large-scale analyses where damage develop only in a confined, but initially unknown area of the structure. This paper is ordered as follows. In the next section, we will give a concise review of the isogeometric continuum shell element. Since out-of-plane stresses in the standard element are of a rather poor quality, an algorithm that reconstructs these components based on element-wise post-processing of the smooth stress fields [6], is presented here as well. Subsequently, the accuracy of the computed stresses in the element with different discretisations throughthe-thickness of the shell is assessed in Sect. 3. The validity of the reconstructed stresses is also presented in this section. In Sect. 4, a strategy to perform the automatic knot-insertion to enhance the interpolation in thickness direction is presented. The performance of the automatically adapting element is studied by means of three benchmark problems in Sect. 5. Finally, the paper is closed with some conclusions and an outlook to future developments.

Isogeometric continuum shell element
The kinematic relations and the discretisation of the isogeometric continuum shell element discussed in this section, are presented in detail in Hosseini et al. [15]. Figure 1 shows the undeformed and the deformed configuration of the continuum shell element. The reference surface of the shell is denoted by S 0 and S for each respective configuration. The variables ξ and η are the local curvilinear coordinates in the two independent in-plane directions and ζ is the local curvilinear coordinate in the thickness direction of the shell. Furthermore, the position of a material point within the shell body in the undeformed configuration is written as a function of the three curvilinear coordinates: where X 0 (ξ, η) is the projection of the point on the reference mid-surface of the shell and D(ξ, η) is the thickness director perpendicular to the mid-surface S 0 in this point. In any material point, a local reference triad can be established. The covariant base vectors are then obtained as the partial derivatives of the position vectors with respect to the curvilinear coordinates i = [ξ, η, ζ ]. We define a set of basis vectors on the reference surface in the undeformed configuration as:

Kinematics and equilibrium equations
where t is the thickness of the shell. Now, using Eq. (1), the covariant triad for any point within the undeformed shell body is obtained as: where the comma subscript • ,α denotes partial differentiation with respect to the directions in the undeformed configuration.
With this relation, we can define the covariant triad in the deformed shell element g i at any material point as: This reference triad is used to construct the Green-Lagrange strain tensor γ : where I and F respectively are the unit tensor and deformation gradient tensor. The latter one is given by: The virtual internal work δW int in a Total Lagrangian setting (referring to the undeformed configuration 0 ) is formulated as: where S denotes the Second Piola-Kirchhoff stress tensor that arises from a linear constitutive relation between the material stiffness tensor C and the Green-Lagrange strain tensor: The strain field in Eq. (9) is expressed in the framework of the contravariant base vectors G i , i = 1, 2, 3. These where, for a fibre-reinforced orthotropic material, T 1 is the fibre direction and T 2 and T 3 respectively are the in-plane and out-of-plane perpendicular directions. The virtual external work δW ext is defined as: In this expression the contribution is split in first Piola-Kirchhoff tractions t acting on the undeformed outer surface t and body forces q, see Fig. 2. The principle of virtual work states that the work done by internal forces equates the work done by external forces: This yields a system of non-linear equations that can be solved by an incremental-iterative procedure, wherein the total loading is applied in a finite number of load steps. For more information on this, the reader is referred to [15].

Discretisation
The mid-surface of the shell is constructed using NURBS basis functions as shown in Fig. 3. In the thickness direction, the displacements are discretised using higher order B-spline basis functions [14]. A B-spline volume patch is generated by multiplying the bivariate (planar) splines S i with the univariate spline function in thickness direction, H j .
The total displacement field u can then be approximated as: where a I is the vector containing the displacement degrees of freedom and N cp = n × m is the total number of control points (or shape functions), with n being the number of planar shape functions and m the amount of shape functions present along the shell thickness. The combined shape functions N I are now given by: Consequently, the component k of the displacement field can be obtained by: where the index k indicates either the x-, y-or z-directions. The displacement field u in Eq. (13) can be of any arbitrary order, but is chosen higher than one to obtain a B-spline function (for lower orders, B-splines depreciate to the traditional finite element shape functions). This means that the strain field varies at least linearly over the thickness, which is important to avoid thickness locking. This is similar to the standard continuum shell formulation, where an internal stretch term is added to obtain a quadratic term in the displacement field through-the-thickness of the shell [1].

Fundamentals of B-splines
B-spline basis functions are defined along a knot vector , which contains knots ξ i symbolizing coordinates in the parameter domain: In this definition, n denotes the number of basis functions that get defined and p represents the degree of the functions. Univariate B-splines are generated via the Cox-de Boor recursive algorithm [8], where basis functions of a certain degree are are expressed as functions of one order lower until a piece-wise constant ( p = 0) is obtained: The piece-wise constant is given by: It should be noted that each constant function defines a single element along the parameter domain (between ξ i and ξ i+1 ), analogous to an element in conventional finite element methods (FEM). B-splines of degree p are C p−k continuous at a knot with multiplicity k [19]. This means that we are able to control the continuity of the basis functions at a knot by arbitrarily selecting its multiplicity. We benefit from this property in two ways. On one hand, we use it to compute the Bézier extraction operator, required to obtain an element-wise data structure for isogeometric analysis analogous to traditional finite elements. On the other hand, this property is used to model layered structures with C 0 continuity (weak discontinuities) between the layers [16] as well as traction-free cracks or cohesive interfaces by enforcing C −1 continuity at the interfaces (strong discontinuities).

Bézier extraction
A characteristic of B-splines or NURBS of order p is that they are present along a knot interval ξ i , ξ i+ p+1 . It indicates that basis functions exist along multiple elements and thus, elements support different shape functions (see shape functions along the two elements at the left side of Fig. 4). This property is in contrast to the traditional FEM, where identical shape functions are defined within all elements. In FEM software, integral equations are typically solved on a single parent element using quadrature rules. In order to retain this finite element structure for the discretisation of the mid-surface of the isogeometric shell element, we make use of the Bézier extraction technique [20].
The principle of Bézier extraction is explained in Fig. 4. It shows that identical third order, p = 3, displacement fields can be obtained using both bases (bottom). The left basis is defined along = 0, 0, 0, 0, 1 2 , 1, 1, 1, 1 , whereas the right basis is developed along = 0, 0, 0, 0, 1 2 , 1 2 , 1 2 , 1, 1, 1, 1]. The second knot vector contains two additional knots that reduce the continuity of the basis functions halfway the parameter domain to C 0 . Via this knot insertion process, the original basis is decomposed into one consisting of two C 0 continuous Bézier elements with identical (Bernstein) shape functions B(ξ ). Upon insertion of additional knots, also new shape functions and control points are defined. To ensure that the displacement field remains unaltered, the positions of the new set of control points P new must be determined by: In this expression, P old is the array with old positions and C is the total Bézier extraction operator, composed from all operators resulting from single repetitive knot insertions: where m (in this context) denotes the total number of inserted knots. The process of finding this Bézier extraction matrix C is explained in detail in [20]. The extraction operator can then be used to map the Bernstein basis functions B(ξ ) onto the original shape functions N(ξ ): By defining N e = L e N and B =L e B e to extract only the shape functions belonging to element e, Eq. (21) can be rewritten into an element-wise mapping: This final expression allows us to use Bézier extraction to map a regular Bernstein basis B e onto the element specific shape functions N e . By doing so, the Bernstein basis functions (Bézier elements) can be used as the finite element representation of splines and NURBS. The Bézier extraction technique also holds for multivariate NURBS (surfaces or solids), which in the current work is exploited to discretise the mid-plane of the continuum element. Figure 5 shows the steps in order to create discontinuities through-the-thickness of the shell structure. Assume that quadratic B-spline basis functions H j (ζ ), defined along a knot vector = [−1, −1, −1, 1, 1, 1], have been used as the univariate functions in thickness direction. Three basis functions are then present along the total thickness of the element, which in this configuration will be called the lumped element.

Introducing discontinuities in the displacement field
Now suppose that we want to have a composite shell consisting of two layers of equal thickness. The deformation of layered composite structures requires a unique displacement at the interfaces and different strain fields in the adjacent layers. In the example of Fig. 5, this is simply achieved by having a displacement field which is C 0 continuous at the interface at ζ = 0. This leads to a new knot vector = [−1, −1, −1, 0, 0, 1, 1, 1], with five basis functions through-the-thickness. Henceforth, we will denote this element as the layered element. Subsequently, complete separation of the layers is obtained if we insert another knot, such that: = [−1, −1, −1, 0, 0, 0, 1, 1, 1], which gives six basis functions along the shell thickness. An element in this state is denoted as the discontinuous element and at the moment it obtains this configuration, a cohesive zone model is introduced to govern the delamination failure process. In Fig. 6 a discontinuous element is shown in the undeformed and deformed configuration. The cohesive crack surface that arises from splitting the element is defined in the deformed configuration as the average between the top ( + ) and bottom ( − ) crack surface.
In the figure, u + and u − are the displacements of the original material point on the top ( + ) and bottom ( − ) crack surface. The displacement jump v between two layers can be written as: The Cauchy traction vector at the interface is acquired from the cohesive zone model and is defined as: where σ is the Cauchy stress tensor and n is the vector normal to , defined below. The internal virtual work is established in the global coordinate system, yet cohesive material laws are often expressed in a local reference framework. The coordinate systems in the implementation are therefore transformed by a rotation matrix Q that relates the global directions (e x , e y , e z ) to the local ones (n, s 2 , s 3 ): where n denotes the vector normal to the surface and s 2 and s 3 indicate the shear directions. This local reference system is assumed to be the average of the respective directions on the top ( + ) and bottom ( − ) surface: The displacement jump and traction vector in the local framework are then given by: The interface behavior is accounted for when elements are in the discontinuous configuration. In that case, the virtual internal work first defined in Eq. (8) is expanded to: where the latter term is related to the cohesive tractions.

Stress enhancement scheme for lumped elements
The lumped continuum shell element yields accurate inplane stress results for layered composite structures, as will be shown in the next section. However, using this element configuration to calculate transverse out-of-plane stresses in layered composites, generally leads to a poor prediction due to oversimplified out-of-plane kinematics. To still be able to obtain reliable predictions of these stresses we adopt a strategy similar to Kant and Manjunatha [21] where improved values are recovered from the three-dimensional momentum balance equations. Here, the first and second order derivatives of the in-plane stress components need to be extracted from the solution. In traditional finite element shell models, the in-plane stress components are predicted with good accuracy, however, the stress derivatives are not well predicted because of only C 0 in-plane continuity of the shape functions.
Since the derivatives of these traditional shape functions are discontinuous across element edges, the resulting stresses are non-smooth also for elastic problems. The benefit from the IGA approach is that stresses indeed are smooth when crossing element edges, which means the derivatives of each component can be computed element-wise with good accuracy. For zero body forces under quasi-static conditions, the balance equations are given by ∇ · σ = 0, which can be written in components: We can reconstruct the transverse stress variation throughthe-thickness of the shell by rewriting Eqs. (29) into integrals over the shell thickness: where z is the local transverse direction, z (n−1) and z (n) denote the lower and upper thickness coordinate of ply n and N represents the total number of composite layers. Further, the first and second order derivatives with respect to coordinate α = [x, y] are indicated by • ,α and • ,αα respectively and the recovered values are denoted by•.
As can be seen above, the integration of the threedimensional equilibrium equations yields integration constants which in general have to be determined from the traction conditions at the top and bottom shell surface, cf. Främby et al. [22]. In that paper it is also shown that the integration constants can be used to average the integration error 3 over the thickness, such that the discrepancy between the predicted stresses at the surfaces and the applied tractions is minimised.
In our procedure to solve Eqs. (30) and (31), we make use of the continuity of the stresses and project per element the stress variation in each plane of integration points on a second order Lagrangian basis (using conventional second order finite element shape functions). This projected stress field can then, for each layer of integration points, be used to evaluate first and second derivatives of the planar stress components in the element centre point at a given position through the thickness. As a final step, we integrate the stress derivatives according to Eqs. (30) and (31) to obtain the more accurate recovered stress profiles.
Note that, for a sufficient level of accuracy of the second order stress derivatives in Eq. (31), at least a C 2 in-plane continuity is required, i.e. cubic NURBS are necessary for the in-plane displacement approximation. If a lower approximation degree is to be used, the stress projection onto the second order Lagrangian basis needs to pursued in a non-local manner, e.g. [22].

Element performance assessment
The performance of the continuum shell element and stress enhancement scheme is first studied in the simulation of the deflection of a multi-layered composite panel. For this problem, zero-thickness shell elements would generally suffice to accurately calculate the displacements, but would not be adequate to compute the fully three-dimensional stresses and strains in the individual layers accurately.
We consider the laminate shown in Fig. 7. The laminate is a = 600 mm long and b = 400 mm wide. It consists of five equally thick layers, stacked in a 0/90/0 s sequence. The laminate with a total thickness of t = 10 mm, is loaded at the top surface by a bi-harmonic pressure load; p(x, y) = −q 0 sin( π x a ) sin( π y b ). The amplitude of the pressure load at the centre is q 0 = 0.01 N/mm 2 , which leads to small deformations of the laminate. Because of symmetry, only a quarter of the model is considered.
The total model is simply supported (S.S.) along its four mid-plane edges, which corresponds to two lines of the quarter model. Along the long line, u y = u z = 0 and along the short line u x = u z = 0. At the symmetry (Symm.) surfaces, u y = 0 (long boundary) and u x = 0 (short boundary) over the complete thickness.
The analyses are all conducted using third order Bézier elements. 4 In thickness direction, fourth order B-splines have been employed, which results in a third order stress profile through-the-thickness. All numerical results in this work are mesh converged.

In-plane stress comparison
The stress results in thickness direction are compared at x = 290 mm and y = 190 mm (element centre), see Fig. 8. In Fig. 9 the stresses in fibre direction obtained by lumped and layered elements are compared to outcomes of the Classical Laminate Theory (CLT). It is evident that the results of both lumped and layered elements are in excellent agreement with the CLT solution. Similar statements can be made for the remaining in-plane stress components σ 12 and σ 22 (not included).

Out-of-plane stress comparison
The out-of-plane element stresses are compared in the global coordinate system. The analytical solutions proposed by Pagano [23] are used as a reference. Figure 10 shows that out-of-plane shear stresses are predicted very well when the elements are in the layered configuration. The resulting stress profile is symmetrical and ends up indicating zero at the top and bottom face of the shell element. According expectations, the lumped shear stress profile deviates from the reference solution. The stresses are smeared out over the entire thickness of the panel and are independent of the composite lay-up. Identical conclusions  can be drawn upon investigation of σ yz results and are therefore not shown. Figure 11 shows the through-the-thickness normal stress component along the panel thickness. Layered elements capture the physical stress profile accurately. The value of the pressure load is accepted at the top surface, while stresses vanish at the bottom surface. In the centre layer a small stress oscillation appears. This is a boundary effect caused by enforcing simply supported boundary conditions at the mid-plane [14]. The lumped stress profile does not fulfill the stress boundary conditions at the top or bottom face of the shell. Moreover, the appearing stress jumps violate the interface continuity conditions. In conclusion, elements in a lumped configuration are able to capture in-plane stress components accurately, but indeed return poor out-of-plane stress predictions.

Reconstructed stress comparison
In Fig. 12 the reconstructed out-of-plane shear stress component σ xz is presented. Here, the values of the integration constants are C 1 = C 2 = 0 as the top and bottom surfaces are traction free. From the results it is clear that the stress enhancement scheme somewhat underestimates the stress magnitudes, yet the stress profiles are in good accordance. Figure 13 shows the reconstructed result for the stress component in thickness direction of the panel. The reconstructed stress profile nearly coincides with the analytical solution. The integration constants C 3 and C 4 are determined to fit the traction boundary conditions at the top and bottom face of the elements.
The stress enhancement scheme confirms a significant improvement of the out-of-plane stress predictions for lumped elements.

Adaptive discretisation
In the previous section, different element configurations were introduced based on the discretisation in thickness direction.
In principle, the discretisation through-the-thickness of the shell element can be changed during the simulation such that elements can adopt different configurations based on an interface stress criterion. To overcome the problem of a poor out-of-plane stress resolution acquired from lumped elements, the more accurate reconstructed stresses can be used in this stress criterion when upgrading from a lumped to a layered element configuration. However, to enable such an automatic update of the discretisation, some essential problems have to be addressed. First, the initial values of the new degrees of freedom of the element need to be determined. Second, the Dirichlet boundary conditions need to be updated as well.

Initialisation of new degrees of freedom
Upon enhancement of an element, the displacement field in a control point is discretised using more degrees of freedom.
The initial values of these new degrees of freedom are not equal to zero, in contrast to most X-FEM procedures. Instead, these initial values must be calculated using the values of the existing degrees of freedom (which in turn will obtain other values too). For this purpose the Bézier extraction operator C, similar to the one formulated in Sect. 2.2.2, can be utilized (see examples in Fig. 14): In this expression, the new vector of degrees of freedom,ã i k , is related to the old vector of degrees of freedom a i k , where k indicates the x-,y-and z-directions and i denotes a planar control point.

Updating Dirichlet boundary conditions
During the upgrading of the elements, the prescribed degrees of freedom must be updated in order to maintain the same conditions at the boundaries. There are multiple ways of updating these, depending on the type of boundary condition. We will use a simply supported boundary condition to demonstrate this, see Fig. 15.
In the lumped element state, only one control point in thickness direction is constrained (at the mid-plane), such that the boundary rotation is free. During an upgrade, the original constraints are imposed on the new control point(s) located at the mid-plane. Additionally, the old constraints on the original control point(s) are removed, such that the same boundary conditions are preserved.

Extended crack tip and mixed elements
Adaptive elements adjacent to each other can be in different configurations (discretisation states). When an element upgrades, all control points that are linked to that element are enhanced. Since isogeometric elements share a larger number of control points or nodes compared to conventional finite elements, it is worthwhile investigating the inter-element compatibility.
The element configuration is determined by the control point with the least refined through-the-thickness kinematics that belongs to that element. When all control points corresponding to an element have equal discretisation states, then this is a pure element in that respective state. If this is not the case (due to neighbouring elements in a different state), then  the element is denoted as mixed and requires more consideration.
To illustrate this, we consider a propagating crack in one dimension in Fig. 16, together with the in-plane control points and shape functions. Due to the higher order continuity of the shape functions, they have support over multiple in-plane elements. In the figure, all discontinuous control points belong to element IV, yet it is clear that the crack extends all the way into element II, which is a lumped one. This originates from the fact that the brown in-plane shape function belonging to control point number four has support in element II. Since this control point is in a discontinuous configuration, the interface (out-of-plane) displacement jump also prolongs into the lumped element. The shape function vanishes exactly at the boundary between element I and II. The actual crack tip therefore is present at this element boundary, instead of in-between the layered and discontinuous element.
Due to this 'extended' crack tip, interface displacement jumps can also occur in layered and lumped elements. However, this can only occur when these are mixed elements and thus, the cohesive material law is considered also for such elements.

Numerical examples
In this section the capabilities of the element with automatically adapting through-the-thickness discretisation are assessed by means of three benchmark cases. Simulations using both adapting and non-adapting 5 elements are conducted in order to compare computational size (in terms of degrees of freedom) and run times. In all cases, the delamination behaviour is described by the cohesive zone model proposed by Xu and Needleman [24]. Simulations are compared through force-displacement characteristics and upon an identical adaptive and non-adaptive response, the improvement of the computational efficiency is quantified. In the benchmark simulations, initial cracks are modelled by initializing specific elements in the discontinuous state. For these elements, the cohesive constitutive behaviour is neglected, which yields a traction-free interface.

Double cantilever beam simulation
The first benchmark case is the double cantilever beam (DCB) shown in Fig. 17. The DCB model has a length of L = 100.0 mm, a width of b = 5.0 mm and a thickness of t = 0.5 mm. The material is considered to be elastic and isotropic, with a Young's modulus of E = 10 5 N/mm 2 and Poisson ratio of ν = 0.3. This model does not contain any initial cracks, in order to test whether the adaptive element is able to capture crack initiation. For the cohesive zone model an ultimate traction of t ult = 10 N/mm 2 has been used in combination with a mode I fracture toughness At the left boundary (x = 0), the DCB model is clamped, i.e. u x = u y = u z = 0. At the right boundary (x = L, z = −h, +h) the model is loaded by by applying positive and negative displacements in z-direction until a final value of u z,limit = ±2.0 mm is obtained. The displacement is applied in non-uniform steps as u z = λ(t) × u z,limit , where λ(t) ∈ [0, 1] is the load factor as a function of the equally sized pseudo time steps t ∈ [0, 1]. In order to describe the crack initiation with a higher resolution, the load factor is chosen as λ(t) = t 3 .
The DCB model is built up from 200 × 1 (L × b) Bézier elements with second order bivariate in-plane B-splines and univariate second order B-splines through-the-thickness of the elements. No material non-linearities other than the cohesive constitutive relation are considered.
In this example, the initiation and propagation of a delamination crack is driven by out-of-plane normal tractions at the composite interface. These stresses are used in a stress criterion to change the element configurations. The elements upgrade from lumped to layered when the normal out-ofplane traction (at any sample point at the interface) exceeds 0.05 × t ult . Similarly, the transition from layered to discontinuous occurs at a stress level of 0.1 × t ult . It is necessary to upgrade elements a certain distance in front of the fully cracked region (traction-free part in the cohesive material law) to get a correct response. On the other hand, the delami-nation driving stresses are in this case very concentrated near the crack tip, which indicates that many (small) incremental load steps are required to prevent elements from upgrading too late. Hence, it should be noted that the stress thresholds are a trade-off between accuracy and efficiency.
The numerical results are compared to the analytical prediction of linear elastic fracture mechanics (LEFM) in the damage regime: Due to the strong stress concentrations at the right boundary when the simulation starts, there is no difference in the predicted location of first failure with or without the stress reconstruction model. In any case, a mode I delamination crack initiates at the right boundary. For this reason, the stress reconstruction model is disabled for this benchmark case.

Results
In Fig. 18, the force-displacement curves obtained from the adaptive and non-adaptive simulations are shown. A small elastic regime is visible, where the layers are extended in thickness direction. Upon a significantly large out-of-plane normal traction, a delamination crack initiates at the interface and the curves adapt towards the LEFM solution. An explanation for the deviation in the first part of the damage regime is that the LEFM solution is based on clamped boundary conditions for both layers at the crack tip, while rotations are present in the simulations (material is deformable). These rotations are more significant when the crack length is small, hence the larger deviation for small displacements. The figure shows that both simulations predict the same response and this allows us to compare run times and simulation sizes later. In Fig. 19, the element states are shown for a number of time steps. In the first time step, all elements are in the lumped configuration, indicating perfectly attached layers. In subsequent steps, the two layers are opened, causing elements to upgrade to the discontinuous state via the intermediate layered configuration.

Efficiency gain
For this benchmark case, a total number of 80 incremental time steps has been used. The non-adaptive simulation contains a fixed amount of 10,908 degrees of freedom and runs for 2566 s. The simulation that utilizes adaptive elements, starts and ends with 5454 and 6471 degrees of freedom respectively and takes 403 s to complete. This yields a speed-up of 6.4 for this particular simulation when adap-

End notch flexure simulation
In the second benchmark case, an end notch flexure (ENF) specimen is analysed, see Fig. 20. The model has a length of L = 120 mm, a width of b = 20 mm and is t = 4 mm thick. The specimen contains an initial crack of a 0 = 46.9 mm (this value assures a stable crack growth) and is present along the full width at z = 0. The area moment of inertia is defined as I = 1 12 bh 3 . The ENF sample consists of two unidirectional composite plies of equal thickness, both having the fibers aligned with the global x-axis. The material properties are listed in Table 1.
The left (x = 0, z = −h) and right (x = L, z = −h) bottom edge boundaries are constrained in the y-and z-direction,  i.e. u y = u z = 0. These boundary conditions allow the ENF specimen to rotate at the supports and set it free to contract in the x-direction. At the top layer in the centre of the specimen (x = L 2 , z = +h), a step-wise prescribed displacement u z = λ(t) × u z,limit up to u z,limit = −4.0 mm is applied. In this case, the load factor is given by λ(t) = (t−1) 3 +1 2 , with λ(t) ∈ [0, 1] and t ∈ [0, 2]. This expression for the load factor results in reduced load increments when sudden crack growth is expected (to prevent the upgrading sequence from lagging). At the same time, the displacements in x-direction are constrained at this point to avoid rigid body movements.
The numerical results are compared to analytical solutions of Euler-Bernouilli beam theory in the elastic regime and LEFM in the damage regime: In order to get a physically correct behaviour, a cohesive penetration stiffness is used to prevent intrusion of the two fractured layers. The friction forces arising between the cracked layers are not considered in this benchmark case and no non-linear material behaviour is included in the simulations. The in-plane displacement field is approximated by second order bivariate B-splines and the displacements in thickness direction are also discretised by second order Bsplines. 175 × 1 (L × b) Bézier elements have been used in this example. The interlaminar shear stresses that are present drive the crack propagation. In order to upgrade the elements, the utilized stress criterion couples the shear stress at the interface to the ultimate cohesive shear traction. The transition from lumped to layered elements occurs at a stress level of 0.07 × t ult and the evolution from layered to discontinuous configurations takes place at a stress magnitude of 0.09×t ult . It is noted that the stress recovery model is not enabled in this simulation.
The cohesive zone model is provided with an ultimate traction of t ult = 50.0 N/mm 2 and the mode II fracture toughness has been set to G c = 0.50 N/mm.

Results
In Fig. 21, the force-displacement results are presented. The adaptive and non-adaptive outcomes are identical and are in excellent agreement with the analytical predictions, both in the linear elastic and damage regime.
In Fig. 22 the element configurations and interface damage parameter are shown for a number of time steps. The initial crack is represented by the discontinuous elements in the first time step. In subsequent time steps the elements are refined to allow the crack to propagate through the specimen. In this example, the stress thresholds that govern the upgrading sequence mostly affect the results at larger displacements (when the crack passed the center). Initially the response is linear, so the elements are not required to upgrade. During first crack growth (a ≤ L 2 ), the simulation is controlled by very small load increments, which assures a smooth and correct way of upgrading. However for larger displacements (a > L 2 ), the load increments increase and this makes the elements sensitive to lagging. Raising the adaptive settings (stress thresholds), leads to deviations that first appear in the final time steps.

Efficiency gain
During the ENF simulations, 100 equally sized time steps have been used. The non-adaptive simulation contains 9558 degrees of freedom and takes 971 s to complete. The adaptive one start with 6579 and ends with 7902 degrees of freedom and runs for 681 s before it successfully terminates. These results yield a speed-up of 1.4 when adaptive elements are used in the ENF simulation. The efficiency gain is less compared to the DCB case for a number of reasons. First, a significant amount of elements are initialized in the discontinuous configuration. These elements are already in the most expensive computational state and therefore this has a negative influence on the efficiency improvement. The second reason is very similar; a large part of the domain eventually gets delaminated, which means many elements must obtain the discontinuous configuration and this decreases the speedup.

Simply supported thick beam simulation
The final benchmark case is the simulation of a simply supported thick beam subjected to a bending load. The aim of this example is to first predict the correct location of interlaminar damage initiation by using the stress reconstruction model and subsequently investigate the damage propagation. The bi-layered beam model is shown in Fig. 23. It has a length of L = 9 mm, a width of b = 0.9 mm and is t = 1 mm thick. The laminate is built up from an elastic composite material, where the fibers are aligned with the global x-direction. The properties of the composite material are listed in Table 2.
The boundaries are simply supported at z = 0, which means that u x = u y = u z = 0 at the left boundary and u y = u z = 0 at the right boundary. At the top layer, the  beam is loaded by a triangular distributed load given by: with p 0 = 15 N/mm 2 . The maximum force exerted on the top surface equals F z = 1 2 p 0 b L. During the non-linear simulation, the Neumann boundary condition is applied according to p z = λ(t) × p z,limit , where the load factor λ(t) ∈ [0, 1] is controlled by a Riks arc-length solution procedure. This method is capable of accurately capturing instabilities that would lead to snap-back (displacement-controlled) or snapthrough (force-controlled) behavior in a total Lagrangian calculation. With the current choice for the thickness-tolength ratio of the beam, a mode II delamination crack is expected to initiate simultaneously at both supports. For more details, the reader is referred to [6].
The total applied force by the distributed load, F z , is related to the deflection of the centre of the beam u z . The numerical results are compared to Timoshenko's beam theory: where A = bt is the cross-sectional area and κ = 10(1+ν 12 ) 12+11ν 12 is the shear coefficient according to Cowper [25]. The second moment of inertia is given by: The in-plane displacements are described by third order bivariate B-splines, whereas the out-of-plane displacements are discretised by univariate second order B-splines. Due to the relatively small geometry of the beam, a total number of 30 × 1 (L × b) uniform Bézier elements is sufficient.
In order to assess whether the stress reconstruction model captures both the correct location of damage initiation and corresponding delamination failure mode, the failure criterion proposed by Ochoa et al. [26] has been used to evaluate whether elements should change in configuration: where · + represents the Macaulay bracket. The upgrading thresholds f I have been chosen, such that the transition from lumped to layered occurs at f I = 0.45 and that the next change, from layered to discontinuous, occurs at f I = 0.55. The fracture toughness for this case has been set to G c = 0.20 N/mm and the ultimate traction t ult = 25 N/mm 2 .

Results
The force-displacement plot is shown in Fig. 24. The adaptive and non-adaptive outcomes are not identical, yet they are in agreement. It is visible that the initial response of both the adaptive and non-adaptive simulations is not entirely elastic. The deviation between these results and the analytical solution of Timoshenko is due to accumulation of some interlaminar damage already here. For comparison, a linear elastic simulation with the same beam geometry but without any interface was conducted. As can be seen, the force-displacement results from this simulation are in good accordance with the Timoshenko elastic stiffness (blue dashed line in Fig. 24). Moreover, the adaptive simulation is stiffer at the beginning where lumped elements are present for the greater part. The non-adaptive simulation adds more compliance, due to the predefined cohesive interface. The sudden load drop, which is conceivable due to the arc-length method (negative incremental load factor λ(t)), indicates a sudden crack growth. At the point where the beam exhibits a significant amount of damage, the top and bottom layers are separated and start to slide over each other. The numerical results then  converge towards the analytical Timoshenko solution that holds for a beam with a fully cracked interface.
In Fig. 25, the adaptive beam is shown in the deformed configuration. Please note that the load factor λ(t) is now also able to decrease during the calculations. The simulation starts with lumped elements, where the stress reconstruction model is active with the intention of predicting the out-ofplane stresses accurately. The first elements upgrade exactly at the left and right boundaries as this is where the transverse shear stresses are dominant. This leads to the initiation of a mode II fracture at the supports. Further on, the elements refine towards the center of the beam, matching nicely with the crack propagation. In the final load steps, all elements are upgraded to the discontinuous configuration and the crack extends almost throughout the complete patch.

Efficiency gain
Although the number of time steps for both simulations are not equally chosen by the arc-length method, there is no significant difference preventing us from comparing the simulations in terms of efficiency. The non-adaptive simulation, counting 2376 degrees of freedom, finishes successfully after 1266 s. The adaptive one starts with 1188 and ends up with 2376 degrees of freedom and takes only 634 s to complete. These outcomes result in a run time speed-up of 2.0 when adaptively refined elements are considered.

Conclusions
In this paper we presented an isogeometric continuum shell element in which the displacement interpolation through-thethickness of the shell can be modified automatically. Based on an interface stress criterion, the considered element can be upgraded while the simulation is ongoing in order to improve the accuracy of the element under high stress states or to model delamination cracks.
The use of isogeometric shape functions is essential here. It enables to introduce weak and strong discontinuities at the layer interfaces by means of knot insertion. In addition, the higher order continuity of the B-spline shape functions allows for an element-wise reconstruction of the rather poor out-of-plane stresses of lumped elements. These enhanced stress results can then be utilized in the stress criterion.
Simple numerical benchmark problems containing a single interface have already shown significant reductions in the simulation time without compromising in precision. The element currently accounts for single interface models and still needs to be extended to account for the initiation and progression of multiple delamination cracks. This approach remains to be demonstrated in large scale analyses of layered composite structures.