An adaptive shell element for explicit dynamic analysis of failure in laminated composites Part 2: Progressive failure and model validation

To enable modelling of the progressive failure of large, laminated composite components under crash or impact loading, it is key to have a numerical methodology that is both efficient and numerically robust. A potential way is to adopt an adaptive method where the structure is initially represented by an equivalent single-layer shell model, which during the analysis is adaptively transformed to a high-resolution layer-wise model in areas where higher accuracy is required. Such a method was recently developed and implemented in the commercial finite element solver LS-DYNA, aiming at explicit crash analysis (Främby, Fagerström and Karlsson: An adaptive shell element for explicit dynamic analysis of failure in laminated composites Part 1: Adaptive kinematics and numerical implementation, 2020). In the current work, the method is extended to the case of interacting interand intralaminar damage evolution. As a key part, we demonstrate the importance of properly regularising the intralaminar failure described by a smeared-crack model, and show that neglecting to account for the crack-versus-mesh orientation may lead to significant errors in the predicted energy dissipation. We also validate the adaptive approach against a four-point beam bending test with matrix-induced delamination growth, and simultaneously show the capability of the proposed method to – at lower computational expense – replicate the results from a refined, non-adaptive model.


Introduction
In their ambition to reduce vehicle weight, the automotive industry is aiming to increase the share of high-performance carbon Fibre-Reinforced Polymers (FRP) in structural parts. Carbon FRP are interesting because of their superior specific properties in terms of stiffness, strength and energy absorption (in crushing), cf e.g. Carruthers et al. [1], Jacob et al. [2] or Park et al. [3]. However, the latter property is strongly dependent on that the proper failure modes are triggered, and that severe delamination is mitigated (cf e.g. Hull [4] and Grauers et al. [5]). As the development of new cars is driven by simulations, and is strongly constrained by safety regulations, increasing the share of carbon FRP needs to be assisted by accurate finite element (FE) modelling tools that are efficient enough to allow full car crash analyses. To enable such FE simulations at larger scale, most often with an explicit time stepping algorithm, we have in an accompanying paper (Part I) [6] proposed an adaptive refinement method for modelling the failure kinematics of laminated structures. We acknowledge that similar methods have been proposed for progressive laminate failure e.g. Shor and Vaziri [7], McElroy [8] and Adams et al. [9], although we believe that this is one of the first method that is adapted to a dynamic setting and which deals with the associated numerical challenges. In our previous paper, we also detailed the implementation of the proposed method in the commercially available FE software LS-DYNA, commonly used in crash analysis. The focus of the Part I was on the interlaminar fracture. In this second part of the work, we expand the adaptive refinement method with models to describe the intralaminar behaviour and handle the interaction between delamination and intralaminar failure. Such interaction mainly involves delamination initiation caused by interface stress concentration at the tips of matrix cracks. Ideally, this requires modelling of both the intralaminar and the interlaminar cracks as strong discontinuities, cf e.g. Fang et al. [10], Reiner et al. [11], Lu et al. [12] and Zhi and Tay [13]. However, for simplicity we propose to model intralaminar fracture using a smeared-crack material model, where no explicit intralaminar crack is present. Therefore, to model the intra-and interlaminar interaction we adopt an approach suggested by Yun et al. [14], although herein adapted for increased numerical robustness. This interaction works by degrading the critical fracture energy of a cohesive zone (CZ) element by the damage state of the adjacent intralaminar elements.
To properly describe the energy dissipation due to intralaminar failure with a smeared-crack model, the fracture energy must be regularised with respect to the element size and crack orientation. Already 1985, Bažant [15] showed that for an isotropic 2D case and a finite element mesh of right angled rectangular elements, not only the element size, but also the orientation of the localising strain band (representing the crack), needs to be properly accounted for. Despite this, the crack orientation is often neglected, oras we argue in this paperis inadequately accounted for. Therefore, in this paper we extend Bažant's idea to the case of generally shaped continuum shell elements describing fibre-reinforced laminates with different ply orientations. We show the importance of considering both element shape and crack orientation for various cases and element shapes.
Finally, we validate the adaptive concept including interacting interlaminar and intralaminar failure against a four-point bending test reported in the literature. The validation shows that our adaptive refinement method can model the experiment equally well as a reference (non-adaptive) simulation.
By expanding our work from [6] we believe that we take a large step towards realising accurate and computationally efficient industrial crash simulations. This is, as mentioned above, crucial for structural composites to have a widespread use in future cars [16].

Outline of paper
In a first section we make a short summary of the adaptive refinement method followed by presenting intralaminar-based refinement indicators. Next, we will focus on the intralaminar material model, the proper regularisation of this and the interaction of intra-and interlaminar fracture. We then present some numerical examples to show the performance of our proposed method. We end the paper by summarising our work and try to set it into context of other contributions in the field.

Method for adaptive fracture modelling of fibre-reinforced polymers
As mentioned in the introduction, in [6] we presented the implementation of an adaptive refinement method into the software LS-DYNA. The method, which is based on our work in [17], consists of the following steps: i The laminated structure is initially represented by a equivalent single-layer (ESL) of solid shell elements through the thickness. ii Refinement indicators are used to locate elements where the fracture process cannot be resolved. iii These elements are then locally refined through the thickness by enriching the element kinematics to account for material interfaces (weak discontinuities). iv At interfaces prone to delaminate, CZ elements are inserted such that delaminations (strong discontinuities) can initiate and propagate.
This means that every adaptive element in the model can be in one of three stages: a) Unrefined; b) Refined with one or several discrete material interfaces; c) Discontinuous, having one or several refined interfaces separated and CZ elements inserted. This is illustrated in Fig. 1. The first refinement step, where material interfaces are added through the thickness, we refer to as a weak refinements and the secondary refinement step, where CZ elements are inserted, we refer to as strong refinements.
In order to limit the computational costs, the method does not model intralaminar failure with separate DOF. That is, only throughthe-thickness (and not in-plane) refinements are made. Instead, intralaminar fracture is modelled using a smeared-crack material model, which will be the focus of Section 3. Nevertheless, the through-the-thickness refinements will improve the prediction of the stress state in the element and thus also improve the prediction of intralaminar fracture.
The idea of the through-the-thickness refinement is that we can make a two-stage relaxation of the kinematic constraint of the ESL shell. In the first weak refinement stage, when we introduce weak discontinuities, we transform the ESL shell to a layerwise (LW) element whereby the stress prediction can be improved. In the second strong refinement stage, the weak discontinuities can be separated and CZ elements are inserted. This way, delamination initiation and propagation can be modelled.

Refinement indicators
In [6] we presented an interlaminar criterion f CZ , which indicates when strong refinements (with CZ elements) should be introduced to model potential delamination growth. Here, we will complement this with intralaminar criteria, which indicates when the intralaminar state requires element refinement in order to be accurately resolved.

Stress-based indicator
The first indicator is based on the intralaminar failure initiation criteria f I , where I represents the different intralaminar failure modes, i.e. matrix tension and compression (MT and MC) and fibre tension and compression (FT and FC). That is, an element is weakly refined whenever a criterion reaches a certain limit To avoid sudden jumps in the stress state upon refinement, the refinements should be introduced before any failure is initiated both in the unrefined and the refined state of an element. The issue with this is that the full stress state is not well captured by the ESL shell. In general, only the in-plane stress components are predicted with sufficient accuracy. Thus, the stress state of the refined element is not available. However, as we showed in [17], a stress recovery technique can be used to estimate the full stress state in the refined element. Therefore, the refinement criterion in Eq. (1) can be evaluated using both the stress state of the unrefined element as well as a recovered stress state in order to avoid sudden failure of the refined element. As mentioned in [6], the recovery technique has not yet been added in the current LS-DYNA implementation.  1. The different stages of the adaptive elements: In the unrefined stage the entire laminate is represented by one element through the thickness (a); The element can be through-the-thickness refined in order to model weak discontinuities, here exemplified with three refinements (b); The refined interfaces from the previous stage can be separated (and have CZ elements added) in order to model strong discontinuities (c).

Damage-based indicator
Delaminations are initiated and propagated by high out-of-plane interlaminar stresses. These stresses are mainly the results of either high external loading conditions or sharp discontinuities like edges or ply cracks. In the latter case of ply-crack-induced delaminations, intralaminar cracks will create a stress concentration in the adjacent interfaces which can cause delaminations to initiate and propagate. This is illustrated for matrix cracks in Fig. 7. Thus, even though the out-of-plane stresses are predicted to be low (meaning that the interlaminar criterion f CZ from [6] is not met), CZ elements must be introduced. To achieve this also the damage levels d I of all plies are monitored. If the damage in a ply reach a certain limit where once again I represent the different failure modes of the plies, strong refinements with CZ elements are added in the adjacent interfaces. Remark: The CZ elements introduced at an intralaminar crack will likely directly go to a damaged state. Thus, we will violate the requirement to refine before failure is initiated, however, this error will be limited to only a few elements.

Limiting the amount of refinements
In general, weak and strong refinements are possible at any interlaminar interface of an element. This may lead to many refinements through the thickness and thus also high computational cost. As a mean to limit the number of refinements we have added an optional feature, which suppresses refinements in interfaces where the pitch angle 1 Δθ, is low. That is, if a refinement is indicated at any point in the laminate, a search is made down-and upwards for interfaces with pitch angle which are then triggered for refinement. Even though refinement suppression is mainly motivated from a computational efficiency perspective, there is some evidence that the most critical delaminations will occur in interfaces with large pitch angles. For example, in the work on Bouligand CFRP laminates by Mencattelli and Pinho [18] they concluded that for small pitch angle (|Δθ i |⩽10 ∘ ) there is an increasing amount of diffused subcritical damage which do not lead to large delaminations.

Fig. 2.
A smeared-crack model uses the characteristic length L e to translate an intrinsic surface traction-separation law (defined by the maximum traction S 0 , the artificial penalty stiffness k and the critical fracture energy c ) into a volume stress-strain relationship (defined by the maximum traction S 0 , the elastic stiffness E and the scaled critical fracture energy c /L e ). In Table 1 we have listed all parameters which control the behaviour of the adaptive refinement, together with suggested default values. Most of these values have been found during verification and validation in [6], however, we will discuss the appropriate values of f ref,lim and d lim in Section 4 below.

Models for progressive failure
In Part I [6] we presented how the interlaminar fracture is modelled with CZ elements. Here we expand the work by focusing on the intralaminar fracture, which we choose to model using a smeared-crack material model [19]. Below we will briefly present the particular model used in this paper, including detailing the proper regularisation needed to make the model mesh objective. Furthermore, we address the challenge of modelling interaction of intra-and interlaminar fracture and how this is treated in our adaptive refinement method.

Intralaminar material model
The idea of a smeared-crack model is that a crack surface, described by an intrinsic traction-separation law (TSL), can be smeared over a representative finite element volume to define a stress-strain relation. This smearing, which is achieved using a length parameter L e that relates the crack volume to the crack area, is illustrated for a bi-linear TSL in Fig. 2. Since the maximum traction S 0 and the elastic stiffness E are material properties, the critical fracture energy c of the TSL must be scaled by L e , such that the correct energy is dissipated by the stress-strain relation. In practice, this means that the failure strain ε f is dependent on L e as The proper calculations of the so-called characteristic length L e will be the focus of Section 3.1.1 below.
In this paper we have used a simplified intralaminar material model which only consider tensile matrix cracks using a single damage variable d. In a large deformation setting, the material model relates the second Piola-Kirchhoff stress S to the Green-Lagrange strain E. To calculate the Green-Lagrange strain, represented as the matrix E in the reference lamina frame ℒ 0 , the model follows the same procedure as the LS-DYNA option IHYPER = − 1 [20]. The base of the deformation gradient matrix F, given in global coordinates, is rotated to the lamina reference frame as where Z 0 is the rotation from reference global 0 to reference lamina frame ℒ 0 defined in Appendix A. The deformation gradient F, which thereby maps from reference lamina to current global, is used to calculate the lamina Green-Lagrangian strain as The elastic lamina second Piola-Kirchhoff stress, on Voigt form represented by the vector S V is then determined as where E V is the Voigt vector representation of E and C V is the Voigt matrix form of the material stiffness with the inverse The fibre direction normal stress is not included in damage evolution model and has therefore been omitted here.
In Eq. (8), E i , G ij and ν ij are the Young's and shear moduli and Poisson's ratios with respect to the lamina longitudinal (1), transverse (2) and out-of-plane (3) directions.
In order to evaluate failure initiation and damage evolution the model follows the same procedure as in [21], where the elastic (trial) transverse and out-of-plane lamina stresses and strains are rotated to the potential matrix crack plane, inclined an angle ϕ to the out-of-plane direction, as illustrated in Fig. 3. The failure initiation criterion is then evaluated using the stresses on the crack plane as where 〈 • 〉 + are the positive Macauley brackets and Y + , S TT and S LT are the maximum matrix tensile and shear stresses. That is, in each time step the fracture plane angle ϕ which maximises the criterion in Eq. (9) is sought. If the criterion in Eq. (9) reaches unity the crack plane is considered fixed and the crack angle is stored as ϕ 0 . Furthermore, the failure onset resultant stress Ŝ 0 and strain Ê 0 are saved along with the angle between the shear stresses λ 0 and the angle between the resultant shear stress and the normal stress ω 0 . The effective failure strain is then computed as where ̂c ,MT is the effective fracture toughness for matrix tension and L em is the matrix crack characteristic length defined in Eq. (22) in the Section 3.1.1 below.
In each time step following the failure onset the elastic trial stresses (evaluated using Eq. (7)) are rotated to the crack plane ϕ 0 and the effective damage-driving stress and strain are computed using the stress angles λ 0 and ω 0 as where | • | represent the absolute value, and Using these, the instantaneous damage is calculated as and to account for irreversibility the damage variable is defined as The damage d is used to degrade the trial stresses on the crack plane together with the Poisson's ratios. The degraded crack stresses are then rotated back to the lamina frame. The procedure above will make the effective stress Ŝ on the crack plane follow a bi-linear law similar to that in the right side of Fig. 2, i.e. the damage variable has a value of 0 at onset of failure (Ê =Ê 0 ) and 1 at final failure (Ê =Ê f ).
Since LS-DYNA expects the Cauchy stress to be returned we use the deformation gradient F to push forward the updated lamina second Piola-Kirchhoff stress where

Proper regularisation of fracture energy
According to Eq. (4), a smeared-crack material model needs to be regularised with a characteristic length L e such that the correct amount of energy is dissipated. To calculate this length there are primarily two strategies. L e can be considered as a mesh property, where you calculate the volume of the failed elements as the crack passes through the mesh. Bažant [15] derived this type of meshbased characteristic length for a crack passing at an angle θ through a 2D mesh of rectangles with side lengths L ξ and L η as as illustrated in Fig. 4a. The expression in Eq. (19) assumes that the crack is perpendicular to the mesh surface. For cracks inclined an angle ϕ to the mesh surface (cf Fig. 5), Eq. (19) can be expanded to where h is the thickness of the element. Alternatively, L e can be considered as an elemental property, where the imagined crack surface area A c in each element is smeared over the element volume V e as which is illustrated in Fig. 4b. An example of this can be found in the work by Chiu et al. [22]. The potential drawback of the mesh-based expression in Eqs. (20) and (19) is that in its original form it assumes a mesh of right angled (rectangular) elements. However, a clear benefit is that it considers a problem which arises when the crack growth is not aligned to the mesh lines. In this case, overlapping elements are activated in order to model a crack through the mesh, as can be observed in Fig. 5. That is, the fracture energy is regularised both with respect to the element size and the orientation of the crack  (20)). In this particular case with uniform and regular mesh, the characteristic lengths will be the same in all elements. NB: The dashed line is only added for illustration purposes and does not describe the actual position of the crack. (b): Alternatively, the characteristic length can be treated as an elemental property where the overlap of elements is not taken into account. Instead the characteristic length is simply the ratio of the element volume with respect to the crack area in each individual element (Eq. (21)).
growth in relation to the mesh lines. This last part is ignored when treating L e as a pure elemental property in Eq. (21), even if the elemental crack area A c is calculated with the correct crack orientation in mind. Notably, if using the element-based expression in Eq. (21) to estimate the characteristic length, the energy dissipation per crack area is over predicted when the crack growth path does not align with the mesh lines (up to doubled 2 ). Therefore, we will below expand the mesh-based characteristic length expression in Eq. (20) with additional modification for non-rectangular elements.
In an FRP, matrix cracks can form parallel to the fibre direction while fibre fracture will primarily occur perpendicular to the fibres. Thus, two corresponding characteristic lengths L em and L ef need to be calculated for each unique ply orientation. Furthermore, for solid shell elements the bottom and top surfaces are not necessarily identical and the side lengths can vary through the thickness. Thus, each element ply requires a separate set of characteristic lengths L em,i and L ef,i . To estimate these, the following assumptions are made: • In each ply of an element there can only be one matrix crack and one fibre crack.
• The mesh is fairly regular, which includes that the in-plane shape of each element layer is fairly rectangular and the thickness of the element is fairly constant. • The geometry of the element changes little prior to fracture.
This way the characteristic lengths in each element ply can be calculated 3 following Eq. (20) based on the undeformed element geometry. The characteristic lengths in ply i then become   6. Illustration of how a general quadrilateral surface (left) is simplified to a rectangle with side lengths L ξ and L η (right). 2 When the crack is perpendicular to the mesh surface and crosses the diagonal of the element. 3 Even though we consider the characteristic length as a mesh property.
where L ξ,i , L η,i and h i are the side lengths and the thickness of the ply and θ i is the ply fibre angle with respect to the ξ-axis of the element (cf Appendix A for definition). In Eq. (22) we have assumed that fibre cracks always are perpendicular to the midsurface, i.e. ϕ = 0 4 . The matrix crack angles ϕ i are the results of the stress state of the element as described in Section 3.1. Thus, they are determined during the simulation. However, L ξ,i and L η,i can be calculated at the start of the simulation.
Since the in-plane geometry of the elements is not necessarily rectangular, the side lengths need to estimated, meaning that there will be an error for elements which deviate from this shape. For a general quadrilateral element surface the side lengths are estimated by projecting the surface bimedians a 1 and a 2 on the element coordinate axes ξ and η, respectively, as shown in Fig. 6. Since each element layer is an interpolation of the bottom and top surfaces, only the side lengths at the bottom and the top of the element are estimated, i.e.
where ‖ • ‖ represent the vector norm and ψ is the angle between the bimedians. Then the side lengths of each ply in the element can be interpolated as where ζ i is the position of the ply (mid) surface and ζ bot and ζ top are the bottom and top coordinates along the element normal ζ. Next, we need to calculate the angle between the fibre direction and the mesh lines. Following the definitions of the element and the lamina coordinate systems in Appendix A this is simply a rotation around the midsurface normal where β i is the ply fibre direction with respect to the material a-axis (defined from analysis input) and ω is the angle between the material and element frames, calculated from the element-to-material rotation matrix Q as To end this subsection, we want to emphasise that even with proper regularisation of the fracture energy, smeared-crack material models will show a spurious sensitivity to develop damage growth along the mesh lines of the FE-model. That is, this problem is not avoided by using a mesh-based characteristic length calculation. Instead, to overcome this sensitivity, the material model needs to be extended with some type of crack-tracking functionality which controls the growth direction, e.g. that by Mukhopadhyay and Hallett [23]. While crack-tracking does not solve the fundamental ill-posedness of the smeared-crack approach, it forces the crack to follow certain directions, e.g. that a matrix crack must be parallel to the fibre direction. In the results section we will elaborate on the consequence of using crack tracking or not together with the different options for calculating the characteristic length. For completeness we have added the calculation of the element-based characteristic length in Appendix B.

Interaction between smeared ply crack and delamination
Matrix cracks, which can form at comparatively low loads, are common initiation points for delamination. In Fig. 7 the process of Fig. 7. The process of matrix-crack induced delamination: First, a transverse matrix crack starts to form (a); When the matrix crack is fully evolved it will lead to a shear stress concentration at the interface (b); If the shear stresses are high enough delaminations can propagate (c). Adapted from [24]. so-called Matrix-Crack Induced Delamination (MCID) is illustrated, where high interlaminar shear stresses at the vicinity of a matrix crack tip lead to the initiation of delaminations.
Capturing the interaction of intra-and interlaminar failure is challenging when using a smeared-crack approach to model the intralaminar fracture. This since there will be no explicit intralaminar crack present to generate a stress concentration. The problem is illustrated in Fig. 8 where we also compare two methods for modelling intralaminar cracks and how they interact with adjacent interlaminar CZ elements. If the intralaminar crack is smeared, the neighbouring CZ element will not exhibit any significant shear stresses.
To remedy this problem we have followed the approach proposed by Yun et al. [14], where they argue that a significant interaction exists between delamination and intralaminar failure. To capture this in an FE model, where a smeared-crack intralaminar damage model is used, they suggest to reduce the critical fracture energy of a CZ element by the damage state of the adjacent intralaminar elements, i.e.
In Eq. (27) d intra is the maximum intralaminar damage of the adjacent ply elements. As the intralaminar damage evolves the fracture energy of an adjacent CZ element will decrease (potentially down to zero). To achieve a reduction in cohesive fracture energy the penalty stiffness is held constant, which will result in a decreased failure separation. The consequence of this reduction is illustrated in the left side of Fig. 9. To avoid snap-back in the response we have, in contrast to [14], added a reduction of the cohesive strengths if the reduced fracture energies in Eq. (27) would result in a failure separation lower than the initiation separation: This is illustrated in the right side of Fig. 9. An alternative to reducing the failure separation, would be to decrease the penalty stiffness (while keeping the onset and failure separation constant) as in Quintanas-Corominas et al. [26]. However, due to implementation reasons this approach was not used. In practice, the reduction in Eq. (27) is a smooth removal of the CZ element adjacent to an intralaminar crack. While there might exist more advanced methods yielding better predictions of the intra-and interlaminar interaction, we have adopted the approach of reducing the fracture energy since it is very easy to implement.
This approach will to a large extent neglect the fracture energy dissipated by the delamination in the element at the intralaminar crack. However, as we pointed out in the introduction, even though delamination is one of the governing failure mechanisms in the crushing of laminated FRP it is not a primary energy absorbing mechanism. Instead, delaminations influence the overall deformation pattern and thereby the occurrence of other failure mechanisms. We therefore conclude that for general engineering applications it is more important to capture the interaction between intra-and interlaminar failure than the correct energy dissipation in an individual CZ element.

Results
In this section we will present some numerical examples to show the performance of our proposed method. The method is still a proof of concept and some key aspects necessary for performing simulations of crash problems, e.g. including the refined surfaces in the global contacts, are out of scope for the present work. Furthermore, crash problems are naturally chaotic and we therefore instead show examples where different aspects of fracture modelling are isolated. The intent of the numerical examples is to show that the concepts presented in the present article can be used in conjunction with our adaptive modelling method, presented in Part I [6]. We will first focus on how the choice of characteristic length calculation can affect the results and then we will validate the adaptive approach against experimental test.

Regularisation of fracture energy
In Section 3.1.1 we covered how a smeared-crack material model needs to be regularised with a characteristic length L e such that the correct amount of energy is dissipated. Furthermore, we explained that we treat the matrix and fibre characteristic lengths as mesh properties (as opposed to element properties) and showed how we can calculate the lengths in each element.
In this section we will show how the numerical results are affected by the characteristic lengths. The calculation of the matrix characteristic length in Eq. (22) is dependent on the crack inclination ϕ, which is a result of the loading conditions and thus not available in the beginning of the simulation. Therefore we will in this part only consider (matrix and fibre) cracks perpendicular to the midsurface (ϕ = 0).

Comparison to numerically calculated characteristic length
In the first example we want to show the ability of the mesh-based characteristic lengths calculation in Eq. (22) to obtain a correct estimation. To achieve this we run a series of analyses on 2D meshes of 101 by 101 rectangular elements with different side length ratios L ξ /L η (L η is held constant), where the angle θ between the mesh and cracks is varied. The correct (numerical) characteristic lengths are computed by dividing the total volume of all elements needed to represent an inclined crack path by the total area of a corresponding crack (length × thickness). At the same time, in each individual element, we compute the mesh-based and the elementbased characteristic lengths using Eq. (22)    factor of two. For example for the square element mesh in Fig. 10, a crack passing at θ = 45 ∘ will yield a numerically calculated characteristic length of ̅̅̅ 2 √ L ξ . The mesh-based characteristic length (Eq. (22)) is the also ̅̅̅ 2 √ L ξ , while the element-based characteristic To also show the performance on a mesh of non-rectangular elements, we vary the shape of the elements by distorting the element nodes in opposite diagonal directions, which results in a mesh of identical non-rectangular elements. As a consequence of nonrectangular elements, the side element lengths L ξ and L η need to be estimated according to Eq. (23).
To quantify the element shape distortion we compute the determinant of the Jacobian matrix J of the element isoparametric mapping using 2 x 2 Gauss quadrature points (see e.g. Ottosen and Petersson [27]) and take the minimum over the maximum value of all points: The distortion value in Eq. (29) is 1.0 for a rectangular element and the element will become concave when r det(J) < 0.27 (due to the choice of the Gaussian quadrature scheme). In Figs Finally, to show the performance on a non-repeating mesh, we randomly distort the nodes up to half an element side length in any in-plane direction and compute the characteristic lengths. In Fig. 14 we have compared the average mesh-based 5 and element-based characteristic lengths for the cracked elements to the correct characteristic lengths. Even though we are comparing average values it is clear that the mesh-based characteristic lengths in Eq. (22) matches the correct lengths also in this case.
In summary, if the characteristic lengths are treated purely as element properties the energy dissipation would be over predicted by up to a factor of two. However, by treating them as mesh properties and estimating the side element lengths for non-rectangular elements, very accurate predictions can be obtained.

Mesh direction sensitivity
The calculation of the characteristic lengths calculation in Eq. (22), assumes that the crack growth will depend only on the fibre direction. However, as we have already mentioned the smeared-crack material models can suffer spurious sensitivity to develop damage growth along the mesh lines. Thus, if this is not handled, e.g. by a crack-tracking technique, the results might be inaccurate.
To evaluate the mesh direction sensitivity we consider an off-axis tensile test where a plate, with length L = 20 mm, width W = 10 mm and thickness h = 1 mm, is subjected to uniaxial displacement Δ x along the length direction. The fibres are oriented at an angle of θ = 45 • to the loading direction and, as illustrated in Fig. 15, oblique ends of 63 • are used to enable a uniform stress field, cf Sun and Fig. 15. Geometry for off-axis tensile test used to evaluate the calculation of characteristic length and mesh direction sensitivity.

Table 2
Ply material properties used in off-axis tensile test example. Please note that isotropic fracture toughness (and strength) is assumed. Chung [28]. The material properties used for this example are given in Table 2.
The plate is meshed with quadrilateral shell elements (solid shell element formulation 3 in LS-DYNA) and the mesh lines are progressively skewed towards the ends to enable a regular mesh in the central 10 × 10 mm 2 part of the plate. Three different models are used: i) With a fixed crack path, which corresponds to using a crack-tracking technique 6 ; ii) With a free crack path, where the central 2 × 2 elements have reduced strength (30 MPa) to trigger fracture; and iii) A reference model where a fixed crack is modelled using cohesive elements (solid elements formulation 19 with material model MAT138). The FE-models are illustrated in Fig. 16.
The first two models are run twice: first using an element-based characteristic length, corresponding to Eq. (B.3), and secondly using the mesh-based characteristic length in Eq. (22). Since the correct crack growth path is 45 • to the mesh lines, the element-based characteristic length will be half of the mesh-based (cf Fig. 11), which should lead to twice as much energy dissipation. To evaluate the results we will compare load versus displacement of the five different cases and the level of the dissipated energy in the crack.
In Fig. 17, the load-displacement response is plotted for all five cases. We can notice that the cases with an element-based characteristic length will result in higher maximum load and displacement, compared to the corresponding mesh-based characteristic length cases. In addition, the results in Fig. 18 show that the element-based cases dissipate twice the amount of energy in the cracks compared to the mesh-based cases (as expected). This figure also confirms that the reference model dissipate the same amount of energy as the analytical value, As shown in Fig. 19, when no crack tracking is employed the crack grows perpendicular to the loading direction along the mesh    lines, regardless of the choice of characteristic length type. This leads to an underprediction of the maximum force and displacement levels. Furthermore, the amount of dissipated energy is approximately 40% too low for the mesh-based characteristic length while the element-based matches the correct dissipation level better (11% to high). If crack tracking is used, it is evident that the model with mesh-based characteristic length match the reference results well (within 4%) while the element-based model overpredicts the force response and especially the energy dissipation (112% to high). The results above are in line with the reasoning in Section 3.1.1, i.e. when the crack growth does not follow the mesh lines a meshbased characteristic length calculation should be used. However, if the correct crack growth cannot be ensured and growth is expected to always be along the mesh lines the mesh-based characteristic length can underpredict the force and energy dissipation. Then, it might be better to use an element-based characteristic length instead, since this will predict the correct response better. The choice of characteristic length is therefore largely dependent on which phenomena the analysis is to be able to capture correctly. If the correct crack path is important to capture, a crack-tracking algorithm will likely be necessary, which must then be combined with a meshbased characteristic length.

Matrix-crack induced delamination
In this example we demonstrate the ability of our proposed adaptive refinement method to simulate fracture of laminated FRP by simulating matrix-crack induced delamination in a four-point beam bending. The model set-up is the same as in the experimental study performed by Mortell et al. [24].
The beam, which has a layup of [90 4 /0 7 /90 4 ], is illustrated in Fig. 20. The bottom 90-layers of the beam are modelled with the smeared-crack material model given in Section 3.1, with the limitation of only allowing cracks perpendicular to the beam midsurface (ϕ = 0). Since no inclined cracks were identified in the experiments this is considered a reasonable simplification. The rest of the beam is modelled as an orthotropic linear elastic material, i.e. final failure of the beam when the 0-layers fail is not included. Furthermore, only delamination of the bottom 90/0 interface is considered and cohesive delamination failure is prevented at the ends of the beam (up to 1 mm from the edge).
At the moment the implementation of the adaptive method does not include any sub-element deletion technique. This has the implication that if cracks form in two adjacent elements their bottom shared nodes will be unconstrained, which can lead to negative element volumes in the analysis. If this happens, LS-DYNA will terminate the simulation -a feature that cannot be deactivated. Therefore we have enforced a minimum crack-spacing distance of two times the element length, which will prevent negative volumes. In [24], Mortell et al. observed crack spacing of approximately 2 mm, thus we are (for this reason) limited to use elements with in-plane side lengths below 1 mm to minimise the effect of the crack-spacing enforcement.
The ply material is HTA/6376 and the material properties, given in Tables 3 and 4, are taken from Reiner et al. [29]. However, please note that the interlaminar strengths have been reduced by 50% in order to allow the use of elements with an in-plane length of 0.25 mm and still have a good representation of the cohesive delamination failure (the suitable strength reduction is found through simulations and strategies from [30]). Furthermore, the beam density is assumed to be ρ = 1, 600 kg/m 3 and a cohesive density of ρ CZ = 0.8 kg/m 2 has been added to increase the critical time step.
The beam is modelled with one user element through the thickness. However, the adopted solid shell kinematics is not suitable to model bending with one element through the thickness, cf details in [6]. Therefore, an initial weak refinement has been activated at the eighth interface (from the bottom). In practice, this means that the beam actually consists of two solid shell elements through the thickness, although described internally in the user element. Besides the initial weak refinement, three different cases are examined: 1) an initial strong refinement of the bottom 90/0-interface; 2) an adaptive refinement of the bottom 90/0-interface without stabilisation and 3) an adaptive refinement of the bottom 90/0-interface with stabilisation. Here stabilisation refers to the procedure presented in Part I [6], where non-physical oscillations related to the introduction of refinements are removed. For the weak refinements, this is achieved by applying a mass-weighted nodal damping (DAMPING_GLOBAL) for a short period while for the strong refinements, an approach by Menouillard and Belytschko [31] is used. For the adaptive cases the 90/0 weak refinement is activated when the matrix tension failure initiation criterion f MT reaches 0.3 and the strong refinement when it reaches 1.0 (i.e. when intralaminar damage initiates). By setting Δθ del > 0, refinement activation in any of the bottom 90-layers will result in a refinement of the bottom 90/0-  Table 4 Interface material properties used in the MCID example. interface. The load is applied by displacing the two top load pins downwards at a constant speed of 0.4 mm/ms (with a linear ramp over the first 1 ms) and the total simulation time is 21 ms.
To stabilise the weak refinement a small amount of damping is applied over a short time period [6]. By first analysing the oscillations in the non-stabilised adaptive simulation, it can be found that the smallest (dominating) angular frequency is ω n = 556 rad/ms.
Choosing a damping ratio of 5% thereby means adding damping with a damping constant of C d,ref ≈ 56 rad/ms, according to the procedure proposed in [6]. To have as little effect of the damping as possible on the global beam response, we let the damping degrade linearly 7 over the damping duration Δt ref = 0.2 ms. The strong refinement is stabilised by gradually releasing the nodes of the CZ elements [6]. The separation is done over 100 time steps, which in this case yields a duration of Δt CZ = 0.0024 ms.
To validate the simulation results we compare the load vs. displacement curves from the experiment and the so-called Local Delamination Ratio (LDR) [24]. LDR compares the length of an individual delamination e to the through-thickness height h crack of the transverse intralaminar crack which initiated the delamination (cf Fig. 21): Fig. 22. Load-displacement curve for the non-adaptive reference simulation, where a strong refinement is activated in the bottom 90/0-interface from the start of the simulation, compared to the experiment [24]. At the experimental failure displacement (Δ z = 6.6 mm), the load is approximately 11% lower and the total work is 7.6% lower in the simulation compared to the experiment. Fig. 21. Illustration of through-thickness height h crack of the intralaminar crack and the delamination length e used to calculate Local Delamination Ratio in Eq. (30) as defined in [24]. ©2014 Elsevier Ltd. 7 If the damping is turned of abruptly when the damping time is over, the slenderness of the beam will yield vibrations.

Load-displacement response
In Fig. 22 we compare the load-displacement response of the reference (non-adaptive) simulation to the experimental counterpart. Up to the point of delamination propagation at Δ z = 2.2 mm the simulation matches the experiment. Thereafter, the response is softer. Since fibre failure is not included in the simulation we cannot compare the failure load, but at the experimental failure displacement Δ z = 6.6 mm the simulation load is approximately 11% lower and the total work is 7.6% lower. However, we want to emphasise that stiffening effects such as matrix bonding at transverse cracks and fibre bridging, which were observed in the experiment, are not included in the simulation.  [24]. The weak refinement is stabilised with a damping of C d,ref = 56 rad/ms linearly degraded over Δt ref = 0.2 ms. At the experimental failure displacement (Δ z = 6.6 mm), the load is approximately 9% lower and the total work is 3.9% lower in the adaptive simulation compared to the experiment.  [24]. No stabilisation of the refinements are made. At the experimental failure displacement (Δ z = 6.6 mm), the load is approximately 11% lower and the total work is 4.3% lower in the adaptive simulation compared to the experiment.
Furthermore, after delamination propagation, high-frequency vibrations can be observed in the simulation response. These are a side effect of the element-wise propagation of the delaminations, and thus hard to avoid without decreasing the in-plane element size (or decreasing the interlaminar strengths). Besides this, the dynamic effects in the simulation are small, and the kinetic energy is found to be three orders of magnitude lower than the total energy.
The load-displacement response of the non-stabilised adaptive simulation is plotted in Fig. 23, where we can observe oscillations after the weak refinement is activated at Δ z ≈ 1.0 mm. Initially, the magnitude of the oscillations are comparable to the vibrations  [24]. The weak refinement is stabilised with a damping of C d,ref = 56 rad/ms linearly degraded over Δt ref = 0.2 ms and the strong refinement is gradually separated over 100 time steps (Δt CZ = 0.0024 ms). At the experimental failure displacement (Δ z = 6.6 mm), the load is approximately 9% lower and the total work is 3.9% lower in the adaptive simulation compared to the experiment. associated to the delamination propagation in the reference simulations. However, the magnitude increases at the end. This type of non-physical oscillations, which are related to the refinement and not the dynamic problem itself, was one of the topics in our previous paper [6] and they can disrupt the solution, causing unwanted damage propagation.
The result of stabilising the weak refinement is shown in Fig. 24, where we can observe strong qualitative similarities to the reference simulation in Fig. 22. The load response is a bit too high throughout the simulation, however, this is mostly related to a too stiff response of the unrefined beam. To reduce the initial incorrect stiffness, the beam must be modelled with more element through the thickness, which defeats the purpose of adaptive refinements.
When the strong refinement is activated in the adaptive simulation with weak refinement stabilisation (can be seen as a small disturbance at Δ z ≈ 1.3 mm in Fig. 24) there is an immediate (but small) damage initiation in many of the CZ elements. As we mentioned already in [6], this is a strong indicator that the refinement needs to be stabilised. By activating stabilisation of the strong refinement this non-physical CZ damage initiation can removed. As can be seen in Fig. 25, the effect of this is negligible on the load response, however, the refinement-induced damage initiation in the CZ elements is removed.

Damage evolution in the cohesive interface
In the experiment 15 matrix cracks were observed along the beam, with five of these found in the central load span (between the top load pins). Interestingly, the same intralaminar crack pattern was also formed in the current non-adaptive reference simulation (in total 15 matrix cracks out of which five are in the load span). These resulting matrix cracks are indicated as dashed lines in Fig. 26a, where also the CZ element damage at a displacement of Δ z = 4.0 mm is plotted. Furthermore, the minimum crack distance is 1.5 mm, which is three times the enforced crack-spacing distance. We can thereby conclude that the crack-spacing enforcement has not restricted the formation of matrix cracks in the simulation.
Compared to the non-adaptive reference simulation, we observe the same crack pattern and similar cohesive damage development also for the adaptive simulations with stabilisation (both only weak and strong cases), cf Fig. 26c and d. However, when no stabilisation is applied, only 14 matrix cracks form and the CZ element damage pattern appears very different with non-smooth evolution along the beam, as shown in Fig. 26b. The jagged shape in Fig. 26b indicates that much of the CZ damage is related to the coupling of the intralaminar damage to the reduction of the CZ critical fracture energy, described in Section 3.2. That is, there is non-physical intralaminar damage present when the weak refinement is not stabilised. In addition, compared to the reference simulation, there are more CZ elements with damage, even out towards the ends of the beam.
In Fig. 27 we have compared experimental and simulation LDR-curves (cf Eq. (30) for definition) for the five load-span matrix Since the delamination cracks in the simulations evolve element by element (when full CZ damage is reached) the simulation curves grow in a stepwise manner.
cracks. Overall, the simulations overestimate the delamination growth, which is in line with the softer load response of the simulations compared to the experiment. However, the adaptive simulations with stabilisation ( Fig. 27c and d) qualitatively compares well with the reference simulation (Fig. 27a), although with a slightly lower delamination growth prediction. In contrast, the adaptive simulation without stabilisation (Fig. 27b) predicts a very different and generally larger delamination growth compared to the other simulations.

Interaction between smeared ply crack and delamination
To evaluate the effect of the proposed coupling between intra-and interlaminar damage in Section 3.2, the non-adaptive reference simulation was rerun with this coupling deactivated.
The load-displacement response, plotted in Fig. 28, is similar to the reference response in Fig. 22. However, instead of the expected 15 matrix cracks along the beam, 25 cracks appear of which nine are in the load span (compared to five in the experiment and the reference). Furthermore, the CZ element damage, plotted in Fig. 29a, is underestimated compared to the reference case in Fig. 26a. The LDR-curves for the nine load-span matrix cracks are plotted in Fig. 29b, which shows that no delamination cracks form in the load-span for the duration of the simulation.

Fig. 28.
Load-displacement curve for the non-adaptive simulation without coupling between intra-and interlaminar cracks (cf Section 3.2), compared to the experiment [24]. At the experimental failure displacement (Δ z = 6.6 mm), the load is approximately 11% lower and the total work is 6.7% lower in the simulation compared to the experiment.

Discussion
If we compare the computational time of the initially refined and the unrefined models up to a load displacement of Δ z = 1.0 mm, we find that the unrefined model is approximately 8% faster. While this is only a moderate computational save, we would like to emphasise that all (including the non-active) extra DOF used to define the adaptive refinements are updated by the LS-DYNA solver during the entire simulation. Thus, the improved computational efficiency is only related to the internal computations of the user element.
Section 4.2.3 shows that the coupling between intra-and interlaminar cracks, proposed in Section 3.2, does not affect the resulting load-displacement response to any large extent. However, when the coupling is deactivated the damage extent in the CZ elements is underpredicted whereas the number of matrix cracks is overestimated.
Finally, in Fig. 30, the load-displacement curves are shown for all simulations. To facilitate a comparison, all curves have been filtered with an SAE 20 kHz filter. It can be observed that all simulations compare quite well to the experimental curve, however, the adaptive simulation without stabilisation (red line with circular markers) show a different vibration pattern compared to the other simulations. This is particular evident towards the end of the simulation. We can therefore conclude that the simulations with adaptive refinement can model the experiment, however, the weak refinement may have to be stabilised. More importantly, however, it is clear that the adaptive simulations can achieve a similar prediction as the reference simulation.

Concluding remarks
In a sequel of two connected papers (Part I and II), we have addressed the challenge of developing accurate modelling tools that are efficient enough to allow full car crash FE simulations. In Part I [6], we presented the implementation of an adaptive refinement method for modelling of multiple and arbitrarily located delamination cracks using an ESL shell model in the LS-DYNA explicit solver, including necessary stabilisation techniques. The method enables the laminated structure to initially be represented by a single layer of solid shell elements through the thickness. Refinement indicators are then used to activate so-called weak refinements (weak discontinuities), which facilitates a better prediction of stress state and intralaminar damage, and strong refinements (strong discontinuities) needed for a proper representation of growing delamination cracks. In this second part of the work we completed the interlaminar-based refinement indicator, presented in [6], with indicators based on intralaminar criteria in order to model both intra-and interlaminar progressive failure and the interaction thereof.
In order to limit the computational costs, the proposed method does not model intralaminar failure with separate DOF. That is, only through-the-thickness and not in-plane refinements are made. Instead, we propose to use a smeared-crack model to represent intralaminar fracture (in this paper limited to matrix cracking). In particular, we discuss how the interaction between intra-and interlaminar failure can be considered by coupling the cohesive fracture energy to the intralaminar damage state of the adjacent layers. If this coupling is not present, the damage evolution of the CZ elements is restricted and the correct intra-and interlaminar crack propagation cannot be predicted. Moreover, we demonstrate the importance of considering both the element size and the crack Fig. 30. Load-displacement curves for all simulation filtered with an SAE 20 kHz filter, compared to experiment [24]. orientation when properly regularising such smeared-crack models.
In our numerical examples we first show that the presented mesh-based regularisation scheme can accurately compute the characteristic length of the cracked elements. In contrast to a more traditional element-based regularisation, for which the energy dissipation can be overpredicted by up to a factor of two. However, as showed in the second example, the choice of characteristic length type does not control the growth direction of the crack. When the crack-growth direction show clear mesh-orientation dependency, the element-based characteristic length actually yields better energy dissipation predictions. This is because the error in crack-path prediction is compensated by another error in the energy regularisation. Nevertheless, when the crack follows the correct crack path (either by chance or by forcing it) a mesh-based characteristic length must be used for correct energy dissipation prediction.
The intralaminar model we presented in Section 3.1 is intended to be used primarily for development purposes and we acknowledge that a more advanced model needs to be adopted in order to properly model the intralaminar failure of FRP. Thus, the next step is to use a state-of-the-art material model, e.g. similar to Maimí et al. [32], Gutkin and Pinho [33] or Costa et al. [21], in the future.
In the four-point beam bending example it is evident that the adopted solid shell kinematics (LS-DYNA element formulation 3) is not fully suitable to model bending problems. Therefore, in future work it should be investigated if alternative shell formulations further increase the accuracy (or the computational efficiency) of the proposed method such that bending problems can be modelled with one element through the thickness. This example also shows that computational efforts can be saved by employing adaptive refinements. However, as we mentioned already in Part I [6], we are currently not fully exploiting the computational efficiency of the adaptive refinement approach. To better achieve this the refinement should be limited to an in-plane patch which then can expand if necessary and a less expensive element formulation should be employed. By implementing these features the computational efficiency is expected to increase as larger parts of the structure can remain in an unrefined and computationally efficient state during the simulation.
Finally, we conclude that our adaptive refinement method, presented in Part I [6] and expanded in this paper, can reproduce similar results as a high-fidelity model while saving computational effort. This will help to enable computationally efficient crash simulations of laminated structures, which in the long run will help to develop automotive structures made of carbon FRP.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
and ζ is then These basis vectors are then used to define the transformation from reference global to current element frame ( 0 → ℰ t ), expressed in matrix form as In Eq. (A.1) it is assumed that the material and lamina frames rotate with the element such that Q and R can be defined in the reference configuration and be constant throughout the simulation. The material directions are defined in global coordinates on the input deck and an element to material transformation matrix can therefore be calculated as where Y 0 are the input material directions and A 0 is the global to element transformation matrix at time zero (i.e. x = X in Eq. (A.2) to (A.3)). We have chosen to define the material directions using LS-DYNA AOPT = 3 [34], which like the element coordinate system is constructed using the midsurface. This means that Q will be a rotation ω around the midsurface normal where β is the angle between the material first axis and the lamina fibre direction as illustrated in Fig. A.32. Each layer of the element need to have the angle β given as model input.

Appendix B. Element-based characteristic lengths
To estimate the element-based characteristic length, we make the following assumptions: • A crack pass through the centre point of the element; • The geometry of the element changes very little prior to fracture; • We can separate the length into one part representing cracks perpendicular to the midsurface and one to represent the inclination of the crack (this assumes that the element thickness h is fairly constant).
The definition of the element-based characteristic length in Eq. (21) can then be approximated as where Fig. A.32. The material coordinate system (a, b, c) is a rotation (ω) of the element coordinate system (ξ, η, ζ) around the midsurface normal (left) and similarly the lamina coordinate system (1, 2, 3) is a rotation (β) of the material coordinate system (right).  is the characteristic length for cracks perpendicular to the midsurface, calculated by relating the ply midsurface area A e to the crack length L c (cf Fig. B.33). To derive the expression in Eq. (B.1), we have used that V e ≈ A e h and A c ≈ L c h/|cos(ϕ)|, together with the fact that the characteristic length must be equal to the thickness h for crack angles of ϕ = 90 ∘ Fig. 5. It is noteworthy that the only difference between the element-based expression in Eq. (B.1) and the corresponding mesh-based in Eq. (20) is how the parts representing cracks perpendicular to the midsurface L □ e⊥ and L ⊞ e⊥ are calculated. As mention before, matrix cracks can form parallel to the fibre direction while fibre fracture will primarily occur perpendicular to the fibres. Thus, two corresponding element-based characteristic lengths L □ em and L □ ef need to be calculated for each unique ply orientation. Furthermore, for solid shell elements the ply midsurface areas A e,i can vary through the thickness. Thus, each element ply requires a separate set of characteristic lengths L □ em,i and L □ ef,i , however, we will drop the ply index i in the following. Using Eq. (B.1) the element-based characteristic lengths for a matrix and a fibre crack become and we have, once again, assumed that the fibre cracks always are perpendicular to the midsurface, i.e. ϕ f = 0. The angle ϕ m is a result of the loading of the element as described in Section 3.1. However, we can calculate L □ em⊥ and L □ ef⊥ based on the element geometry and the fibre direction. This calculation requires to estimate the ply midsurface area A e and the lengths of the potential matrix and fibre cracks, L cm and L cf , illustrated in Fig. B.33. Since we assume small geometrical changes of the element prior to fracture this is done using the reference nodal coordinates at the beginning of the simulation. We have in the following restricted our calculations to general hexahedral elements and do not allow for wedge elements to be used.
We first define the ply fibre and transverse direction vectors E 1 and E 2 , following the procedure in Appendix A. A potential matrix crack will be parallel to the fibre and thus along the E 1 direction. Similarly, a fibre crack will be perpendicular to the fibre, i.e. along the E 2 direction. This is illustrated in where P l and P ij are the points with the closest distance between the lamina direction vectors and the edge vectors (which are likely skew): Finally, the area of the midsurface is calculated using the diagonal vectors of the midsurface where × denotes cross product.