A special ﬁnite element method applied to off-axis tunnel cracking in laminates

Off-axis oriented tunnel cracking in a laminated structure is modeled using a special 2D off-axis ﬁnite element formulation, thus replacing full 3D ﬁnite element simulations with much less demanding 2D simulations. The 2D off-axis element is formulated as a user element in the commercial ﬁnite element code ABAQUS, and the user element code is made available for the reader. The ﬁnite element formulation is used to predict the steady-state energy release rate and the mode-mixity for either one isolated crack or multiple interacting tunnel cracks in the central ply of a [0 /θ/ 0 / − θ ] s layup. The laminate is loaded uni-axially or bi-axially, and analyses are presented for glass and carbon ﬁber composites for arbitrary layup angles θ . The 2D ﬁnite element model is found to give precise predictions when comparing the detailed simulations made by a full 3D ﬁnite element model. It is demonstrated that despite being a 2D method, the model rigorously provides detailed solutions for the out-of-plane deﬂection and strain. With a speed-up of more than four orders of magnitudes, it is thus possible to perform new, detailed, and very accurate studies of the dependency of laminate thickness as well as neighboring cracks on off-axis tunnel cracking.


Introduction
Tunnel cracks are often seen in layered materials where less ductile layers under tensile loading will develop transverse-oriented tunnel cracks when surrounded by more ductile materials, Reifsnider and Jamison (1982); Tong et al. (1997). During fatigue loading, tunnel cracks can act as damage initiation points for the subsequently fiber fatigue failure in the load carrying fibers, Talreja (1981); Jamison et al. (1984). Hence, understanding and controlling the initiation and growth of such tunnel cracks can potentially have a large influence on the subsequent fatigue damage evolution in the load carrying fibers and therefore on the overall stiffness degradation and fatigue life-time of the laminate.
A large number of studies have addressed 90 • transverse cracks where numerical studies typically rely on 2D plane strain finite element models. In this way, similarly to Beuth (1992) for channeling surface cracks, Ho and Suo (1993) used an energy balance method based on the crack opening profile of a fully developed tunnel crack and the stress-state of the non-cracked case to obtain steady-state energy release rates. Hutchinson and  extended this 2D plane strain analysis to crack interaction studies. Lundmark and Varna (2011), used the 2D plane strain model to investigate the stiffness reduction due to the tunnel-crack density while Carraro et al. (2019) investigated tunnelcrack-induced delamination. Various other investigations of tunneling cracking at 90 • in laminates exist, such as Cox and Marshall (1996) discussing several aspects of the steadystate tunneling cracks, or ; , treating both tunneling cracks and delamination caused by these cracks. These investigations are typically conducted with generalized plane strain models and thus limited to cracks at 90 • .
For orthotropic materials in the presence of weak interfaces, tunnel cracks can develop in other directions than the transverse 90 • direction. An important case is multi-oriented fiber-reinforced polymer plies in composite laminates where the tunnel cracking typically will follow the weak interfaces between the fiber and matrix material. For non-90 • tunnel cracks, the crack-tip will be loaded in a mixed-mode (a combined normal and shear loading). Quaresimin et al. (2014) proposed a generalisation of the energy balance method for predicting the energy release rate and mode-mixity, which was applied by Maragoni et al. to determine the energy release rate using a 3D finite element model. A prediction based on a full 3D model, Mikkelsen et al. (2020) found that the method was giving accurate predictions of both the energy release rate and the mode-mixity even when based on a 3D finite element analysis of a fully developed crack without explicitly modeling the crack-tip stress-field. The numerical solutions showed that the far-field self-similar solutions (far behind the crack front) could be used to determine the crack-tip load based on an energy balance, thus essentially utilizing 2D features of the solution which may be directly modeled using specialized off-axis 2D finite elements.
In the present study, such an off-axis finite element is developed and implemented in the commercial finite element code ABAQUS. The element is used for predicting the energy release rate and mode-mixity for a specific symmetric laminate layup [0/θ/0/ − θ] s with a central tunnel-crack in the symmetry ply. The material is assumed to be a linear elastic orthotropic material representing typical glass and carbon fiber composites. Even though only such specific cases are considered, the off-axis finite element framework is relevant for modeling a wide class of problems, including material non-linearity in the form of plasticity as well as micro-mechanical effects such as inhomogeneous stress states between the reinforcing fibers, and effects of curing induced residual stresses -all potentially in combination with crack growth simulations based on 2D energy balances.
The new 2D element formulation has been validated against full 3D simulations by comparison of both the 3D stress-state in a symmetric laminate, the energy release rate and the mode-mixity for tunnel-cracking under uniaxial loading, Mikkelsen et al. (2020), as well as under bi-axial loading, Bangaru et al. (2021). The user element can be downloaded from Mikkelsen et al. (2021). Fig. 1 shows two different 3D models which can be used for calculating the average steady-state energy release rate,Ḡ ss , and the average mode-mixity,ψ, for a tunnel crack front positioned in the central layer of the composite layup [0/θ/0/ − θ] s . Thus, except for the ply thickness studies presented in Fig. 9, the middle layer has the thickness 2h whereas all other layers have the thickness h with a resulting total laminate thickness 2H = 8h. All laminates considered, are then symmetric and balanced. A Cartesian coordinate system, x i , is aligned with the sides of the laminate. The laminate shown in Fig. 1 is to be loaded in uni-axial tension in the x 1 -direction by the average membrane stress N 1 , but bi-axial loading with N 3 in the x 3 will also be investigated. The crack front is loaded in mixed Mode I and Mode II. A standard material coordinate system aligned with material directions is also introduced as (x L , x T , x T ). Here, x T will always be parallel to x 2 and the angle between x 1 and x L is θ. Fig. 1a shows the full cracked laminate. The tunnel crack front in the middle layer is positioned at an arbitrary position in between but far from the laminate edges. The crack front shape will generally vary with x 2 , but under steady-state conditions, the average

Crack front analysis:
J(x 2 ); K I (x 2 ); K II (x 2 ) (a) Full-thickness model including the crack front.
Energy balance method: x 1 x 3 x 2 x T Half-thickness model without the crack front. Figure 1: Two 3D tunnel crack models with definitions of load, geometry and two coordinate systems (x 1 , x 2 , x 3 ) as well as (x L , x T , x T ).
steady-state energy release rate,Ḡ ss , along the crack front can be found as the average value of the J-integralḠ and the corresponding average mode-mixity,ψ, can be defined as Here, K I (x 2 ) and K II (x 2 ) are the Mode I and Mode II stress intensity factors, respectively.
It is noted that K I (x 2 ) and K II (x 2 ) as well as the J-integral will depend on x 2 in general.
As the layup considered is symmetric and balanced, only half of it needs to be analyzed. Furthermore, in Mikkelsen et al. (2020), it was validated by 3D analyses, that the average steady-state energy release rate,Ḡ ss , and the average mode-mixity,ψ, along the crack front can be found from an off-axis generalization (see Quaresimin et al. (2014)) of the energy balancing method described in Beuth (1992); Ho and Suo (1993). Hence, the crack front itself does not need to be represented in the model in order to evaluate the average energy release rate,Ḡ ss , and mode-mixity,ψ. This is illustrated in Fig. 1b where a full-width crack without any crack front is shown for the half-thickness laminate. Then the average energy release rate,Ḡ ss , and mode-mixity,ψ, will be taken as a Mode I and a Mode II contributionḠ with (Bao et al., 1992) Suo et al., 1991). Exploiting symmetry in Fig. 1b, the average energy release rates in Mode I and Mode II, respectively, are calculated as where δ n (x 2 ) and δ t (x 2 ) denote the normal and tangential crack opening profiles, respectively, of the full-width crack (see Fig. 2b), while σ T T (x 2 ) and τ LT (x 2 ) are the transverse normal and shear stresses, respectively, before initiation of the tunnel crack. It is noted, that since the crack considered is fully aligned with the longitudinal material direction denoted by subscript L the normal Mode I opening denoted by subscript n follows the material transverse direction denoted by subscript T and tangential Mode II opening denoted by subscript t follows the material longitudinal direction denoted by subscript L.
The stress-state used to evaluate energy release rate is that existing in the absence of the specific tunnel crack analyzed. For a tunnel crack growing between neighboring cracks, this stress-field is found from a finite element solution, while for isolated tunnel cracks, the constant far-field stress in the off-axis ply is used. In this latter case, Eq. (5) can be written asḠ where the average crack openings can be calculated from as done in Mikkelsen et al. (2020).

2D formulation
Neglecting edge effects, the blue box in Fig. 2a indicates a domain in which the deformation state does not vary along the direction of the crack. Therefore, instead of finding the crack opening shape using the computationally expensive solution strategy in The solution will be obtained on a 2D domain indicated by thick lines in Fig. 2a. In x 3 x 2 (a) Half-thickness model with an indication of the 2D model.    The computational 2D domain is initially oriented along the main tensile direction x 1 , and hence, it is oriented at an angle θ relative to the tunnel crack. The size of the domain is given by the length in the x 1 -direction, L and the half laminate height, H, as shown in Fig. 2c. Note that when referring to the x i coordinate system, the displacement field within the 2D domain, u i , will be three-dimensional except for specific values of θ, for which the solution will reduce to a planar displacement field.
The external far-field loading is prescribed as an elongation, ∆, in the x 1 -direction along the edge x 1 = −L as illustrated in Fig. 2c for x 3 = 0. The in-plane symmetries are imposed along x 1 = 0 for x 2 ≤ −h by u 1 = 0 and along x 2 = 0 by u 2 = 0 (see Fig. 2c).
In addition, anti-symmetric boundary conditions in the out-of-plane direction need to be prescribed by zero displacements in the x 3 -direction, u 3 = 0, along x 1 = 0 for x 2 ≤ −h and along x 1 = −L. This is illustrated by red crosses in Fig. 2c.
To accommodate such 3D deformations, a new coordinate axis perpendicular to the crack plane is defined as as also shown in Fig. 2b. For θ = 90 • one has x 1 = x 1 where the formulation will simplify to a generalized plane strain formulation. For θ ∈]0 • ; 90 • [, the unit of the x 1 -axis is stretched relative to that of the x 1 -axis. For θ = 0 • , eq. (8) is not defined . The total 3D displacement field, u i , only depends on x 1 and x 2 and will here be expressed by the sum of that originating from a homogeneous far-field displacement gradient,ū i, j , and a fluctuation field,ũ i , that vanishes in the far-field The far-field small strain and rotations, resulting from the imposed far-field loading, are given asε Without loss of generality, and consistent with the boundary conditions discussed above, we takeω i j = 0, implying thatū i, j =ū j,i . The gradients of the total displacement field, u i , are then given by

Element formulation
A special plane element with three translational displacements degrees of freedom, u i , in each of the M nodes is developed. The element is derived in a general form to model the full six components of a 3D strain field. The finite element solution will be based on the displacement field on the plane x 3 = 0 in Fig. 2, at which From Eq. (8) it follows that for x 3 = 0 one has x 1 = x 1 and u 3 =ũ 3 . For the finite element implementation, the total in-plane displacements, u α [x 1 , x 2 ], and the fluctuation of the out-of-plane displacement,ũ 3 [x 1 , x 2 ], are defined as the unknown fields. Following the decomposition of the displacement field, Eq. (9), and using Eq. (11) the small strain components are expressed by these unknown field variables as ε 11 =ε 11 +ũ 1,1 = u 1,1 ε 22 =ε 22 +ũ 2,2 = u 2,2 γ 12 = 2 ε 12 + 1 2 (ũ 1,2 +ũ 2,1 ) = u 1,2 + u 2,1 Here, the engineering shear strains are introduced as twice the tensorial shear strains, i.e.
γ 12 = 2ε 12 , γ 13 = 2ε 13 and γ 23 = 2ε 23 . Using Voigt notation, the strain vector, ε, can be found based on the element displacement vector, D, for an element with M nodes, as with the element strain-displacement matrix Here, N i for i = 1 . . . M are the shape functions and c = cot(θ) for short. It is noted, that B the element stiffness matrix, K, can then be established in the usual manner as where t is the unit thickness of the element, while E is the 6×6 constitutive matrix determined by Young's moduli, shear moduli and Poisson's ratios, Here, and, T 1 (θ) and T 2 (θ) are the transformation matrices between the material coordinate system (x L , x T , x T ) and x i , which depend on θ in each individual layer. The transformation matrices are given as where c = cos(θ) and s = sin(θ). The constitutive matrix defines the relation between stresses, σ, and strains, ε, through It is noted, that for θ = 90 • a generalized plane strain element is recovered where the out-of-plane strain ε 33 will have a constant value equal toε 33 . Beuth (1992); Ho and Suo (1993) based their analysis on a plane strain model with ε 33 ≡ 0.

User element implementation
The developed element is implemented in ABAQUS through a user element (UEL) subroutine. The results can be visualized conveniently using shadow elements based on a user material (UMAT) subroutine (see Mikkelsen (2007)). The preprocessing is based on a Python script, which creates the geometry, model properties, the load, and displacement boundary conditions. Then an input file for ABAQUS is generated and automatically modified to include the user element subroutine. The input file is submitted and finally, the Python script continues with the post-processing to compute the crack profile and the energy release rate. The Fortran code for the UEL and the UMAT based "shadow element" as well as the Python scripts used for the simulations, can be obtained from Mikkelsen et al. (2021).

User element subroutine
The UEL subroutine is called by each element and requires the definition of the stiffness matrix K and the residual vector. The stiffness matrix K can be found by evaluating the integral (17) using Gauß integration. The residual vector is then found by multiplying the stiffness matrix with the displacement and far-field strain vector D. Finally, the stresses and strains at the Gauß points are calculated with Eqs. (14) and (20).

Strain prescription to all nodes
The presented element formulation poses implementational challenges in commercial codes, due to the additional and unconventional degrees of freedom regarding the far field strains,ε. In the formulation of the special element strain-displacement-matrix, B in Eq. (15), the strains,ε, of the homogeneous far-field are to be applied as additional degrees of freedom common to all the elements.
In the case studied here of the symmetric and balanced laminate under in-plane loading in the (x 1 , x 3 )-plane the implementation can be simplified. Here, the only non-zero strains of the homogeneous far field areε 11 andε 33 . As the tunneling cracks inside the laminates have a negligible influence on the far-field, the homogeneous far-field strains can be approximated by classical laminate theory and the corresponding additional farfield degrees of freedom can be eliminated. This step is done to simplify the element implementation. Note, that adding the additional degrees of freedom concerning the farfield in the element formulation instead of estimating the far field beforehand will yield identical results.
Thus the far-field average strains are calculated during pre-processing with classical laminate theory, and are thus known. The strain-displacement matrix can then be included in the UEL subroutine in the following form Furthermore, an element displacement vector is defined with size (3M + 2) × 1, where the last two rows represent the far-field homogeneous strain field. This is then used to calculate the stiffness matrix, K in Eq. (17)

Shadow elements
The visualization of the results of a UEL subroutine requires an additional step. By default, ABAQUS is only able to illustrate the nodal displacements of a UEL but is not able to interpret the field values inside the element, as neither the interpolation nor the location of the Gauß points is known by the post-processing tool in ABAQUS.
A solution to this is using similar standard elements from the ABAQUS library -here referred to as shadow elements (see Mikkelsen (2007)). The standard elements must have the same Gauß point locations and numbering as the element developed. A secondary identical mesh has to be created using the standard elements. The standard elements use the same nodes as the UEL, so that the nodal results in terms of displacements and reaction forces, can be transferred. This allows for a visualization of the displacement field.
The stresses and strains are transferred with a UMAT subroutine through a Fortran common-block from the UEL subroutine. In the UMAT subroutine, a zero stiffness constitutive model leads to a zero stiffness matrix. This ensures, that the material subrou-

Input file
The input file defines the model that is to be solved. It is generated by ABAQUS after preprocessing and can then be submitted for solving. The inclusion of the UEL requires a modification of the input file, as the UEL has to be defined and assigned to the individual elements, together with the specific UEL inputs in terms of far-field strains. Subsequently, copies of all UEL have to be created in relation to the specified shadow elements. These operations are built into a Python script.

Python script
The Python script automates the entire process and can be run through ABAQUS to calculate the energy release rate and crack profiles. It includes the pre-processing to create the model and calculation of the homogeneous far-field strains from laminate theory. Furthermore, the generation and modification of the input file, solving, and post-processing for obtaining the energy release rate are included. The script requires the geometrical dimensions, material properties, and loading. This enables an efficient method for carrying out systematic parametric studies. Material Three glass fiber composites (GlassFRP1, GlassFRP2, and GlassFRP3) and two carbon fiber composites (CarbonFRP1 and CarbonFRP2) will be considered in the analyses of a The out-of-plane element coordinate direction, x 3 , in the 2D off-axis finite element formulation is aligned with the orientation of the center off-axis ply. Thereby, the stress components in the finite element solution will be given directly in these off-axis material orientations. Fig. 3 shows the stress components in the normal direction to the fiber orientation, σ T T , and the shear stress, τ LT , for cases without any cracks. For the cases with cracks to be studied, these two stress components are required in Eqs. (5) and (6) for calculating the average steady-state energy release rate and mode-mixity for a tunnel crack in this layer using the energy balancing method. For validating the ABAQUS UEL implementation of the off-axis 2D finite element formulation, the predicted stress components have been compared with a Classical Laminate Theory (CLT) in Fig. 3, and a very good agreement for all the cases analyzed is found.

Results
The material systems studied here (see Tab. 1) were also considered by Quaresimin et al. (2014); Mikkelsen et al. (2020) and Bangaru et al. (2021) in 3D. Hence, direct validation of the novel 2D formulation against the full 3D model is possible. In the following, this will be done for a number of cases. The normalized far-field out-of-plane strain found using the Classical Laminate Theory isε 33 /ε 11 = −0.452 . Fig. 4 shows the out-of-plane displacement, and the out-of-plane normal strain com-ponent found from the 2D finite element solution near the tunnel crack, for the material system GlassFRP1 with the layup angle θ = 45 • . Fig. 4a shows the full 2D model with contours of normalized out-of-plane displacement, u 3 /ε 11 , Fig. 4b shows a closeup view with the mesh of Fig. 4a near the crack, and Fig. 4c shows the same close-up view, but with contours of normalized out-of-plane strain, ε 33 /ε 11 . As the simulations are based on a linear finite element model, a scaling factor is chosen for the in-plane deformation field in order to have a visible crack opening in Fig. 4. Despite being a 2D finite element model, the off-axis finite element formulation is able to correctly predict the out-of-plane displacements, u 3 , which are included as additional degrees-of-freedom in the finite element formulation, as well as the out-of-plane strain ε 33 =ε 33 − cot(θ) ·ũ 3,1 (see Eq. (13)), whereε 33 (andε 11 ) is found from Classical Laminate Theory. In this specific case,ε 33 /ε 11 = −0.452 . Both the displacement field and the strain field exhibit a localized deformation around the tunnel crack while they represent homogeneous deformation in the major part of the model. The solution obtained has been found to agree well with the full 3D finite element prediction with a maximum out-of-plane deflection of between the strain field from the 2D and 3D solution along the crack front (the top-ply).

Uniaxially loaded off-axis tunnel-crack
This is due to that the strain values from the finite element solution in the finite element integration points are extrapolated out to the crack surface where the element size in the 2D mesh, L el × h el = 0.004h × 0.004h, is 1/10 th of the element length used in the 3D mesh and therefore give a much more accurate extrapolation in the 2D model. The jump in the strain prediction between the second and third ply in the 3D solution is due to a use of a much coarser mesh in the two lower plies joined to the two upper plies with a tie constrain between the two incompatible meshes. Fig. 6 shows the detailed crack opening profile given by the normal opening, δ n , and transverse opening, δ t , using the definition given in Fig. 2b, which is found by a projection onto the crack normal in the (x 1 ; x 3 )-plane. In the contour plot of Fig. 4  maximum normalized crack opening displacement was found to be (u max 1 ; u max 3 )/(2hε 11 ) = (−0.691; −0.242). This corresponds to a maximum crack opening of (δ max n ; δ max t )/(2hε 11 ) = (0.635; 1.320), which coincides with the end-points of the curves in Fig. 6. In Fig. 6, the crack opening profile obtained from the 2D analysis is also compared to the profile found from the full 3D solution given in Mikkelsen et al. (2020). The lines in Fig. 6   (7) based on a full 3D simulation and the new 2D off-axis finite element simulation are compared. The difference is found to be less than 1.5%. Using Eq. (6), the average energy release rate in Mode I and Mode II can be calculated from the average crack opening and the stress-field in the un-cracked laminate. These values are also shown in Tab. 2 together with the average steady-state energy release rate,Ḡ ss , and average mode-mixity angle,ψ, calculated using Eq. (3). Very good agreement has been found between all the values with less than 0.5% difference between the 2D and 3D results for the average steady-state energy release rateḠ ss and the average mode-mixityψ. calculations. An excellent agreement is obtained, whilst saving more than four orders of magnitude in computational time using the novel 2D off-axis element.

Biaxial loaded off-axis tunnel crack
An off-axis tunnel crack loaded biaxially was analyzed by Bangaru et al. (2021). The layup analyzed here is the same and the material system studied here is GlassFRP3 with igure 7: The total average energy release rate,Ḡ ss , and average mode-mixity,ψ, for the four material cases. The lines are from the 2D off-axis finite element, and the circular and cross marks are from the 3D simulations in Mikkelsen et al. (2020). the properties shown in Tab. 1. The biaxial loading is given by the two loads N 1 and N 3 , and the effect of the biaxial loading ratio N 3 /N 1 is investigated. Bangaru et al. (2021) performing such an analysis based on a full 3D model applying a superposition of the average stress intensity factors K I and K II . Using the 2D off-axis finite element formulation, those corresponding simulations are sufficiently fast that all the load cases considered here are simulated directly applying the different load ratios N 3 /N 1 . The out-of-plane loading is applied in the 2D off-axis finite element model through theε 33 value which together withε 11 is determined from the classical laminate theory. Note that the applied axial deformation ∆ along the outer boundary should correspond toε 11 calculated from Classical Laminate Theory.  Fig. 8 shows the average steady-state energy release rate,Ḡ ss , and the average modemixity,ψ, predicted for biaxial load ratios in the range N 3 /N 1 ∈ [0; 1] for three different layups with θ = 30 • , θ = 45 • and θ = 60 • . The load ratio N 3 /N 1 = 0 corresponds to the uniaxial loaded off-axis tunnel crack analyzed in the previous sub-section. Similarly to Fig. 7, the curves in Fig. 8 are the predictions using the 2D off-axis finite element, while the circular and cross marks are predictions taken from the full 3D finite element simulations performed in Bangaru et al. (2021). As for the uni-axial load case shown in Fig. 7, a good agreement is also obtained using the 2D off-axis finite element method.
Comparing the predictions from the three different layup angles in Fig. 8, the blue curves for the smallest layup angle, θ = 30 • , show the largest sensitivity to a change in the biaxial loading ratio on both the energy release rate and the mode-mixity. For the uniaxial loading case, corresponding to N 3 /N 1 ≡ 0, this layup angle has the lowest predicted energy release rate together with a high mode-mixity angle corresponding to a Mode II (ψ 85 • ) dominated loading. This implies that a tunnel-crack is resilient to growth as the required energy release rate normally increases significantly when approaching Mode II crack tip loading, see e.g. Sørensen and Jacobsen (2009). Nevertheless, applying an additional transverse load, N 3 , to the laminate, both a significant increased energy release rate,Ḡ ss , and a significant decreased mode-mixity,ψ, approaching pure Mode I (ψ 0) loading for N 3 /N 1 ≈ 0.45, is achieved. For a further increased transverse load, the modemixity becomes negative according to Fig. 8, implying that K II in Eq.
(2) is negative as K I will always be positive. This effect follows from the change in direction of the shear stresses in the cracked laminate introduced by the N 3 loading as discussed by Bangaru et al. (2021). The change of sign indicates a shift in the direction of the shear loading between the crack surfaces going from a shear stress dominated by the N 1 -loading to a shear stress dominated by the N 3 -loading. In the numerical analysis, both the shear stress τ LT and the shear crack opening δ t will change sign while the Mode II contribution to the energy release rate G II will remain positive. As the mode-mixity in the present work is calculated from the inherently positive energy release rates using Eq.
(3), this change of sign must be defined differently, and hence the mode-mixity shown in Fig. 8 is defined where the sign ofψ is defined by the sign of τ LT similarly to Bangaru et al. (2021).
For a small θ-interval, τ LT andδ t have been observed to have opposite signs resulting in a negativeḠ II contribution when using equation (6). Obviously, a negative energy release rate is non-physical and will also result in an undefined mode-mixity value-based in Eq. can not be established rigorously, but was found to give a good approximation of the mode-mixity obtained using an accurate crack-tip analysis.

Thickness dependency of the off-axis tunnel-crack
The effect of the ply thickness has been investigated by Crossman and Wang for tunneling cracking at 90 • using generalized plane strain models. Also, in Mikkelsen et al.
(2020), a limited study of the influence of the ply thickness on the energy release rate was performed under uniaxial loading with arbitrarily oriented cracks. However, due to the increasing numerical demands in the full 3D model, especially for an increasing H/h aspect ratio, only a few cases were analyzed. Now, using the novel 2D off-axis finite element formulation, an extensive parameter study can be performed.  Mikkelsen et al. (2020) are also included for comparison in Fig. 9 as blue cross-marks. A good agreement between the full 3D and the new 2D finite element predictions can be observed. Despite the fact that the thickness dependence is found to be rather limited for the glass fiber cases, the carbon fiber case shows a significant variation for the layup angles θ = 30 • and θ = 45 • with a reduction of the steady-state energy release rate of about 80-90% when comparing very thin (h H) off-axis plies with equally thick plies. Both those two layup angles are cases where the crack-tip loading is pure Mode II with compressive normal loading over the crack surfaces. The decreasing energy release rates for thin off-axis plies with small offaxis angles is of high technical relevance when designing backing layers in non-crimp fabrics-based composites as tunnel-cracks in the off-axis backing bundles have been observed to be fatigue damage initiation sites in high cycle tension-tension cyclic loading (see Jespersen and Mikkelsen (2017); Wang et al. (2018)). Fig. 9b shows the corresponding change in the mode-mixity for the cases analyzed. All the cases show a reduction of the mode-mixity angle when reducing the off-axis ply thickness. A reduction, which is most significant for the GlassFRP1 case with a predicted mode-mixity decrease of up to 15 • for thin plies in a θ = 30 • off-axis layup. As a decreasing mode-mixity typical will decrease the critical energy release rate for the material (see e.g. Sørensen and Jacobsen (2009) required energy release rate for crack growth will decrease. Hence, whether the decreasing energy release rate, which has been predicted for decreasing ply thicknesses, actually will result in a lower tendency of the laminate to form tunnel cracks will depend on how sensitive those cracks are on the corresponding decrease in the mode-mixity.

Effect of neighboring cracks
For all the analyses performed so far, the tunnel-crack is assumed to be isolated from neighboring cracks. Typically, this will only be the case during the initial phase of tunnel cracking. Later, when multiple tunnel-cracks are present, they will shield the stress-field in a neighboring region, thus lowering the energy available for nearby tunnel-cracking. Hutchinson and Suo (1992)  (2020)). Due to symmetry, the cases studied here pertain to two pre-existing cracks and a tunnel crack growing in between at a distance, , to the existing cracks resulting in the crack spacing in the 2D off-axis finite element model given as / sin θ (see insert in Fig. 10). The presence of the neighboring cracks will result in an inhomogeneous stressfield in between the pre-existing cracks so that the stress-state cannot (as for the case of a single crack) be found from the Classical Laminate Theory. Hence, Eq. (6) cannot be used, but rather Eq. (5) with a through-thickness integration must be used. In order to perform the integration in Eq. (5), both the stress variation through the thickness before a crack has developed and the crack opening profile after the crack is fully developed is needed. This is obtained by solving the finite element model twice with two different boundary conditions along the crack-front at x 1 = 0. One model with u 1 ≡ 0 represents a configuration with no central crack for determining the stress-field, and one model with u 1 unconstrained represents a fully developed central crack for determining the crack opening. Therefore, each point on the curves shown in Fig. 10 is based on the results from two numerical solutions. Nevertheless, as each of these 2D solutions can be obtained in about one minute using a single CPU, all the curves can be established during a few hours of simulation time. Fig. 10 shows the dependency of the crack spacing, , on the available energy release rate and mode-mixity. The energy release rate in Fig. 10a is normalized by the value obtained for an isolated crack,Ḡ =∞ ss , while the crack spacing is given by the ratio h/ . For the GlassFRP1 case with θ = 45 • , the 3D finite element predictions obtained in Mikkelsen et al. (2020) are included as cross-marks. A reasonable agreement has been obtained where the difference for small crack-spacing is caused by the coarser mesh used in the 3D model. Regarding the energy release rate shown in Fig. 10a, a rather similar behavior is seen for all the analyzed cases where an influence from the neighboring crack is observed for crack spacing < 10h with a lowering of the available energy release rate of about 10% for a crack spacing of h/ ∈ [0.2; 0.35]. Similarly, the mode-mixity in Fig. 10b is seen to be significantly influenced for crack spacings less than about 10 ply thicknesses ( < 10h). With closely spaced neighboring cracks, the mode-mixity is seen to increase for all the cases analyzed, with the largest increase of almost 20 • observed for the material case GlassFRP1 for the layup θ = 60 • . The results quantify a decreasing tendency for tunnel cracks to appear in between pre-existing tunnel cracks based on both the increasing mode-mixity and the decreasing energy release rate with h/ .

Conclusion
In this study, it has been demonstrated for a symmetric and balanced laminate, that it is possible to obtain accurate predictions of the 3D deformation and strain-fields around a tunnel-crack in an arbitrarily oriented ply layup using a specialized off-axis 2D finite element model. The main advantage of the 2D formulation is the significantly reduced computational time of more than four orders of magnitude, enabling thorough parameter studies. The novel 2D model was implemented in the commercial finite element code ABAQUS using the user element interface. The accuracy was demonstrated for a number of different tunnel-crack cases. The code can be obtained from Mikkelsen et al. (2021) where a GitLab link to the newest version of the code can also be found. In addition to being able to reproduce results from 3D simulations, the computationally much more efficient method has enabled detailed studies of the dependence of the tunnel-cracking on both the ply thickness and crack distance.
In the present study, the new 2D model has been limited to analyzing the available energy release rate and mode-mixity for the growth of an arbitrarily oriented off-axis tunnel crack, using a linear elastic material model. However, the framework has a broad range of applicability. For instance, the modeling technique can be extended to include the effects of the randomness in the actual micro-structure on the stress and strain-fields by a detailed representation of the fiber reinforcement and voids as in Vajari et al. (2014).
The formulation could also be extended to cover non-linear cases, including elasto-plastic deformations in the off-axis ply as well as delamination predictions using a compatible off-axis cohesive element.