1 Introduction

Carbon-fibre-reinforced plastics (CFRPs) have been used in structures and also in electronics housings of satellites due to the weight efficiency achieved in extreme light-weight concepts [1, 2]. Naturally, composite-based housings can be applied to future 5G technologies with highly optimized attenuation windows in otherwise protected enclosures [3, 4].

In the case of electronic housings, the typical means to realize high enough protection against various radiation from the (space) environment, and to also control the thermal and electrical conductance of the CFRP-based parts, is to laminate metal foils as part of the CFRP lamination. Although these enclosures are not always primary load-carrying components, delamination of CFRP and metal foil would lead to significant deviation of the heat flux, radiation attenuation, and geometry. Therefore, predictive analyses of delamination are required in the design of CFRP-metal laminates for aerospace applications (Fig. 1) [5].

Fig. 1
figure 1

Three main design phases in the hybrid enclosure design for satellites [1, 5]

The main numerical methods for finite element delamination analyses are the Virtual Crack Closure Technique (VCCT) and Cohesive Zone Modelling (CZM). Both methods have their pros and cons in a practical design process since VCCT is primarily applied for structures with an initial flaw [6] and CZM requires fitting parameters other than pure fracture toughness and standardized fitting procedures do not exist. VCCT does not rely on energy dissipation by a cohesive zone after nodal release whereas CZM can be input various models of residual stiffness [7] before full release of the contact. Consequently, VCCT is typically applied for brittle crack propagation [8, 9]. The scientific challenge is to understand the limits of brittle crack propagation for VCCT so that dynamic effects remain insignificant.

The hybrid CFRP-metal laminates utilize steel or tungsten foils [10], which result in extremely high residual stresses during the manufacture [11, 12]. Due to the very high Young’s modulus of the unidirectional CFRP plies and, say, tungsten foils, even the smallest difference in thermal expansion leads to build-up of high interfacial loading between the dissimilar layers. Due to the fact that a separate crack onset phase is typically not simulated, any initial flaw tend to onset unstable crack propagation during a VCCT analysis. What follows is severe convergence problems in the numerical solution iteration or a need to apply artificial damping leading to an unclear error to the solution.

In this paper, we study the application of VCCT in the delamination analysis of CFRP-tungsten (CFRP-W) laminate. A cracked lap-shear (CLS) specimen is 3-D modelled using the finite element method and the effects of VCCT and CZM crack models are compared in terms of the crack onset stresses and the crack-tip loading. Additionally, a combined analysis with a separate crack nucleation model using a CZM-zone along with a VCCT zone is studied.

2 Methods

2.1 Reference Experimental Data

The experimental reference, i.e., the validation data, is based on a cracked lap shear (CLS) testing in this study. CLS testing is used for determining mixed mode fracture toughness values for interfaces and composites [13, 14]. The basic concept of the cracked lap shear specimen and the reference test setup in this study are illustrated in Fig. 2. The case material system in this study is based on a hybrid system of tungsten foil and CFRP. The detailed description of the laminate configuration and the test procedure can be found in a previous publication [5]. The fracture during the CLS testing of CFRP-W specimens is highly brittle and sudden. Several analytical methods have been developed to analyse CLS testing. However, accurate solutions of interfacial fracture toughness (Gcr) require taking into account 3-D residual stress distributions and related residual strain energy making numerical methods necessary [15].

Fig. 2
figure 2

Description of the reference test specimen and tensile test setup to be numerically simulated in this study. The cracked lap shear (CLS) specimen consists of a strap part and a lap part. The pre-crack tip is located at the interface between the strap and lap

2.2 Finite Element Method

A 3-D finite element (FE) model of the hybrid CLS specimen was generated by using ABAQUS® (standard, 2017, Dassault Systèmes). The strap and lap parts’ laminate thicknesses were 1.55 mm and 1.5 mm, respectively. The lay-up of unidirectional CFRP and tungsten (W) was [05/W] and [05] for the strap and lap, respectively. The total specimen length was 180 mm, width 20 mm, and lap length 115 mm. The FE model is shown in Fig. 3.

Fig. 3
figure 3

The finite element model of the CLS specimen, applied coordinate system, and CZM element zone of the combined method along the CFRP-W layer interface

The regions of the specimen that match with the test machine gripping point in reality, were subjected to boundary conditions (BCs). For all the BC surfaces, out-of-plane (Z-coordinate) displacements were fully set to zero (initial value). The thick end of the specimen (without lap opening) was fixed as regards the longitudinal displacement. In turn, the other end was subjected to enforced displacement to simulate real test machine loading.

All the parts of the geometric model were meshed using reduced integrated continuum elements (C3D8R) with a nominal size of 0.5 mm \(\times \) 0.5 mm. The solution procedure was divided into two separate steps to include thermal residual stresses. During the first step, a thermal load of 100 ℃ was applied over the model to simulate the cool-down phase of laminate cure. After the first step, the defined enforced displacement followed. CFRP and tungsten were modelled as linear material and the material constants are given in Table 1. In this study, three different methods were used to study the crack onset and propagation:

  1. 1.

    Cohesive zone modelling over the entire strap-lap interface;

  2. 2.

    Virtual crack closure technique applied over the entire strap-lap interface;

  3. 3.

    Combined method where CZM elements are applied at the interface edges and VCCT for the propagation over the rest of the specimen.

Table 1 Material data for CLS specimen modelling; for directional properties refer to Fig. 3

2.2.1 VCCT Method:

VCCT relies on the idea of virtual crack closure, where the force required to hold a nodal connection and the deformation due to a virtual crack opening are used to compute the energy release rate (ERR) related to the crack opening process. The computed ERR value is compared to a critical value to justify possible release of the nodal connection. In the event of release, the force-displacement behavior is presumed linear at the crack tip. The basic concept of VCCT is show in Fig. 4a.

Fig. 4
figure 4

a The VCCT method and the related nodal representation; b the combined method of this study

2.2.2 Combined VCCT-CZM Method:

CZM techniques for various fracture analyses are versatile and numerous different applications exists. Previously, it has been reported that CZM is needed to simulate the interfacial metal-CFRP debonding process in CFRP-W laminates under mixed mode fracture conditions [15]. Jokinen and Kanerva formulated a procedure for determining the critical fracture toughness value based on the crack onset during testing and, subsequently, to fit a CZM model with two separate critical traction levels based on crack propagation. This procedure is accurate and reliable, yet is rather element mesh-dependent due to the CZM zone that covers most of the fracture plane. Therefore, a combined method is considered in this study, where simultaneous application of CZM and VCCT is analysed. In this combined method, CZM is simulating the crack onset process (Fig. 2) and VCCT handles the crack propagation after the process zone of the real crack tip has fully developed.

The concept of the combined VCCT-CZM method is show in Fig. 4b. The length of the process zone, i.e. the modelled crack nucleation zone, is \(L = 1.0 \) mm in this study (CZM element size 0.5 mm \(\times \) 0.5 mm), to analyse crack nucleation modelling also at the specimen sides. For the CZM elements with zero thickness, the initial stiffness of the nodal bond was set to 1015 N/m3 for all three fracture modes. To analyse the applicability of the combined method for CLS testing and, particularly to CFRP-W interfaces, the fracture parameters are given the values verified in the current literature [15]. Likewise, the stress criterion for damage onset (CZM zone) was applied in a quadratic form as follows:

$$ f = \left( \frac{\tau_{1,a}}{\tau_{1}} \right)^{2} + \left( \frac{\tau_{2,a}}{\tau_{2}} \right)^{2} + \left( \frac{\tau_{3,a}}{\tau_{3}} \right)^{2} , $$
(1)

where \(\tau \) refers to traction, and sub indices (1, 2, 3) refer to the three fracture modes (mode I, mode II, mode III), and the sub index a to a momentary traction value. For energy release rate (ERR), the following linear power law was applied to account for mode interaction:

$$ f = \frac{G_{I}}{G_{Icr}} + \frac{G_{II}}{G_{IIcr}} + \frac{G_{III}}{G_{IIIcr}} , $$
(2)

where it was presumed \(G_{IIc}\) = \(G_{IIIc}\) in this study.

For comparisons, a CLS specimen model with pure CZM interface (0.2 mm elements) was computed according to the description in a previous work [15]. The selected crack models’ parameter values for all the three models are shown in Table 2.

Table 2 Input values for the different crack models of crack nucleation and propagation [15]

3 Results

3.1 Overall Simulation Response

The simulation of a brittle fracture is convenient for methods, which primarily rely on linear elastic fracture mechanics. For example, the CFRP-W hybrid laminate involves CFRP-W layer interfaces, where mode II dominated fracture has brittle response upon loading and the crack propagates in an unstable manner. Hence, the failure of the interface can be estimated to have linear response. The drawback of brittle failure is that the simulation of crack growth becomes computationally challenging due to convergence problems. Therefore, we pursue to compare different crack modelling methods to assess their capability to simulate unstable crack propagation. The overall strain-force ε-F response of the three crack models is shown in Fig. 5. It can be seen that the slope \(dF/d\varepsilon \) is essentially linear until the crack onsets and passes the strain measuring point. After passing the strain measuring point, \(dF/d\varepsilon \) increases significantly and the overall strain energy continues to grow in the specimen. This behavior is typical for unstable crack propagation where the crack is unable to extensively release strain energy and the surplus energy speeds up the crack propagation.

The behavior of the crack onset and propagation are clearly dependent on the residual stress state [15]. The omitting of thermal load and respective residual strains, the crack onset is postponed and the crack development into a delamination occurs clearly more intensively, i.e. the force range \({\Delta } F\) needed for the crack to pass the strain measuring point decreases (corresponding to \({\Delta } time\) in a real test).

Fig. 5
figure 5

Strain-force response of simulated CLS testing. Strains are recorded on the strap side of the specimen model and at the bonded section (average over 10×2 elements corresponding to a 5 mm strain gauge in reality)

To evaluate the performance of the crack models at the moment of crack onset, the crack tip loading in terms of the mode-mixity must be known for different crack lengths. The \(G_{I}\) and \(G_{II}\) distributions over the specimen width are shown in Fig. 6 for pre-crack lengths of 20, 30, 45 and 60 mm. It can be seen that the crack-tip loading is essentially constant over the 0...60 mm-range of crack length. Finally, near the test machine gripping, the crack-tip loading turns into mode II (dominated) shearing. The mode-mixity is given by \(\phi \):

$$ \phi = tan^{-1} \left( \frac{G_{II}}{G_{I}} \right) . $$
(3)

The mode mixity remains nearly constant over the specimen width and only at the free-edges, mostly due to very high residual stresses, turns into mode I dominated loading, as shown in Fig. 7. When residual strains are not accounted for, the dominance of mode II at the specimen free edges increases (Fig. 7b).

Fig. 6
figure 6

aGI distributions in the specimen width direction for four different pre-crack lengths; bGII distributions in the specimen width direction for four different crack lengths (computed using pure VCCT and ΔT = 86 ℃)

Fig. 7
figure 7

Mode-mixity distribution at crack tip: a for four different crack lengths (computed using pure VCCT and ΔT = 86 ℃); b for a case without any residual strains

3.2 Crack Onset Computation

The differences in the crack-tip stresses affect the formation of the 3-D crack-front, i.e. delamination shape after the onset of crack propagation. The initial phase of the crack onset represents the crack nucleation process at a real crack tip. Basically, the different crack models should produce similar (realistic) crack-front shape in addition to the valid load-strain response of the specimen or structure. It should be noted that the sheer crack-tip loading (mode mixity) cannot be directly compared between VCCT and CZM since the homogenized region (i.e. the simulated process zone) is different for the two methods. CZM is based on the concept of process zone whereas VCCT differentiates only between bond and full nodal release. Figure 8 shows the opening stresses per crack model and it can be seen that the CZM and VCCT-CZM models produce essentially similar stress-states around the crack tip; the Combined method (with \(L =\) 1 mm) is slightly stiffer due to the VCCT zone and, thus, the stress peak is slightly higher and shifted towards the crack tip. If the crack-tip stress-state is compared for a simulation case omitting residual strains (i.e. \({\Delta } T\) = 0 ℃), all the crack models produce exactly the same crack-tip stress-state, as expected.

Fig. 8
figure 8

Crack-opening stresses (σ33): a distributions when thermal residual strains are included; b distributions when thermal residual strains are omitted. Stresses are recorded at a tensile force level of F = 5 kN

The pure VCCT model faces severe convergence problems during the crack propagation simulation. The crack onset occurs during the residual stress step (simulating the cool-down after cure in reality) and the crack propagation at the specimen sides (free-edges) is extensive (see Fig. 9a). Finally, in the beginning of the tensile load step, the simulation either must be stabilized heavily or the tolerance of critical G values must be increased enormously to continue running the computation. In turn, the pure CZM crack model is capable of propagating the crack throughout the simulated test. The main difference between these two types of simulations is the loading (stresses) amidst the crack tip. The VCCT model does not allow any relief of strain energy until the full release of crack-tip nodes while the CZM elements deteriorate (see Fig. 9b) and the bond stiffness decreases after the damage onset criterion is satisfied (Eq. 1). By using a combined VCCT-CZM method, specific nucleation can be simulated (see Fig. 10).

Fig. 9
figure 9

Crack-front growth (delamination) for pure VCCT and CZM crack models: a VCCT delamination at ΔT = 86 ℃ before any mechanical load; b CZM element damaging (SCRT) after full thermal load (ΔT = 100 ℃) and axial load of 5 kN

Fig. 10
figure 10

Crack-front growth (delamination) for the combined VCCT-CZM model: a delamination after the thermal load step; b crack onset damage (SCRT) of the CZM elements after the thermal load step

To summarize, the VCCT model is validated (via its critical G values) for a tensile load level according to experiments, i.e. residual strain are included. When the delamination is allowed to propagate (set by operator for software) along the specimen edges (in addition to the pre-crack tip), the simulated crack propagation begins already during the residual stress step. This means that:

  • Simultaneous free-edge stresses and the pre-crack tip singular stresses are problematic for the VCCT model validated for the CFRP-tungsten laminate;

  • For the real crack, a high stress state without a singularity point does not lead to crack propagation (since it was not observed at edges);

  • For the real free-edge, already slight deformation (opening along the layer interfaces) release efficiently residual (thermal) strain energy.

The CZM model is validated (via its critical tractions) also by using the experimental data, in addition to its critical G values, which are based on the VCCT validation. Here, for the CFRP-tungsten laminate, the damage initiation criterion of the bi-linear traction-separation law prevents the ’pre-mature’ delamination, although the critical ERR level is achieved, because the stress-based criterion is not yet fulfilled (only at the very corners, Fig. 9). Also, the capability to relieve the stress peaking prior to full debond is shown in the lateral stress distributions and when compared to VCCT with hard nodal ties. VCCT simulation ends up to release the nodes and the delamination front halts only after a significant decrease of stresses (σ33 and \(\sigma _{32}\)), as shown in Fig. 11.

Fig. 11
figure 11

Lateral stress distributions when comparing pure CZM and VCCT simulation results at delamination (i.e. along the original pre-crack line): a opening stress (σ33); b shearing stress (σ32). All the stresses are recorded at ΔT = 86 ℃

4 Discussion

In the current literature, different means have been introduced to simultaneously model both crack nucleation and the crack propagation after full development of the crack. In general, the question of constant materials properties related to fracture (in most cases fracture toughness) depends on the length-scale of interest [16]. In this study, the process zone is negligible compared to other dimensions of the structure. For this type of simulation case, the development of the crack could be taken into account using a set of fracture toughness values (rising R-curve method) [17], coupled models of fracture and stress criteria [18], or alternatively by using the traditional CZM [19].

The fundamental challenge of CZM, which also offers possibilities, is the wide variety of parameters included in the damage onset criterion, traction-separation law, and in the criterion for mode mixity. The crack nucleation process in polymeric materials involves a complex set of micromechanical phenomena, such as crazing, cavitation, shear banding and mechanical interlocking [20,21,22], thus, it is expected that a crack nucleation model involves several material parameters. In our previous study, we proposed a method to fit a CZM model for the simulation of CFRP-W interfacial delamination for design purposes. In this study, we applied a combined method that uses critical G values for a VCCT zone according to an analysis of a non-propagating crack at the experimentally determined critical load level (Fc in a previous report [5]). The CZM zone of the combined method was given input parameter values according to the fitting procedure presented previously [15]. For the combined method, CZM elements (zero thickness) ahead of the VCCT zone are harnessed by a bi-linear traction-separation law and a power law is applied to account for ERR per fracture mode [23]. Originally, the critical tractions (τ) were validated based on local strain fluctuation at the strap and lap parts of the CLS specimen (due to crack-front propagation over the strain recording point).

The comparison of the three methods revealed that CZM and combined VCCT-CZM model can simulate the crack onset correctly. Pure VCCT model results in over-sized delamination due to severe propagation from the specimen edges towards the mid-line—the addition of external tensile load easily shuts down the simulation for reasonable values of control parameters. It should be noted that there are challenges to use the combined VCCT-CZM model to compute the entire delamination process through yet it is not necessary for analysing CLS testing. In the event that the simulated process zone is larger or in the order of the spatial size of stress gradients, CZM can dissipate strain energy and halt the crack in a realistic manner. In the combined method, the required process zone size is related also to the length of the transition zone between the methods (length L in Fig. 4). For the CFRP-W interface, the need for dissipation exists due to high strain energy induced by internal residual strains. The computational transition, after the crack reaches the VCCT zone in the combined method, tend to speed up the crack propagation because the mode II dominance and residual stresses remain essentially constant after each nodal release. At the VCCT zone, any degree of crack opening is due to the deformation of bulk material elements (CFRP/W) whereas the CZM zone allows crack opening due to the deformation of the interface elements (traction-separation model) and adjusting its momentary stiffness.

5 Conclusions

This study focused on the analysing of the application of pure VCCT, pure CZM, and a combined VCCT-CZM crack model in the simulation of highly brittle and high-energy intensive mode II dominated fracture. As a case study, hybrid CFRP-W radiation shielding laminate was simulated in a CLS test setup. It is important to note that the crack onset modelling using VCCT-CZM for CFRP-W was initially given parameter values based on the procedure defined by Jokinen and Kanerva [15] for a pure CZM-based fitted CLS model. In order to simulate the entire delamination process over wide complex shapes, the effects of CZM element size and the transition zone length must be studied in the future.

The simulation results were analysed from the point of view of crack onset stresses and the simulated crack nucleation process, namely delamination area in the FE model. The main outcomes of the simulation results are summarized as follows:

  • Pure VCCT crack model is computationally challenging to apply for the CFRP-W laminate with very high internal thermal residual strains;

  • Combined VCCT-CZM model can be used to simulate crack onset at the CFRP-W interface of the radiation protection laminate;

  • Crack onset modelling using VCCT-CZM for CFRP-W requires a minimum of \(L = \) 1 mm process zone (CZM elements + transition zone).