Experimental and numerical study of rate-dependent mode-I failure of a structural adhesive

ABSTRACT We present an experimental and numerical study of the rate dependence of the mode-I failure of adhesive joints, focussing on aluminium plates bonded with Araldite® 2015. For the experimental part, we tested 24 double-cantilever beams (DCB) at six different prescribed speeds, from 0.1 to 5000 mm/min. The numerical simulations use a previously proposed cohesive-zone model (CZM) based on fractional viscoelasticity and a novel finite element combining a Timoshenko beam and an interface element. The CZM had previously been validated for a rubber interface, so here we present a procedure to identify its input parameters and validate its capability to predict the failure of joints made with an epoxy adhesive. An effective procedure is also developed to evaluate the dependence of the fracture energy on the crack speed without experimentally measuring the crack speed. The adhesive response was found to be markedly rate dependent. Within the range of tested speeds, the fracture energy of the adhesive more than doubles its value and the shape of the ‘fracture energy-crack speed’ curve resembles a sigmoidal shape, but more tests are needed at higher speeds to better determine the maximum value of the fracture energy and the actual shape of the complete curve.


Introduction
Adhesive joints are nowadays used in a very wide-ranging variety of structural engineering applications, because of the associated reduction of stress concentrations and their suitability for joining lightweight components. Although a number of procedures are available to characterise the interface fracture energy in such structures for quasi-static problems, the increasing use of adhesive bonding in the automotive and aerospace sectors makes it particularly important to evaluate the rate dependence of the fracture energy, because of the importance of modelling these structures under dynamic loading, and in particular during impact loading. This is more challenging, which is why a significant amount of research is still conducted in this area [1][2][3][4][5][6][7][8][9][10][11][12][13].
One problem that arises is how to account for inertial effects in the analysis of experimental result. To this end, Lißner et al. [3,4] tested specimens in mode-I butt joints, mode-II single-lap joints and mixed-mode joints made with end caps bonded on an interface at 45 deg to the load application, and at the highest speed, they used the Split-Hopkinson Tensile Bar (SHTB) rig because the results of these tests are not practically affected by inertia effects. On the other hand, these specimens do not have a predefined crack tip, so that fracture mechanics approaches could not be used. A Split-Hopkinson Pressure bar (SHPB) rig and digital image correlation (DIC) were also used by Lißner et al. [5], but in this case, the SHPB was effectively only used to apply highspeed dynamic loading on different types of specimens, namely wedge-DCB (WDCB) in mode I, end-notched-flexure (ENF) specimens in mode II and single-leg-beam (SLB) specimens for mixed mode, relying on DIC and a bespoke postprocessing technique to measure the position of the crack tip at each time. High-speed video acquisitions were also used by Blackman et al., [6,7] who conducted double-cantilever-beam (DCB) tests at several speeds and included a correction term to the data-reduction scheme when the kineticenergy contribution to the fracture energy was higher than 5% of the conventional quasi-static contribution. More recently, a similar approach was used by Sun et al. [8] to filter dynamic effects from the results of DCB tests.
One common characteristic of these valuable contributions is the use of high-resolution high-speed cameras and other state-of-the art techniques to experimentally determine the speed of crack and/or the crack tip opening. While this is feasible for research purposes, it requires expensive equipment, significant expertise and time-consuming procedures, which does not make them cost-effective for industrial standardisation. Furthermore, even for research purposes, the difficulty of accurate measurements of the crack-tip position, advance and opening displacement makes it difficult to characterise a large number of adhesives and adherents.
Another issue is how to account for rate dependence in computational modelling. Using traditional fracture mechanics approaches, the fracture energy can be taken as a function of the crack speed, regardless of whether the fracture energy is calculated as the critical values of the energy release rate, the stress intensity factor or the J integral. However, the crack speed is not known a priori, particularly for irregular structures, and therefore in the recent two decades cohesive-zone models (CZMs) have gained popularity as an alternative approach.
A number of CZMs have been proposed to account for the rate dependence of fracture in adhesive joints or delamination of composites. Among them, Corigliano et al. [14] proposed two models, one based on softening plasticity and a second one in which a rate-dependent damage evolution law is adopted. Allen and Searcy [15] developed a CZM via homogenisation of the response studied at the micromechanical scale. Xu et al. [16] a developed rate-dependent CZM using a rheological model combining a rate-independent and a rate-dependent Maxwell element. Musto and Alfano [17] proposed a model combinding damage mechanics and viscoelasticity, in which the stress provided by a viscoelastic rheological model, made of one elastic arm and one Maxwell arm, was essentially scaled via a damage variable, whose evolution was taken as rate independent. Essentially, the underlying idea was that there is an intrinsically rate-dependent fracture energy in the interface, associated with the energy of the bonds, so that the rate dependence of the response is only the result of the viscoelastic dissipation occurring during the debonding process.
All these CZMs are found to be effective in modelling rate dependence and to some extent are formulated in an attempt to reproduce the underlying mechanics, but their range of validity is limited to the range of speeds for which the input parameters are determined. More recently, Lißner et al. [3][4][5] used a CZM in which the peak stress and the fracture energy are taken as a logarithmic function of the rate of relative displacement, essentially justified by curve-fitting of the experimental test results, whereas the stiffness of the interface is assumed to be rate independent. Again, given the empirical nature of this approach, it is difficult to expect that the predictive capability of the model extends beyond the range of tested speeds.
These limitations were the motivation for the CZM developed by Musto and Alfano, first in [18] and then revisited and extended in Ref. [19] The general idea is the same as in Ref. [17] but the model is formulated in the context of fractional viscoelasticity, which allows simulating the rate dependence of polymers and elastomers within a very wide range of strain rates with a much smaller number of input parameters than those required by models based on an exponential kernel. This is because for these types of materials, the entire relaxation function is well approximated by a power law, which leads to the derivation of rate forms of the constitutive law that involve derivatives with respect to time of non-integer (i.e., fractional) order. This idea of combining fractional calculus with cohesive-zone models has then been followed by other authors [20,21] and the formulation presented in Ref. [18] with a very similar numerical implementation was combined with the widely used Park-Paulino-Roesler cohesive-zone model in Ref. [22] However, the model proposed in Ref. [18] was only validated against experimental results for a rubber interface joining two steel plates and tested in a double-cantilever beam (DCB) configuration. The validation showed that the model was able to capture very well, and with a unique set of seven parameters, the load-displacement curves at different loading speeds, ranging from 0.01 to 500 mm/min, leading to a monotonic 'fracture energy-crack speed' relationship, with a sigmoidal shape. In Ref. [19], its formulation was recast within a thermodynamic framework, which led to consider a wider set of options for the damage evolution law. It was found that, depending on which part of the free energy is assumed to drive damage, different types of rate-dependent behaviour are obtained. In particular, if only the immediately available part of the free internal energy is assumed to drive the damage, according to the model, the fracture energy as a function of crack speed turns out to be monotonically increasing, with the sigmoidal shape of the graph found in Ref. [18] If the rest of the energy is included, the same function is found to first increase and then decrease, with a bell shape of the graph. It was suggested that the latter case could be more representative of the behaviour of polymers, but no quantitative validation against experimental results was provided.
Experimental results obtained by other researchers do not provide a common qualitative trend of the relation between fracture energy and crack speed for different adhesives. For the rubber-toughened single-part epoxy adhesive Betamate XD4600, using different specimens and substrates, Blackman et al. [6] found practically no rate dependence at low speeds and a decrease in the fracture energy with testing speed at very high speeds, in which dynamic effects were significant. However, they could attribute the decrease in fracture energy to the change in material properties caused by heating, which in turn is due to the transition from an isothermal to an adiabatic process. For the thermosetting epoxy adhesive AF 163-2OST and adherents made of titanium alloy Ti-6Al-4 V, Lißner et al. [3] found a decrease in both the strength and the fracture energy with increasing loading speed. Sun et al. [8] found an increase in the fracture energy with loading speed for the toughened single-part epoxy adhesive SikaPower-497, whereas for the twopart polyurethane adhesive Araldite-2028 a non-monotonic trend was found with an initial increase followed by a decrease in fracture energy with increasing loading speed. This last result could resemble the bell shape postulated by Alfano and Musto in Ref. [19] but Sun et al. [8] could attribute the final decrease to the transition from a combined interfacial-cohesive failure at lower speeds to a purely interfacial failure at the highest speeds.
In this paper, we present a numerical and experimental study of the rate dependence of the failure of adhesive joints based on DCB tests. For the experimental part, we tested 24 adhesive joints made of aluminium Al 2082-T6 bonded with the epoxy adhesive Araldite® 2015, in a DCB configuration with a prescribed cross-head displacement rate ranging between 0.1 and 5000 mm/min, with 4 specimens tested at each speed.
To avoid the challenges and the lack of accuracy related to real-time measurements of the crack tip position, we develop a post-processing procedure that provides an excellent approximation of the relation between crack speed and fracture energy, without the need for measuring the crack length and the speed of its advance during the test. The method is an application of the data-reduction scheme proposed in Ref. [23] which is based on the enhanced simple beam theory (ESBT), and uses the concept of equivalent crack length, introduced by De Moura et al. [24] in their compliance-based beam method (CBBM). In this respect, the proposed method is similar to the approach used by Dagorn et al. [25] which proved its accuracy by comparing the indirectly computed crack speed with that measured by means of a high-resolution video camera. Unlike the approach proposed in Ref. [25], our implementation does not require the assumption of a constant fracture energy and does not neglect shear deformability of the arms. By using the ESBT data-reduction scheme, which accounts for the crack-root rotation in combining Timoshenko beam theory with linearelastic fracture mechanics (LEFM), we increase the accuracy of the procedure and make it particularly suitable for composites. It is worth noting that, with our approach, we characterise the fracture energy using the critical energy release rate G c , justified by the rigorous analysis presented in Ref. [26]. Furthermore, it was shown in Ref. [23] that the ESBT data-reduction scheme is more reliable and accurate than those currently used in standards [27,28], see also [26].
We also present a systematic parameter-identification procedure to determine the input constants for the fractional rate-dependent model proposed by Musto and Alfano [18,19] and experimentally validate the model for the cases tested. Furthermore, the numerical simulations reported in this paper have been conducted by developing a special finite element that combines together a Timoshenko beam finite element with a linear interface element, within one single element. This is a special reduced case of the multi-layer beam model presented by Škec et al. [29] where it has been shown that a combination of beam finite elements and interface elements can be an accurate, fast and robust alternative to the common approach with 2D-solid finite elements. Such properties are extremely important in the context of rate-dependent analyses because an accurate identification of interface parameters requires several runs of multiple simulations (each one with a different load-line displacement speed).
The outline of the paper is as follows. The numerical methods used are presented in Section 2. In particular, the new finite element used in the numerical simulations is first described in Section 2.1. Then, in Section 2.2, for the rate-dependent CZM formulated in Ref. [18], we briefly provide its governing equations and the other essential features that we later refer to, for the sake of self-consistency. The experimental tests conducted and their results are described and discussed in Section 3. In Section 4, we describe the procedure developed to determine the rate dependence of the fracture energy in terms of crack speed and discuss the results obtained in our case study. In Section 5, we describe the identification procedure developed to calibrate the input parameters of the rate-dependent CZM and compare the obtained numerical results with the average load-displacement curves. Some conclusive remarks are finally provided in Section 6.

Combined beam-interface finite element
In the numerical model, a novel finite element, specially tailored for the DCB test, is used. In this element, the standard geometrically linear two-node Timoshenko beam finite element, with one integration point to avoid shear locking, is combined with a 4-node interface element (INT4). Because in a DCB test there are no axial forces in the arms if a geometrical linear model is used, this type of element has only two degrees of freedom per node, for a total of 4-element degrees of freedom, namely the two nodal vertical displacements v 1 and v 2 and the two nodal cross-sectional rotations θ 1 and θ 2 (see Figure 1).
Due to symmetry, only the upper half of the specimen is modelled. The contribution of the interface is simply added to the corresponding elements of the beam's element vector of residual forces and the stiffness matrix at the corresponding positions. Therefore, the main advantage of using this element is that the only degrees of freedom of the DCB model are those of the nodes of the beam elements. Because the two-node Timoshenko beam elements are well known and reported in many FEA textbooks [30], here we focus on the contribution of the interface, as explained in the following subsection.

Contribution of the interface in the mode-I beam finite element
It is assumed that the interface mode-I stresses at coordinate x, denoted by σðxÞ, act on the beam finite element as a distributed load qðxÞ, as shown in Figure 1. It is further assumed that the distribution of interface stresses is constant along the width of the specimen and that (mode-I) tensile stresses are positive. Because the y-axis of the beam's coordinate system, unlike the positive interface stresses acting on the beam, is pointing upwards (as shown in Figure 1), it follows that qðxÞ ¼ À b σðxÞ. Finally, the virtual work W � of q done on the beam can be written as where b is the width of the specimen (interface), L e is the length of the beam finite element, v is the virtual transversal displacement of the beam and σ is the distribution of mode-I interface stresses along the interface, the last two being function of the position x on the interface. The interface stresses are unknown and depend on the formulation of the CZM, which will be discussed in the next subsection. For a two-node beam finite element, linear-interpolation functions are used for the transversal-displacement and cross-sectionalrotation fields so that where the unknown virtual kinematic fields are contained in vector pðxÞ ¼ vðxÞ θðxÞ � � T , the nodal values of virtual transversal displacements and cross-sectional displacements are contained in vector is the matrix containing linear interpolation functions ψ 1 ðxÞ ¼ 1 À x=L e and ψ 2 ðxÞ ¼ x=L e . Because the vector of residual forces is computed as the difference between the vector of internal and the vector of external forces, which can be written as g ¼ q int À q ext , the contribution of the interface to the vector of residual forces is actually positive and, from Equation (1), it can be written as where βðxÞ ¼ ψ 1 ðxÞ 0 ψ 2 ðxÞ 0 � � . The contribution of the interface to the stiffness matrix is obtained by linearising the vector of residual forces (4) with respect to the vector of the nodal element degrees of freedom. As it will be shown in the next subsection, the interface stresses σðxÞ depend on the relative displacement δðxÞ, where due to the symmetry of the DCB test δðxÞ ¼ 2 vðxÞ. Thus, in the residual (4), it is necessary to linearise only σ with respect to v. Finally, the contribution of the interface to the element's tangent stiffness matrix can be written as For the numerical integration of the finite-element, the contribution of the beam element is integrated using the 1-point Gauss quadrature, while for the contribution of the interface, the 3-point Simpson's rule is used. Note that in this subsection, for the sake of simplicity, the explicit dependence of variables on time has been omitted. In the next subsection, however, such dependence will be made explicit because of the rate dependence of the CZM.

Rate-dependent CZM model
The interface is modelled using the mode-I rate-dependent CZM proposed by Musto and Alfano [18] which combines damage mechanics with fractional viscoelasticity and whose governing equations are recalled and discussed here because the physical meaning of the parameters is important for the analysis presented later in this article.
The rheological representation of the model, in the absence of damage, is shown in Figure 2, where σ indicates the interface traction that would be obtained as the response to the prescribed history δ if no damage occurs. Although similar to the standard linear solid (SLS) (or Zener model), the model proposed in Ref. [18] has several additional features that allow for the modelling of progressive rate-dependent damage at the interface. First of all, the dashpot of the SLS model is replaced by the 'Scott Blair' (SB) element, in which the stress is obtained as where α is the relative displacement within the Scott Blair element, η is a material constant and the second term is the fractional derivative of order ν of α with respect to time t, so that ν is a second material constant, with 0 < ν < 1. Note that in the original paper with the fractional rate-dependent CZM [18] the authors used the term 'spring-pot' for what they have referred to as the 'Scott Blair' element since their following publication [19] and we are following the latter terminology here.
We will refer to the bottom and top arms of the model in Figure 2 as elastic and inelastic, respectively. The stiffness of the spring in the elastic arm is denoted by k 1 , while the spring in the inelastic arm has stiffness k 2 . The normal interface traction at time t, σðtÞ, is given by the following non-linear functional: where D is the damage variable ranging between 0 (no damage) and 1 (full damage), and δ in Equation (7) refers to the entire relative-displacement history δ : 0; t ½ � ! <, which is why σðtÞ is a functional. Both springs in the model are assumed to be linear elastic, so that the following governing differential equation is obtained: where λ ¼η=k 2 and γ ¼λðk 1 þ k 2 Þ.
For reasons that will be clarified later, here we use the damage-evolution law used in Ref. [18] instead of one of the alternative options proposed in Ref. [19]. In the latter article, it was shown that adopting the damage evolution law used in Ref. [18] is equivalent to assume that the damage is driven only by the energy stored in the elastic arm. In other words, in this model D grows when the energy in the elastic arm, k 1 δ 2 =2, reaches a critical threshold, which itself is a function of D.
Assuming a bi-linear traction-separation law (TSL) in the slow limit, as defined in Figure 3, the damage variable is computed as follows: where δ 0 and δ c are the values of the relative displacement corresponding to damage (softening) onset and total debonding at the interface, respectively, while δ max is the maximum relative displacement in the previous history.
In the slow limit, the inelastic arm is completely relaxed so that α ! δ and the stiffness of the model becomes k 0 ¼ k 1 . In the fast limit, because α ! 0, the stiffness of the model becomes k 1 where Ω i and σ i , with i ¼ 0 for the slow limit and i ¼ 1 for the fast limit, are the work of separation and the maximum interface traction, respectively, as shown in Figure 3.
It should be noted that, because the relative displacements δ 0 and δ c are rate-independent model parameters (see Figure 3), the evolution law for the damage variable D, as defined in Equation (9), is rate independent too. The overall damage process turns out to be rate dependent though because of the interplay between damage and visco-elastic response of the interface.
Between the slow and the fast limit, a sigmoidal variation of the work of separation with respect to the relative-displacement speed is obtained [18]. In Ref. [19], the authors considered different damage-evolution laws in which, besides the energy in the elastic arm, damage can be driven by the energy in both springs or the entire free energy (including the SB element). Such damage-evolution laws give bell-shaped variations of the work of separation with respect to the relative-displacement speed. However, such formulations have not been used in the present work because there has not been experimental evidence of a non-monotonic change of the fracture energy with respect to the relative-displacement speed, although this could be due to the fact that not sufficiently high speeds were reached in the tests conducted.
It needs to be emphasised that in compression (for δ < 0), no damage is allowed at the interface. Moreover, it is assumed that the rate dependence of the interface in compression is negligible so that the interface stiffness in compression is given by a scalar k c . For the time implementation of the rate-dependent cohesive model, which is based on the Grünwald-Letnikov definition of the fractional derivative, the reader is referred to [18].
It is worth noting that in this model, the CZM input parameters are not rate dependent. Instead, this CZM has fixed parameters that define the fast and the slow limit, while the transition between the two limits, which essentially defines the rate dependence of the model, is governed by two (also rate independent) parameters λ and ν (see Equation (8)).

Parameters of the rate-dependent CZM and their identification
The model has seven independent parameters in total, and it is therefore useful to discuss the role of each of them in the simulation of the rate-dependent behaviour of the adhesive because this will be the basis for the parameteridentification procedure that will be presented in Section 5.
As mentioned earlier, k c is given as a separate parameter in order to ensure that the behaviour of the interface in compression is rate independent. In Ref. [18] it has been shown that, within certain limits, this parameter does not have any noticeable influence on the results of a DCB test simulation. Thus, we take the stiffness of the interface in compression as the mean of the stiffness values in tension at the slow and the fast limit, which gives The next three parameters, namely Ω 0 , σ 0 and k 0 (or k 1 ), are those that fully define the bi-linear traction-separation law in the slow limit. If experimental data at slow speeds suggest that the slow limit might have been reached (or at least closely approached), as it is the case in the experimental results presented in Section 3, the values of the slow-limit parameters can be identified using the force-displacement data of the slowest test. Although for identifying Ω 0 there are several accurate approaches based on the critical energy release rate [23,24] or the J integral [26,31], identifying σ 0 and k 0 is less common and more demanding. In this work, the identification of all slow-limit CZM parameters will be done using a fast and accurate method implemented in the software application DCB PAR [32] (more details will be given in Section 5). The fast limit is defined by parameter �, and its value can be easily identified if the experimental data suggest that the fast limit might have been reached or closely approached. In that case, using the load-displacement data of the fastest test, a rate-independent parameter identification (identical as for the slow limit) can be performed. It is sufficient to identify the fracture resistance at the fast limit, Ω 1 , from which it follows that � ¼ Ω 0 =Ω 1 .
The two remaining parameters, namely λ and ν, define the rate-dependent behaviour between the slow and the fast limit (see Equations (6) and (8)).

Experimental tests and results
Experimental tests have been carried out on aluminium DCB specimens bonded with the structural adhesive Araldite® 2015 over a wide range of test speeds.

Experimental set up
Geometrical properties of the specimens are reported in Figure 4 and Table 1. Because the aluminium plates were not thick enough for the loading arrangement with loading holes drilled through the arms (as suggested in Ref. [27,28,33]), loading blocks screwed to the aluminium plates by means of conical-head screws were used, as shown in Figure 5(a). Thus, countersunk holes on the aluminium plates and threaded holes in the loading blocks were prepared for the conical-head screws. An additional cylindrical hole, used to attach the specimen to the tensile testing machine, was drilled through each loading block, whose height c was dictated by the need of accommodating the threaded holes (see Figure 4). It is worth noting that we did not take into  Table 1. account the presence of the loading blocks in Sections 4 and 5, but we provide a derivation in the Appendix showing that the loading blocks alter the value of the computed fracture energy of about 3%. On the other hand, incorporating the analysis from the Appendix into the procedures developed in Sections 4 and 5 would not allow us to use the closed-form expression employed in Section 4 and in the parameter-identification procedure implemented in the software DCBPAR in Section 5. Given that the loss of accuracy is quite less than the scatter of the results, we prefer to neglect the influence of the loading blocks.
After attaching the loading blocks to the aluminium plates, the aluminium surfaces were roughened using sand paper and cleaned using acetone. As shown in Figure 5(a), some aluminium foil was applied on both plates to create the initial notch at the desired position, at a distance a 0 from the loading line.
The adhesive was applied using a thin cylindrical wooden stick that was first dipped into the adhesive and then rolled over the aluminium surface. By using a spatula to additionally spread and level up the adhesive, a relatively thin and uniform adhesive layer was created on each plate, as it can be seen in Figure 5 (a). The two arms were then carefully put in place and secured with clamps at both ends, as shown in Figure 5(b). In order to create a uniform adhesive layer at the interface and squeeze out the excess glue, the specimens were first pressed together with a force of 60 � 10 N for 60 s, and then loaded by two 1 kg weights that were put on a rectangular wooden bar used to uniformly transfer the load to the plates, as shown in Figure 5(b). After curing at room temperature for 24 h, the specimens were oven-cured at 80°C for 60 minutes. Prior to testing, the specimens were allowed to cool down to room temperature, and then, their thickness was measured with a digital caliper having a resolution of 0.01 mm. The difference with the initial measurement without the adhesive provided the thickness, the average of which was 0.25 mm, with a standard deviation of 0.08 mm. This relatively large deviation from the average is due to the fact that the thickness of the interface was not controlled during the process of preparation of the specimens. Although a dependence of the fracture energy on the thickness has been reported in the literature [4,34], for none of the tested any evidence of non-negligible dependence when comparing the results of the 4 specimens, which is also consistent speeds we found with the experimental results in Ref. [35] for the same adhesive and for values of thickness less than 3 mm. Therefore, it was concluded that the scatter of the experimental data can be primarily attributed to the imperfections at the interface that were similarly distributed in interfaces of different thicknesses.
Debonding test have been performed at six different load-line-displacement speeds, namely 0.1, 1, 10, 100, 1000 and 5000 mm/min. For each speed, four tests have been made. All tests were performed at room temperature using a 30 kN electromechanical Instron testing machine, except for the tests at 5000 mm/min, which were performed using a 100 kN servo-hydraulic Instron machine. The experimental set-up is shown in Figure 6.

Experimental results: load-displacement curves
The raw data (containing the time, the load-line displacement and the applied force) have been post-processed in order to allow a consistent comparison between test results and with the results of numerical simulation, in which selfweight was not included.
First, the weight of the specimen was subtracted from the measured load while at the same time ensuring that the remaining part of the initial preloading was included in the measurement. To this end, the load was initialised to zero before the start of the test, when the specimen was mounted, that is when the load applied included both the effects of the self weight and that of the small additional preloading needed to make arms approximately horizontal. The (negative) load measured at the end of the test, after the DCB was completely broken, minus the weight of the half specimen still attached, gave us the total value of the initial preloading. This value, minus the total self weight of the specimen (including the loading blocks) was then added to the measured values of the load during the test.
Next, the zone with minimal compliance (maximum slope), which represents the pure linear-elastic behaviour, was identified in the early stage of every experiment. Due to the inevitable small differences in the initial preloading, the linear fit of these measurements gives a different non-zero intercept with the force-axis. Therefore, all load-line displacements were corrected by subtracting this constant value so that the linear fit of the corrected data passes through the origin of the force-displacement co-ordinate system.
To obtain the final corrected load-displacement curve from each test, the (linear elastic) compliance of the rig, C rig , was measured, and the associated displacements, C rig F, were subtracted from the measured displacements, where F is the measured force.
Moreover, in order to compute an 'average' load-displacement curve, we computed an average force for a sufficient number of displacement points. To this end, the data was modified by defining a constant 0.2 mm displacement increments and interpolating the force values that correspond to the closest displacement below and above the desired one.
The post-processed force-displacement data for all 24 tests with displacement increment equal to 0.2 mm is shown in Figure 7, where the results for each speed are given in a separate plot. Individual tests are labelled as s À i, where s ¼ f0:1; 1; 10; 100; 1000; 5000g represents the speed of the test in mm/ min, while i ¼ f1; 2; 3; 4g is the test number. In addition, the average loaddisplacement curve for each speed is given by a solid red line. The average value of the force is computed only for values of the displacements for which in all of the four selected tests, a force has been recorded before the total debonding of the plates. Therefore, the maximum value of the load-line displacement of the average curve for each speed is the minimum of the four values of the displacement corresponding to failure for each specimen, not their average. This is the reason why in Figure 8 the increase in the maximum load-line displacement with the loading rate, unlike that of the corresponding applied load, is not perfectly monotonic. Such an increase is due to the fact that the fracture energy also increases with the loading rate, as shown in Figure 9(b). From the expressions derived in Section 3.3, it clearly follows that for the same crack-tip position, increasing the fracture energy will also increase the applied load and load-line displacement. This is important because the last measured point for each test is when the crack-tip position is close enough to the edge of the specimen to result in a sudden drop of the force. Therefore, the last measured value corresponds to a crack-tip position which is practically the same for all tests, and therefore, for higher test speeds in general, we normally obtain larger final values of load and load-line displacement.
Although most of the failure was cohesive, in all specimens a small amount of interfacial failure was also found, typically on 5% or less of the bond area. Another small amount of the fracture surface, typically between 5% and 10% or less of the bond area and in small patches distributed in an apparently random way, was also characterised by a smoother and shinier appearance. These may be sites of imperfect bonding and, in conjunction with the small patches of interfacial failure, can partly explain the experimental scatter in the results, not only between different specimens but also within the same specimen, the latter explaining the oscillations in the individual load-displacement curves. Similar scatter and oscillations were found by other researchers for similar epoxy adhesives [36][37][38] and for the same adhesive used in our work [35,39]. The scatter is consistent across all speeds and has been quantified by computing the average value of the coefficient of variation for each speed given in Table 2.

Test results in terms of average load-displacement curves
The average load-displacement curves for each speed (represented by solid red lines in Figure 7) are now compared in a single plot in Figure 8. It is clearly visible that the adhesive has shown a considerable amount of rate dependence because the fracture energy of the adhesive (and therefore the overall bearing capacity of the adhesive joint) obviously increases with the load-line displacement speed. However, it can be noticed that the two lowest speeds (namely 0.1 and 1 mm/min) essentially overlap, which suggests that the slow limit (i.e., the limit below which the rate dependence of the adhesive is negligible [18]) has been reached. Results for the two highest speeds (namely 1000 and 5000 mm/ min) are also very close for displacements greater than 3 mm, but the peak load for 5000 mm/min is significantly higher than for 1000 mm/min. Therefore, despite the results providing a good assessment of the rate dependence over a wide range of speeds, in order to characterise the behaviour when the fast limit is approached more tests would be needed, which could not be conducted in this project because the maximum speed of the testing machine available was 5000 mm/min.

Determination of the fracture energy as a function of crack speed
The rate dependence of fracture resistance is often reported in the literature in terms of fracture energy against crack speed. Producing such type of data accurately is very challenging, time-consuming and expensive. This is because, to measure the crack speed, the crack length needs to be measured in the first place. If this is done optically, the challenge is that the crack tip is not always clearly visible because what appears to be the crack tip is often at a certain distance from the crack, and such distance can often vary for different values of the prescribed load or displacement. For tests at high speeds, an additional challenge consists of acquiring images of the area around the crack tip with sufficient resolution, which requires the use of expensive cameras and significant time for post-processing results. Alternative experimental methods can also be used, but they are time-consuming and still require the use of expensive equipment. Last but not least, the accuracy achieved with these methods significantly depends on the expertise and experience of the person conducting the test.
For the specific case of a DCB, the crack speed significantly varies during a test even if a constant value of the line-load-displacement speed is prescribed, as is shown here as well as in Ref. [25] It is interesting though that, although evaluating the crack speed at each time of the test is an additional challenge, when such measurements are available (directly or indirectly, the latter approach being used here), rate dependence within a non-negligible range of crack speeds can be assessed within a single DCB test performed at a constant value of the line-load-displacement speed.
To overcome the difficulties related to the measurement of the crack speed, in this paper, we determine the variation of the fracture energy as a function of the crack speed, based on the recently developed enhanced simple beam theory (ESBT) data-reduction scheme [23] used to estimate the critical energy release rate G c . The method does not require any direct measurement of the crack length, which makes it very time-and cost-effective and therefore particularly convenient for industrial applications. With our approach, based on the Timoshenko beam theory and the concept of equivalent crack length, the fracture energy, i.e., the critical energy release rate, is computed as (11) where F is the applied load, b is the width of the specimen, EI and μA s are the bending and shear stiffness of each arm, respectively, while α ¼ μA s =EI. The equivalent crack length is defined as a eq ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi is the equivalent crack length that is obtained if the Euler-Bernoulli beam theory is used. Note that in the above equations, δ is the load-line displacement, EI ¼ E � I and μA s ¼ μ � A � k s , where E is Young's modulus, μ is the shear modulus, k s ¼ 5=6 is the shear correction coefficient [40], while the area and the second moment of area for a rectangular cross-section of the arms read A ¼ hb and I ¼ bh 3 =12. Because the actual crack length is not provided in the ESBT data-reduction scheme, the values of the fracture energy are plotted against the equivalent crack length, as shown in Figure 9 (a). Therefore, instead of the resistance curve (R-curve), with the ESBT datareduction scheme the equivalent R-curve can be obtained. However, it is also shown in [26] and [41] that the derivative da eq =da of the equivalent crack length a eq with respect to the actual crack length a is so close to unity that, in practical terms, the equivalent R-curve is simply shifted to the right by a constant value of Δa with respect to the actual R-curve, which has no influence on the crack speed. Therefore, the equivalent R-curve can still be used to accurately represent the values of the fracture energy, as a function of both the crack length (except for a constant shift) and of the crack speed.

Fracture energy as a function of the equivalent crack length (R-curve)
The equivalent R-curves in Figure 9(a) were obtained by taking the average from the four tests at each speed. In turn, for each speed, Figure 9(b) reports the average of the R-curves of Figure 9(a). Both figures confirm the ratedependent behaviour of the adhesive. In Figure 9(b), it can be noticed that the fracture energy between the slow and the fast limit increases more than twice. Interestingly, the minimum mean value of the fracture energy is obtained for 1 mm/min, but the difference with the mean value obtained for 0.1 mm/min is very small and below the experimental scatter reported by the error bars in Figure 9(b). It is worth noting that Figure 9(b) does not provide an accurate representation of the rate dependence of the fracture energy of the adhesive because the latter undergoes different strain rates at different points of the interface, even for a constant rate of the prescribed load-line displacement. One way of addressing this issue was proposed by Nunes et al. in Ref. [42] who prescribed a variable load-line displacement rate determined in such a way to result in a constant strain rate at a single point of the interface. An alternative approach used here is to determine the relationship between the fracture energy and the crack speed, which is done in the next section.

Fracture energy as a function of the equivalent crack speed
Since the equivalent crack length is obtained for known values of time, the equivalent crack speed can now be computed. To this end, the equivalent crack length vs time is interpolated by a power function, whose derivative gives the equivalent crack speed. Such approach, in comparison to the finite-difference approach, gives smoother results with less scatter. Therefore, with our procedure, we can obtain a number of points for the fracture energy as a function of the equivalent crack speed, without any need to measure the actual crack length or crack speed. As discussed earlier, the difference between this function, reported in Figure 10, and the fracture energy as a function of the actual crack speed is practically negligible.
Not only does Figure 10 clearly show the rate dependence of the adhesive, it also gives additional information about the response for the lowest and highest speeds that is not available in the force-displacement plot (Figure 8) or the R-curve (Figure 9(a)). First of all, it can be noticed that the equivalent crack speed varies within a single test (here represented by an average of four tests). In fact, it can be easily shown that for a constant load-line displacement speed, the crack speed decreases during crack propagation, at least in the simple assumption that the arms of the DCB are clamped at the crack tip and using Euler-Bernoulli theory, for which the load-line displacement (or the crackmouth opening) can be written as Replacing the expression of F in this simplified case (see also Equation (11) for a ¼ a eq and neglecting the shear deformation): we get: Taking the time derivatives of both sides and solving for _ a we have:  ffi ffi ffi ffi ffi ffi ffiffi EI G c b r _ δ (17) which shows that the crack speed is inversely proportional to the crack length for constant _ δ. In Figure 10, it can also be seen that for each load-line displacement speed, except for the test speed of 1000 mm/min, the highest value of the fracture energy is obtained for values of the crack speed equal to or close to the highest one (that corresponds to the beginning of the crack propagation). For every load-line displacement speed, except for the tests with loading speed of 0.1 and 1000 mm/min, an increase in the fracture energy between the minimum and the maximum equivalent-crack speed can be noticed. This can also be noticed for 5000 mm/min, which would suggest that we cannot claim that the fast limit has been reached. This is regardless of the actual shape of the complete curve relating the fracture energy to the crack speed. If this curve had a bell shape, that is one of the possibilities suggested in Ref. [19], then the fast limit would certainly be very far away. If instead the curve had a sigmoidal shape, then the results for the highest speed could potentially be not far from those in the fast limit, but there does not seem to be sufficient evidence to claim that the fast limit has been reached. This is also in consideration of the significant scatter in Figure 10 for the highest loading speed of 5000 mm/min.
The results for the test speed at 1000 mm/min can be explained by noticing that the individual load-displacement curves in Figure 7 for 1000 mm/min are considerably smoother than all the remaining results at other speeds. This is not the case for the other speeds and can be attributed to the fact that 1000 mm/min was the highest speed tested in the 30 kN testing machine, and for this speed, the measurement interval used must have resulted in automatic time-averaging of measured data during its acquisition in the tensile-testing machine software.
On the other hand, results for 0.1 mm/min indicate that the slow limit has probably been reached.

Parameter identification using the rate-dependent numerical model
In this section, the rate-dependent numerical model presented in Section 2 is used to simulate the DCB experiments presented in Section 3. The goal is to identify the values of the parameters of the rate-dependent CZM so that the average experimental force-displacement data is accurately captured over six different load-line displacement speeds and with a single set of input parameters.
To this end, an assumption is made regarding the behaviour at the highest tested speed. We noted that, based on Figure 10, there is not enough evidence to support the fact that the complete shape of the 'fracture energy-crack speed' curve is sigmoidal. However, the parameter-identification method used in this section is based on the mean fracture energy, computed for each test speed and reported in Figure 9(b), which does result in a sigmoidal shape of the relationship. We therefore assume, for the sole purpose of identifying the input parameters, that the curve is indeed sigmoidal and, furthermore, we also assume not only that the slow limit is reached for the lowest tested speed of 0.1 mm/min but also that the fast limit is reached for the highest speed of 5000 mm/min. With this assumption, if we obtain good correlation between numerical and experimental results within the range of tested speeds, we can expect that the model keeps its predictive capability for lower speeds than those tested, because the slow limit seems to have been reached and due to the typically observed predictive features of fractional viscoelastic models outside the ranges of tested speeds. However, we have no sufficient evidence to expect that the same happens for speeds higher than those tested.
The identification procedure comprises three steps, reported in the following three subsections.

Identification of the slow-limit CZM parameters
As discussed in the previous section, we assume that the slow limit has been reached. Therefore, the slow-limit CZM parameters (Ω 0 , σ 0 and k 0 ) are identified using the procedure proposed in Ref. [32] and implemented in the DCB PAR software, which is based on an analytical solution with a rateindependent bi-linear TSL at the interface [41]. After loading the experimental load-displacement data (the average load-displacement curve for 0.1 mm/ min), and entering the geometrical and material properties of the plates (as reported in Section 3), the software returns Ω 0 ¼ 0:2605 N/mm, σ 0 ¼ 16:9528 N/mm 2 and k 0 ¼ 4363:4 N/mm 3 . A comparison of the experimental loaddisplacement data and the analytical solution from [41] is given in Figure 11 (a). Excellent agreement before and during crack propagation can be noticed.

Identification of the fast-limit CZM parameter � and of the compressive stiffness k c
Likewise, as earlier explained, it will be assumed that the average load-displacement curve for 5000 mm/min corresponds to the fast limit. This curve is therefore used as an input for the DCB PAR software to obtain a new set of parameters Ω 1 ¼ 0:5682 N/mm, σ 1 ¼ 35:8366 N/mm 2 and k 1 ¼ 7642:1 N/mm 3 .
Because � ¼ Ω 0 =Ω 1 ¼ σ 0 =σ 1 ¼ k 0 =k 1 , using the DCB PAR results for the slow and the fast limit, � can be computed in three different ways, that is as The first value � ¼ 0:458 was used because it gives the best agreement with the experimental data during crack propagation.
All values of the identified (Ω i , σ i , k i ) and computed (δ 0 , δ c ) parameters are given in Table 3, where the chosen value of � is denoted in boldface font. It should be noted that, based on the identified parameters for the slow and the fast limit, δ c is indeed found to be practically rate independent, which is in accordance with the assumption made in Section 2.2. Although δ 0 should be also rate independent, the difference between computed values for the slow and the fast limit is bigger (but still relatively small) due to the difference between slow/fast ratios for Ω i and k i (as given in Table 3 for i ¼ 0; 1).
A comparison of the experimental load-displacement data at 5000 mm/min and the analytical rate-independent solution from [41] is given in Figure 11(f).

Identifying the rate-dependent CZM parameters λ and ν
Having determined the slow-and fast-limit parameters, identifying the two remaining rate-dependent CZM parameters is relatively simple and was done by using a trial-and-error method, first for ν and then for λ , facilitated by the following observations as well as engineering judgment. Increasing ν essentially narrows the band of speeds between the slow and the fast limit. Thus, in the present case, ν has been increased until the model's solution at 5000 mm/min reached the fast limit, which corresponds to ν ¼ 0:85. Only few iterations were then necessary to find the value of λ that can accurately capture the load-displacement curves of the remaining speeds between the slow and the fast limit with a satisfactory level of accuracy. Lower values of λ essentially reduce the force in the force-displacement diagram and finally λ ¼ 0:25 has been determined. The final identified values of all CZM parameters are given in Table 4, while the comparison of the average Table 3. Values of the identified (Ω i , σ i , k i ) and computed (δ 0 , δ c ) CZM parameters for the slow (i ¼ 0) and the fast (i ¼ 1) limit. CZM  experimental load-displacement curve with the predictions obtained by the rate-dependent model for each speed is given in Figure 11. The agreement across all speeds is satisfactory. For all numerical simulations presented using the finite-element model, convergence analyses have been made in terms of both spatial and time discretisation, and the results presented are characterised by a negligible error, which is unnoticeable on the graphs.

Conclusions
We presented an experimental-numerical study to characterise the rate dependence of fracture energy of adhesive joints under mode-I crack propagation, to validate the rate-dependent cohesive-zone model based on fractional viscoelasticity previously formulated in Ref. [18] then, extended in Ref. [19], and to calibrate its input parameters.
In the experimental part, 24 double-cantilever beam (DCB) specimens, made of aluminium Al 6082-T6 plates bonded with the adhesive Araldite® 2015, were tested in displacement control, with constant cross-head displacement speed for each test, and six speeds ranging from 0.1 to 5000 mm/min. For each speed, four specimens were tested, and average load-displacement curves were obtained.
The fracture energy was found to markedly increase with crack speed, with values obtained at the highest tested speed more than doubled with respect to those at the lowest two speeds.
We also presented an effective procedure to determine the 'fracture energycrack growth' curve without the need for measuring the crack length and the crack speed, but only by closed-form post-processing of the measurements of the load and the displacement, immediately available at the end of the tests.
The numerical simulations were conducted using a novel finite-element that combines together a linear interface element with a linear Timoshenko beam element, making use of the symmetry of the configuration, as well as the cohesive-zone model proposed in Ref. [18] and previously validated only for a rubber interface.
To determine the CZM input values, we presented a systematic parameter identification procedure, which assumes that the slow and fast limits of the rate-dependent material behaviour have been sufficiently approached at the slowest and highest tested speeds, respectively. At each of these speeds, the parameters of a rate-independent CZM were computed using the procedure developed in Ref. [32] and the analytical rate-independent solution provided in Ref. [41]. Based on these, the remaining parameters of the rate-independent model were calibrated. The comparison between experimental and numerical results confirms the capability of the CZM of capturing the experimental response, also for a glassy polymer, over a wide range of speeds with the same set of seven parameters, only two of which are related to the model rate dependence. This is unlike models based on experimental kernels, which require a much larger number of parameters.
The proposed numerical model can be used to simulate the response of a DCB joint at any loading rate (constant or variable). The procedure presented in this work is convenient for obtaining the seven parameters of the model for which, obviously, some experimental data is required. However, once the parameters are known, the model can be used for different loading rates (constant or variable) without the need for any additional experiments.
Extrapolation outside the range of tested speeds (and strain rates) is possible, but in this case, as explained earlier, the predictive capabilities of the model might be limited.
The CZM used, unlike the more general formulation presented in Ref. [19], assumes that the 'fracture energy-crack speed' curve is monotonic and its graph has a sigmoidal shape. Although this assumption was useful for the purpose of the calibration of the model parameters, from the analysis of our results, we could neither rule out nor confirm that the 'fracture energy-crack speed' curve is monotonic and with a sigmoidal shape of its graph. Although the curve found is sigmoidal, we do not think we have sufficient evidence to exclude that at higher speeds than those tested in this work the monotonicity could be lost and the fracture energy could decrease. This would in turn result in some sort of 'bell shape' of the curve, such as the results theoretically obtained in Ref. [19] when damage is assumed to be driven by the entire free energy. This final observation suggests that more tests should be done to fully characterise the behaviour of the adhesive, which will be the objective of future work.
It is also worth noting that the models considered here are all visco-elastic, meaning that visco-plastic behaviour has not been considered. Adding viscoplasticity to the modelling framework is likely to affect how the fracture energy changes with crack speed, both qualitatively and quantitatively, and influence the results for non-monotonic loading. This line of research will also be pursued in the future work.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Funding
This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 701032. Appendix A. Derivation of the formula for G c that takes into account the loading block of a DCB specimen

ORCID
Let us consider the case in which the load on a DCB specimen is not applied directly to the arms, but via a loading block whose deformation during a DCB test is assumed to be negligible, as depicted in Figure A12. The force is transmitted from the tensile testing machine to the loading block by means of a pin whose centroid distance from the arm's midplane is denoted by l b (line AB). It will be assumed that the DCB arm is modelled as an Euler-Bernoulli beam, whose left-hand end corresponds to point A (see Figure A12). The stiffening effect of the block on a small portion of the arm around point A is assumed to be irrelevant because the bending of the arm is negligible close to the free end. In a deformed configuration, due to bending of the arm, the block rotates, while the applied load F transmitted from the tensile testing machine in point B remains always vertical. The force is then transferred from point B to A, which implies that the arm with a loading block is loaded not only by a transversal force F, but also by a bending couple M ¼ F e, where e ¼ l b θ for small cross-sectional rotations θ of the arm at point A (see Figure A12). Assuming that the arms are clamped at the crack tip, the following expression for the loadline displacement and cross-sectional rotation of the arm at point A can be derived, respectively: where a is the crack length and EI is the bending stiffness of the DCB arm. The second (negative) term in the parentheses of each expression represents the contribution of the loading block. It can be noticed that the presence of the loading block decreases the deformation of the DCB arm for the same values of the load and crack length. The critical energy release rate can be computed from where for a DCB with loading blocks the total potential energy can be defined as Note that the signs in the above equations are the result of assuming that on the top arm, the force is positive if upward, θ is positive if clockwise and δ is the total cross-head displacement, assumed positive when it is opening. The derivative in (A.3) can be computed assuming that either force or displacement is prescribed, which in both cases leads to the same result for G c . Although DCB experiments are usually displacement controlled, the analytical derivation of G c (especially in this case) is simpler if we assume that the force is prescribed. By doing so, it follows that In the absence of the loading block (e ¼ 0), this expression becomes the well-known expression for DCB, G c ¼ F 2 a 2 =ðb EIÞ. From (A.2), one can write: Solving for e provides: As discussed in detail in Ref. [23,26] by using the equivalent crack-length approach, G c can be computed using only the measured values of load F and load-line displacement δ, while the measurement of the crack length a (much less practical and accurate) can be avoided. The main idea is to express a from (A.1) as a function of F and δ. This is not the actual value of the crack length, but the equivalent crack length (a eq ) that, for measured values of F and δ, accommodates the assumption that the arms are clamped at the crack tip made in (A.1). In the absence of loading blocks (e ¼ 0), obtaining a eq from (A.1) is straightforward and leads to very simple analytical expressions, but for e�0 (where e is defined in (A.7)), this would require solving a relatively complex non-linear equation. An alternative is to numerically determine a eq from (A.1), which is done in the present work using the Newton-Raphson method. Once that a eq is determined, it is substituted for a in expressions (A.5)-(A.8), which finally gives the fracture energy that will be denoted by G c:E . Note that when the equivalent crack length (a eq ) is used instead of the actual one (a), G c:E does not correspond to G c as defined in Equation (A.5) because d� da eq da eq da ¼ G c:E da eq da : (A:9) Determining da eq =da is impossible without measuring the actual crack length, but in Ref. [26] it was shown that the value of this derivative is usually very close to unity. Thus, G c � G c:E can be assumed whenever the crack-length measurement is not available without any significant loss in accuracy, which is also done in the present work. In Figure A13, the comparison between the values of G c computed using the approach that takes into account the loading block (expressions (A.5)-(A.8)) and the one that does not (G c ¼ F 2 a 2 eq =ðb EIÞ) is shown. The average load-displacement data for the load-line displacement speed of 0.1 mm/min (plotted in Figure 8) is used in the analysis. The comparison shows that neglecting the effect of the loading blocks leads to slightly larger values of the fracture energy, with a relative error with respect to the approach that does take it into account of approximately 3%. Note that the main difference between the R-curves plotted in Figure A13 can be attributed not to the difference in the fracture energy but to the horizontal shift that occurs because for the same values of F and δ, a eq is different depending on whether the loading blocks are taken into account or not in Equation (A.1). This implies that for the present set of experiments, the effect of the loading blocks can be neglected without a significant influence on the results and final conclusions. Note that in Section 4, this analytical model (without loading blocks) is further extended to account for shear strains in DCB arms. Figure A13. A comparison between the R-curves obtained using formulae for G c that either take into account or neglect the influence of the loading blocks using the average load-displacement data obtained for tests performed at 0.1 mm/min.