# Steady-state crack propagation during DCB test with viscoelastic-viscoplastic interface – Correlation between rheological behavior and effect of crack propagation rate

Julien Jumel

► To cite this version:

Julien Jumel. Steady-state crack propagation during DCB test with viscoelastic-viscoplastic interface – Correlation between rheological behavior and effect of crack propagation rate. International Journal of Solids and Structures, 2024, 296, pp.112824. 10.1016/j.ijsolstr.2024.112824. hal-05340791

HAL Id: hal-05340791

https://hal.univ-brest.fr/hal-05340791v1

Submitted on 31 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Steady-state crack propagation during DCB test with viscoelastic-viscoplastic interface
## Correlation between rheological behavior and effect of crack propagation rate

Julien JUMEL

IRDL - ENSTA Bretagne, 2 rue François Verny, 29200, Brest.

**Abstract:**

Adhesively bonded joints are known to suffer from creep and relaxation phenomena due to the viscous nature of the adhesive. These effects are detrimental to the durability of structural joints as they can lead to progressive and slow crack propagation rate phenomena and possible delayed failure. Master curve approaches are often used to assess the durability of structural assemblies, but the physical justification for such approaches requires additional investigation to more accurately assess the coupling between structural effects and the possible complex long-term rheological behaviour of the adhesive layer. Following similar approaches developed to study the steady-state crack propagation regime in bulk materials, an Eulerian description of crack propagation along a viscous adhesive layer during a double cantilever beam test is implemented. Taking into account viscoelastic-viscoplastic and various hardening rules, some master curves can be obtained describing the evolution of the crack propagation rate as a function of the stationary loading conditions. With such a model, the role of adhesive layer behaviour and joint geometry on crack propagation conditions can be discussed, as well as the applicability of such a master curve approach.

**Keywords:**

Steady state, viscoplasticity, viscoelasticity, creep, subcritical crack growth, DCB

corresponding author: Julien.jumel@ensta-bretagne.fr
# 1. Introduction

The adhesive bonding technique is now frequently used to manufacture strong and durable structural assemblies in critical applications that also require high levels of reliability. However, adhesively bonded joints still suffer from limited performance when exposed to aggressive environmental conditions with elevated temperatures and humid atmospheres (Broughton, 1999). A drastic reduction in mechanical performance can then be observed due to temperature activation and plasticization effects. Possible irreversible chemical degradation of adhesive and adherend/adhesive interface may further exacerbate these detrimental effects (Armstrong, 1997; Adams, 2009). In this context, the assessment of the long-term durability of bonded joints is challenging, but crucial to avoid the catastrophic failures observed in the past. Stationary in-service loading conditions could be investigated in detail as they are encountered in many situations such as anchorages or repair patches and can lead to dramatic delayed creep failures (McNichol, 2008; Siddika, 2020).

Two approaches can be used to ensure a safe service life. The first, called "safe life", consists in dimensioning the elements so that no decohesion is susceptible to initiate under the influence of operating conditions (Djeumen, 2020; Mizah, 2016; Sahellie, 2015). The second, called "damage tolerance", consists of controlling the development of damage, cracking, delamination, or debonding so that they do not reach a critical extension that could lead to the collapse of the structure (Nakamura, 2023; Ogawa, 2022; Tanulia, 2022). Both approaches are based on the use of test data, supplemented by numerical modeling and simulation to assess more complex scenarios. In the case of bonded joints and delayed failure analysis, creep and relaxation tests on technological or more specific specimens are mainly used. As part of a fail-safe approach, master curves indicating the time to failure or decohesion rate are usually obtained. The specimens used are SLJ, butt, TAST, etc., so that the stress distribution is quasi homogeneous along the joint plane and the overall failure of the assembly is mainly driven by
the initiation of the first decohesion. As part of a damage tolerance approach, crack propagation tests are carried out under stationary loading conditions to obtain “master curves” where the crack propagation rate is reported as a function of the mechanical parameters (Jhin, 2013; Budzik M. J., 2009; Moreira, 2022) which are representative of the loading seen by the crack (G, C<sub>v</sub> (Park, 2006), C* (Elmukashfi, 2017; Ji, 2004), mode mixity (Yang, 2022), etc...). Once the crack is initiated, the service life of the assembly is controlled by the duration of the crack propagation regime.

The use of master curves for designing structures which are significantly different from test specimens is questionable due to the dependence of the mechanical state within the adhesive layer on the geometry of assembled parts and loading conditions. In addition, the analysis of crack propagation condition is often based on a linear elastic or nonlinear framework which may ignore the progressive deformations due to the viscoplastic and/or viscoelastic nature of the adhesive or the one of the adherends (Škec, 2023; Groth, 1990; Schapery, 2022; Nuwayer, 2018). In order to more accurately assess the contribution required for damage development at the crack tip and the energy required for propagation, complementary analyses are required to assess the overall energy balance. This specific aspect is the focus of the present study. The focus is on the steady state case.

While many modeling frameworks exist for the study of crack propagation (Hui, 2022), the focus on delamination and decohesion of bonded joints and DCB test has been limited (Jirásek, 2019; Costa, 2023; Ferracin, 2003; Li, 2016; Mizuno, 2005; Elmukashfi, 2017; Wang, 2016). This particular regime is characterized by a constant propagation velocity and stationary mechanical fields near the crack front, allowing for the use of an Eulerian description. Using this approach, the shape of the mechanical fields surrounding the crack can be evaluated, considering different material behaviors, without simulating the initiation phase.
In addition, as the mesh remains fixed, there is no need for remeshing techniques, extended finite element method (XFEM), or complex phase fields to describe damage evolution while avoiding mesh sensitivity. This approach has primarily been used to study crack propagation in solid materials, including dissipative mechanisms such as plasticity, viscoelasticity, viscoplasticity, damage, and more complex behaviors such as hyperelasticity and cohesive interfaces (Juul, 2019).

The configuration studied here therefore complements the analyses proposed elsewhere, by focusing on the case of DCB type Mode I crack propagation test, used to study delamination in composite materials or adhesive layer decohesion in bonded joints. We are interested here in analyzing the influence of the crack propagation on the fracture process zone configuration (Budzik M. K., 2011) rate and identifying the particular loading conditions that lead to a slow rate crack propagation regime. The DCB configuration is specific because the evolution of strain and stress along the adhesive layer results from the interaction between the deformation of the substrate and that of the adhesive. The substrates are generally slender and can be represented as flexible beams, while the mechanical behavior of the thin adhesive layer is represented by a zero-thickness cohesive interface, which exhibits complex time-dependent rheological behavior.

We then present the computational methodology used to determine the evolution of strain and stress within the adhesive along the joint plane during steady-state decohesion propagation under constant bending moment loading condition (Plausinis, 1995; Wang, 2016). As in many theoretical analyses of delamination tests or crack propagation tests along bonded joints (Szekrényes, 2023; Azab, 2022), a purely one-dimensional analysis is proposed to help evaluate the contribution of viscous parameters to the dependence of $G_c$ on the crack propagation rate dependence. This could later lead to finite element implementations which are needed to investigate the failure of more complex structures or loading conditions (Abdel-
Monsef, 2023; Riccio, 2023). The behavior of the adhesive layer is assumed to be viscoelastic, viscoplastic and damaging (Mohapatra, 2021; Rocha, 2019; Chen, 2023; Chen, 2021). A complex relaxation spectrum is introduced, as well as a strain-hardening law to study the influence of the isotropic/kinematic nature. The effect of each viscoelastic/viscoplastic contribution is studied separately before coupling situations are considered for a wide range of crack propagation rates. From these analyses, master curves are simulated which describe the evolution of crack propagation rate as a function of macroscopically applied loading conditions, and demonstrate the link with local rheological behavior.

In the following, the classical analysis method for DCB test is first presented in a linear elastic framework. It emphasizes the importance of a steady state crack propagation condition to extend this approach to more complex interface behavior. A viscoelastic-viscoplastic cohesive zone model is then constructed from the uniaxial response of a virtual strain rate dependent polymer material. The steady-state hypothesis and Eulerian numerical resolution scheme for crack propagation during DCB test is then introduced. Theoretical evolutions of critical SERR vs crack propagation rate are then obtained for various combinations of the model parameters values, which emphasizes the role of the main component of the rheological model particularly on the critical speed transition between slow and fast crack propagation regimes.

## 2. Reference solutions – Simple beam theory vs. Beam on elastic foundation.
<figure><img src="image_2.png" /><figcaption>Figure 1. Beam on cohesive interface representation of a DCB specimen</figcaption></figure>

## 2.1 Simple beam theory (SBT)

The double cantilever beam (DCB) test (see Figure 1.) is a popular method for evaluating the toughness of bonded joints or composite delamination under Mode I loading conditions. Bonded specimens are obtained by joining two flexible slabs (width: $w$, thickness: $t$) with a layer of adhesive applied on one side over a distance $L$. A length, $a$, is left unbonded along the other side of the specimen, at the edge of which a progressive opening relative displacement, $\Delta$, is applied between the two bonds. During the test, this cleavage loading condition initiates and propagates a crack along the bond line. Two opposing forces, $P$, are applied to each adherend and recorded during the test. From the $P(\Delta)$ evolution recorded during the test, an energy balance analysis is performed to estimate the strain energy release rate (SERR), $G$, that drives crack initiation and propagation during the test. The SERR is obtained using the Irwin-Kies relation:

$$G = \frac{1}{2w} \frac{\partial C(a)}{\partial a} P^2 \quad (1)$$

where, $C(a) = \Delta/P$, is the specimen compliance and $a$ is the instantaneous crack length. For the Irwin-Kies relation to be valid, the entire specimen should behave as a linear elastic body. A robust $G$ evaluation then requires both accurate monitoring of crack propagation and adequate evaluation of specimen compliance. Due to the geometry of the specimen, the
simplest DCB model assumes that the two adherends behave as two flexible beams under simple bending conditions. It is then assumed that the effect of bond line compliance and beam shear deformation is negligible. The compliance is approximated by:

$$C = \frac{\Delta}{P} = \frac{2a^3}{3EI} \qquad (2)$$

So that:

$$G = \frac{a^2 P^2}{wEI} = \frac{M^2}{wEI} \qquad (3)$$

where *M* is the bending moment at crack tip position. The compliance of the adhesive layer is generally not negligible and a plastic zone may develop in front of the crack. To account for such effects, crack length correction and/or effective crack length methods are introduced where the geometric crack length, *a*, is replaced by a corrected or effective value, *a*<sub>eff</sub> = *a* + *a*<sub>c</sub>, where *a*<sub>c</sub> is the crack length correction. *a*<sub>c</sub> is determined either from modeling or experimental calibration. The effective crack length is also estimated directly from the instantaneous specimen compliance using the relation:

$$a_{eff} = \sqrt[3]{\frac{3EI\Delta}{2P}} \qquad (4)$$

During the crack propagation regime, the failure criterion *G* = *G*<sub>c</sub> is generally introduced, where *G*<sub>c</sub> is the critical SERR. This condition is introduced when a steady-state crack propagation regime is observed, i.e. the dissipative mechanisms involved in the damage and cracking phenomena remain stationary. Then the crack length correction value and the *G*<sub>c</sub>
value should remain the same. Before introducing the complex rheological behavior of adhesive bonds, simple elastic behavior should first be considered to validate the steady-state hypothesis.

## 2.2 Beam on elastic foundation analysis

Extensive research has been carried out to evaluate the evolution of stress and strain along the adhesive layer during DCB crack propagation tests. This has been done using either analytical or numerical approaches. Very few of these have specifically addressed the analysis of the steady-state crack propagation regime during the DCB test along the viscous interface. In all contexts, the main governing equations remain essentially the same. The beam on elastic foundation theories are the simplest to introduce the effect of bond line compliance and the presence of a finite fracture process zone (FPZ) size ahead of the crack tip position where large stress and strain gradients develop (see Appendix 1). The process zone extension is mainly controlled by a characteristic wave number:

$$ \lambda_{\sigma} = \frac{\sqrt{2}}{2} \left( \frac{wk_{a}}{EI} \right)^{1/4} \qquad (5) $$

where EI is the bond flexural stiffness, w is the specimen width, and $k_a=\sigma/\Delta\nu$ is the bond line compliance. In this case, the flexural stiffness of the bond is expected to contribute mainly to the specimen deformation in the FPZ region. However, the shear deformation of the bond may also contribute significantly to the evolution of the bond deformation along the FPZ. A second characteristic wavenumber should then be introduced:
$$ \lambda_{\tau} = \sqrt{\frac{wk_{a}}{\kappa GS}} \qquad (6) $$

with $\kappa GS$ being the adherend shear rigidity. This parameter controls the extension of the FPZ when adherends behave as shear deformable beam.

<table><caption>**Table 1.** DCB simulation parameters – elastic case</caption><thead><tr><th>T</th><th>W</th><th>κ</th><th>E</th><th>ν</th><th>t<sub>a</sub></th></tr></thead><tbody><tr><td>5mm</td><td>25mm</td><td>5/6</td><td>70 GPa</td><td>0,3</td><td>200μm</td></tr></tbody></table>

Figure 2 shows the evolutions of the adherend cross-section deflection and rotation along the adhesive bondline as well as those of the cohesive forces (shear force $T$ and bending moment $M$) for three different adhesive layer stiffness values, $k_a = E_a/t_a$, and considering different dimensionless crack length values, $\lambda_{\sigma a}$. All other parameters are given in Table 1. These evolutions are obtained by using both analytical solutions, second-order and fourth-order finite difference numerical resolution schemes to validate the accuracy of the lasts. The details of the finite difference implementation can be found in other contributions and specialized numerical textbook (Jones, 1997; Jumel, 2017). Both second order and fourth order numerical schemes have been implemented with no noticeable difference. Then, the second order integration scheme will be used in the following. For $L \cdot \lambda_{\sigma} > 7$, the results become independent of the left end boundary conditions (clamped or free end) due to the strong exponential attenuation of the mechanical fields from the bondline end. Then, when the number of integration points is greater than $N = 2000$, the numerical finite difference calculations are almost indistinguishable from the analytical ones.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Evolution of normalized adherend deflection (a), cross section rotation (b) and cohesive forces (shear force (c), bending moment (d)) along the adhesive layer.</figcaption></figure>

In addition, appropriate scaling coefficients are applied to these. Firstly, all evolutions are plotted as a function of $\lambda_{\sigma}x$ rather than $x$ demonstrating that the attenuation is effectively controlled by $\lambda_{\sigma}$. The deflection, $v$, is first normalized by the maximum deflection at the edge of the bondline, $v_{max}$. Then, all simulation results lead to mainly 3 evolutions corresponding to respectively short / zero crack length, intermediate crack length value, and long crack length value / pure bending loading condition. In fact, for $\lambda_{\sigma}a > 10$, all evolutions are almost identical and become independent of the dimensionless crack length $\lambda_{\sigma}a$. Assuming that the adherend shear deformation is negligible, $dv/dx \approx \varphi$. Then, again assuming $v_{max} = 1$, the evolution of the rotation of the adherend cross section is multiplied by $1/\lambda_{\sigma}$, again showing the same three regimes. However, a small discrepancy is observed due to the remaining small contribution of shear deformation. Considering $v_{max} = 1$, the bending moment evolution is again normalized by $1/\sqrt{E_a}$. In fact, for the short or long crack length conditions, the SERR is
proportional to $P^2$. Considering linear elastic adhesive bondline behavior $G_c = E_a/2 v_{max}^2/t_a$ so that $1/\sqrt{E_a}$ scaling is appropriate to obtain the same normalized $G_c$ value and corresponding applied critical bending moment. Finally, after applying the same normalization, the shear force evolution is normalized by $1/\lambda_\sigma$ reflecting the effect of relation $T = -dT/dx$ again.

The following important conclusions can be drawn from this preliminary elastic analysis. The extent of the process zone increases as the stiffness of the adhesive layer decreases. The same spatial evolution of static quantities (i.e. $M$, $T$, $\sigma$) and kinematic quantities (i.e. $v$, $\phi$) are found when appropriate scaling is performed, but considering the same dimensionless crack length value $a.\lambda_\sigma$. This parameter generally cannot be kept constant during the DCB test, and transitions from intermediate to long crack lengths are expected, which means that a crack length correction should be introduced for a proper estimation of $G$. However, the present paper focuses on the role of the adhesive layer rheology on $G_c$, then a pure bending loading condition will only be considered in the following.

### 3. Viscoelastic-viscoplastic behavior:

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Viscoelastic-viscoplastic rheological model</figcaption></figure>

Even when adherends exhibit linear elastic behavior only, crack propagation rate dependence and crack propagation under steady-state loading conditions are observed because adhesives
generally exhibit viscous behavior. A more complex viscoelastic-viscoplastic rheological behavior, which exhibits creep and relaxation phenomena, is then introduced below and schematically illustrated in Figure 3. This rheological model is a more sophisticated version of the Nishihara model. Formulations applicable to three-dimensional media can be found in the literature (Ilioni, 2018). The model uses a classical partition of the total deformation into a serial viscoelastic and viscoplastic contributions. The details of these behaviors are described in the following sections. The values of the material parameters are chosen close to those identified in (Ilioni, 2018), but are adapted to obtain macroscopic responses that are both physically realistic and representative of those that lead to delayed failure phenomena and, consequently to durability issues. Finally, it should be remembered that viscous effects are thermally activated. All processes are assumed to be isothermal here regardless of the crack propagation rate value. Possible thermomechanical couplings are worth investigating in the future.

### 3.1 viscoelastic behavior

A simple linear non-ageing linear viscoelastic behavior is introduced. The viscoelastic response of the adhesive is then governed by its relaxation function, $\mathcal{R}(t)$, so that the peel stress evolution, $\sigma$, can be determined from the applied strain time evolution using the Boltzmann superposition principle:

$$ \sigma(t) = \int_{0}^{t} \mathcal{R}(t-u) \frac{d\varepsilon_{ve}(u)}{du} du \qquad (7) $$

Various viscoelastic modeling approaches are applicable to provide an analytical representation of the relaxation spectra (Shitikova, 2022), (Yusoff, 2011) either in the time or
frequency domain. In the following, a classical approach is chosen which consists in approximating the continuous relaxation function with a discrete sum of exponential functions (i.e. Prony series):

$$ \mathcal{R}(t) = E_{\infty} + \sum_{i=1}^{N} E_i \exp(-t/\tau_{ri}) \quad (8) $$

without losing the generality regarding the application of the global approach. A single time constant would be sufficient to introduce some viscoelastic response and the drop in stiffness in the slow deformation rate regime due to the transition between glassy and rubbery behavior. Nevertheless, additional terms can be added to introduce secondary transitions but also to refine the description of each of them. In fact, not only the characteristic time of each transition should be introduced but also its width. A Gaussian distribution of the moduli around the central characteristic time is then introduced:

$$ E_i = E_{norm} \left[ \exp \left( - \left( \frac{\log_{10} \tau_i - \log_{10} \tau_c}{\log_{10}(\tau_c) \cdot \Delta \tau} \right)^2 \right) \right]^b \quad (9) $$

$$ \tau_{ri} = 10^{\left(\frac{\log_{10}\tau_{max} - \log_{10}\tau_{min}}{N-1}\right)^{\frac{i-1}{N-1}} + \log_{10}\tau_{min}} \quad (10) $$

$$ \tau_{min} = \tau_c 10^{-k\Delta\tau} \qquad \tau_{max} = \tau_c 10^{k\Delta\tau} \quad (11) $$

$\tau_c$ is the central characteristic time associated with a given transition, $\Delta\tau$ controls the width of the transition and $b$ is an additional tunable shape parameter. The values of $E_{norm}$, $k$ and $N$ values are chosen to give sufficient precision in the description of the relaxation spectrum and adequate rubbery, $E_{\infty}$, and glassy, $E_0$, moduli values.
A controlled strain loading condition, $\epsilon_{tot}(t)$, is generally applied to the adhesive from which the viscoelastic contribution, $\epsilon_{ve}(t)$, is determined. The applied stress, $\sigma$, is then determined numerically using the Boltzmann superposition principle (see relation (7)). Since the viscoelastic deformation is known for a discrete set of time values, $\epsilon_{ve}(t_k)$, using evenly spaced time intervals, $\Delta t$, a recursive relation:

$$ \sigma(t_{k:N}) = \sigma(t_{k:N}) + E_{\infty}(\epsilon_k - \epsilon_{k-1}) + \frac{\epsilon_k - \epsilon_{k-1}}{t_k - t_{k-1}} \sum_{i=1}^{N} \tau_{ri} E_i[1 - \exp(-(t_k - t_{k-1})/\tau_{ri})] \exp(-(t - t_k)/\tau_{ri}) \quad (12) $$

is used to determine the applied stress at each time step assuming that the viscoelastic strain evolution is approximated by a series of continuous segments.

To evaluate the creep response of the material under stationary loading condition, the retardation function can be evaluated from the discrete relaxation spectrum. The retardation function is again described with a Prony decomposition according to the relation:

$$ j(t) = j_{\infty} - \sum_{i=1}^{N} j_i \exp(-t/\tau_{ci}) \quad (13) $$

$j_i$ and $\tau_{ci}$ are the creep compliances and characteristic creep times of the discrete retardation spectrum, respectively. The method described in (Loy, 2015) is implemented here. First, the creep time constants are obtained by finding the zeros of the relaxation function in the Laplace domain. Then the compliance, $j_i$, are obtained using the relations:
$$j_i = \left( \sum_{k=0}^{N} \frac{E_k \beta_i^2}{(\alpha_k - \beta_i)^2} \right)^{-1} \qquad (14)$$

with:

$$\alpha_k = 1/\tau_{rk} \qquad \beta_k = 1/\tau_{ck} \qquad (15)$$

For the analysis presented in the following sections, reference viscoelastic parameters are chosen whose values are given in Table 2. The resulting relaxation and retardation functions and spectra are shown in Figure 2, as is the theoretical stress-strain response under constant strain rate conditions.

<table><caption>**Table 2.** Viscoelastic relaxation spectrum characteristic parameters</caption><thead><tr><th>E<sub>∞</sub></th><th>E<sub>0</sub></th><th>τ<sub>c</sub></th><th>Δτ</th><th>b</th><th>k</th><th>N</th></tr></thead><tbody><tr><td>500MPa</td><td>2GPa</td><td>10000s</td><td>0,6</td><td>2</td><td>6</td><td>21</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** reference viscoelastic behavior, (a) relaxation & retardation functions, (b) relaxation & retardation spectra, (c) constant strain rate behavior.</figcaption></figure>
With the chosen parameters, a large rubbery to glassy transition is observed over 6 decades, both in terms of creep/relaxation time and strain rate loading conditions.

## 3.2 Viscoplastic behavior

A Norton's law is used to model the viscoplastic contribution, $\epsilon_{vp}$, to the total deformation, as described by equation:

$$ \frac{d\epsilon_{vp}}{dt} = \mathrm{sign}(\sigma) \left( \frac{|\sigma - \sigma_y|}{K_{vp}} \right)^{n_{vp}} \qquad |\sigma - \sigma_y| > 0 \qquad (16) $$

$$ \frac{d\epsilon_{vp}}{dt} = 0 \qquad |\sigma - \sigma_y| < 0 \qquad (17) $$

$K_{vp}$ and $n_{vp}$ are material parameters and $\sigma_y$ is the instantaneous viscoplastic yield stress. A hardening rule is also introduced, to modify the viscoplastic yield stress value according to the time evolution of the viscoplastic deformation. Considering an isotropic hardening behavior, the evolution of the yield stress is described by the relation:

$$ \sigma_y = \sigma_y^0 + H \epsilon_{vp}^{cum} + Q[1 - \exp(-b \epsilon_{vp}^{cum})] \qquad (18) $$

where:

$$ \epsilon_{vp}^{cum} = \int_{0}^{t} |d\epsilon_{vp}| \qquad (19) $$

is the cumulative plastic stain. As can be seen in Figure 2, the load path to which the adhesive is subjected alternates between compression and tension. An alternative type of hardening
behaviour should then be considered. First, assuming kinematic hardening, the following relations are used to evaluate the viscoplastic flow. Norton's law is modified as follow:

$$ \frac{d\varepsilon_{vp}}{dt} = \mathrm{sign}(\sigma) \left( \frac{|\sigma - X| - \sigma_y^0}{K_{vp}} \right)^{n_{vp}} \qquad |\sigma - X| - \sigma_y^0 > 0 \qquad (20) $$

$$ \frac{d\varepsilon_{vp}}{dt} = 0 \qquad |\sigma - X| - \sigma_y^0 < 0 \qquad (21) $$

The evolution of the kinematic hardening parameter, X, is obtained by incrementing the relations along the loading sequence:

$$ X = X_1 + X_2 \qquad (22) $$

$$ X_1 = H\varepsilon_{vp} \qquad (23) $$

$$ dX_2 = bQd\varepsilon_{vp} - bX_2d\varepsilon_{vp}^{cum} \qquad (24) $$

Finally, due to the confinement of the adhesive layer, a strong hydrostatic component is expected in the adhesive layer along the length sustaining negative peel stress. A dissymmetric plastic behavior is sometime observed, also introduced in more traditional cohesive zone models which consists in blocking the irreversible strain when $\sigma<0$. Three hardening behaviors are then considered in the following sections (i.e. isotropic, kinematic, compression blocked). The viscoplastic parameters given in Table 3 will be used as a reference configuration for most of the simulations presented in the following. These “virtual material” parameters have been chosen because the mechanical responses obtained are consistent with those observed for some representative adhesives, as illustrated below.

<table><tr><td>K<sub>vp</sub></td><td>n<sub>vp</sub></td><td>σ<sub>y</sub><sup>0</sup></td><td>H</td><td>Q</td><td>n</td></tr></table>
<table><caption>**Table 3.** Reference viscoplastic parameters</caption><tr><td>45GPa</td><td>2</td><td>10MPa</td><td>4MPa</td><td>2MPa</td><td>50</td></tr></table>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Stress/strain evolution of viscoelastic/viscoplastic material under slow strain rate condition ($\dot{\epsilon} = 10^{-7}$) under (a) pure tension and (b) compression-tension loading conditions</figcaption></figure>

The effect of three types of hardening is illustrated in Figure 5 where are compared the stress-strain curves are compared for monotonic tension and alternating compression-tension loading conditions. The strain rate is chosen to be significantly low so that the adhesive is in the rubbery state and no gap is observed between the applied stress and the yield stress in the viscoplastic regime. As expected, while there is no difference for the monotonic case, hysteresis is observed when plasticity isn't locked under compressive loading. Also, the onset of plasticity and then the level of the plastic plateau is slightly different depending on the type of hardening behavior, which may change the global energy dissipated during the fracture process in the analysis of the DCB test.

### 3.3 Damage-failure criteria
Crack propagation is generally driven at a the macroscopic scale by energetic conditions. At the bondline scale either ultimate stress or ultimate strain failure criteria could be introduced. Both conditions are considered below. First, by limiting the peel stress value with the condition $\sigma < \sigma_{max}$. Second, by introducing a progressive damage function so that the effective peel stress applied to the adherend is given by the relation:

$$ \sigma_{eff} = (1 - D)\sigma \qquad (25) $$

with:

$$ D = 0 \qquad \varepsilon < \varepsilon_{ini} \qquad (26) $$

$$ D = \frac{\varepsilon_{max} \quad \varepsilon - \varepsilon_{ini}}{\varepsilon \quad \varepsilon_{max} - \varepsilon_{ini}} \qquad \varepsilon < \varepsilon_{max} \qquad (27) $$

When considering elastic damage behavior, the classic triangular interface separation law is obtained. For the reference configuration, damage initiation occurs at $\varepsilon = \varepsilon_{ini} = 40\%$ total deformation value, complete failure is set for $\varepsilon = \varepsilon_{rupt} = 50\%$ total deformation value. A dual failure criterion is then introduced here. In fact, at high strain rates, the mechanical behavior often becomes brittle so that the elongation at break is greatly reduced due to limited plastic deformation. In order to reproduce this observed behavior macroscopically, an ultimate stress criterion is also considered so that complete failure is assumed above $\sigma = \sigma_{max} = 30MPa$.

### 3.4 Resulting viscoelastic-viscoplastic behaviour

A 15 parameter viscoelastic-viscoplastic-damage model is then introduced to represent the mechanical response of the adhesive under either strain controlled or creep loading conditions. An implicit iterative procedure is implemented to determine the stress evolution under prescribed strain evolution. The advantages and limitations of different numerical
schemes applicable to either viscoelastic, viscoplastic or coupled viscoelastic-viscoplastic
situations are discussed in details in many contribution dealing more specifically of this topic
especially in three-dimensional configuration. The numerical integration procedure proposed
here is found to be satisfactory for all loading conditions considered in this article and
adequate for the Eulerian description of the fracture processes.

Firstly, the uniaxial creep tension loading conditions are examined by simulating the
mechanical response after the sudden application of a constant stress in the interval from 14,5
MPa to 30 MPa until failure by means of standard explicit integration. Similar trends are
observed for all loading conditions but with different relative viscoelastic/viscoplastic
contributions and time scales. The resulting total deformation versus time evolution is shown
in Figure 6. The obtained power law evolutions are consistent with the one observed in other
contributions and generally fitted with Bailey-Norton equation (Majda, 2009; Peng, 2021).

As is generally observed in such tests, the time to failure increases by several orders of
magnitude as the applied stress decreases. With the chosen parameters, the time to failure is
consistent with the context of the present analysis, which aims to evaluate the delayed failure
state under long duration loading conditions.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** a) Creep response for various applied stress values, b) time to failure vs. applied stress</figcaption></figure>

<figure><img src="image_8.png" /></figure>
**Figure 7.** a) Stress versus strain evolution under constant strain condition, b) dissipated energy density until complete failure.

Finally, a series of simulations were performed to evaluate the influence of the strain rate
loading condition on the stress versus strain evolution and total energy dissipated prior failure
and assuming monotonic loading scenario. The results are presented in Figure 7. The effect of
viscoelasticity is quite negligible, but a slight variation of the apparent modulus during the
initial elastic loading regime can be observed. The impact of viscoplasticity is particularly
evident through the stress plateau increase with strain rate. The consequence of the dual
failure criteria is also seen on both stress versus strain evolutions and evolution of dissipated
energy with the strain rate. For low strain rate loading condition, progressive failure to
complete failure is observed, while avor a critical strain rate value, rapid failure is seen and
the elongation at break diminishes with strain rate. Two regimes are clearly established,
highlighting important implications regarding energy dissipation. Before reaching the ultimate
stress criterion, the dissipation increases due to viscoplasticity. However, after reaching the
ultimate stress criterion, a rapid decrease in the dissipated energy is observed due to loss of
ductility. This phenomenon has now been modelled at the material scale and the coupling
with the adherend has been introduced. The particular feature related to steady-state crack
propagation along viscoelastic-viscoplastic interface can be now be discussed.

4. Eulerian implementation of steady state crack propagation regime analysis
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Eulerian description of the steady state crack propagation regime.</figcaption></figure>

The analysis of the double cantilever test relies on a comprehensive energy balance analysis to determine the critical SERR needed to develop and propagate a crack. However, more advanced techniques are required to account for the more complex mechanical behavior of the adhesive bondline. Cohesive zone models are now broadly utilized as they allow adequate prediction of both the crack nucleation and propagation processes. Nevertheless, evaluating the sole critical SERR, $G_c$, value during the steady-state crack propagation regime using this method can be numerically costly since simulating the crack nucleation stage is required beforehand. Therefore, a Eulerian description of the steady-state crack propagation regime is proposed which aims to exploit the stationary properties of all mechanical fields in a mobile coordinate system attached to the crack tip position (See Figure 8) at reduced computation times.
3.1. Finite difference implementation of steady state analysis of DCB test

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Beam on elastic foundation reference static analysis - Finite difference implementation</figcaption></figure>

To calculate the beam deformation and the evolution strain and stress in the adhesive layer ahead of the crack tip position during the steady-state crack propagation regime of a DCB test, a Timoshensko beam on an elastic foundation analysis is used, as shown in Figure 9. The constitutive equations listed in Appendix 1 are solved using a finite difference resolution technique. Due to the non-linear behavior of the adhesive layer, an additional pressure evolution, $p(x)$, is introduced which corresponds to the difference between the predicted elastic stress versus beam deflection evolution and the actual stress.

<figure><img src="image_12.png" /><figcaption>**Figure 10.** iterative procedure for F.D. resolution of DCB steady state crack propagation regime.</figcaption></figure>
The double cantilever beam specimen is loaded by two opposite forces, P, which are located at a distance, a, away from the crack tip. Therefore, the boundary conditions applied at the crack tip position should be, $T(x = 0) = 0$ and $M(x = 0) = aP$. As explained in section 2, a pure bending loading condition is considered instead of discussing the relative effects of shear force and bending moment on the stress and strain evolution ahead of the crack tip position. This corresponds to a situation with large $a\lambda_\sigma$ values. Due to the strong attenuation of the peel stress and strain evolution ahead of the crack tip, the results of the numerical simulations become independent of the chosen left boundary conditions despite the necessary limitation of the grid length. Therefore, the clamped condition (i.e. $v(x = -L) = 0$; $\phi(x = -L) = 0$) was chosen.

Finally, the strain rate dependent nature of the adhesive layer should be considered in the analysis. As mentioned earlier, the static beam on elastic foundation analysis carried out using the F.D. resolution scheme corresponds to a steady-state regime in a moving coordinate system. Its velocity is equal to the crack propagation velocity. Therefore, the time evolution of the strain and stress can be obtained from the spatial evolution of the same quantities along the x-axis using the equation:

$$t(x) = (L + x) \cdot \frac{1}{da/dt} \qquad (28)$$

The time evolution of the stress applied to the bondline is determined using the relation:

$$\varepsilon(t) = 2 \frac{v(x(t))}{t_a} \qquad (29)$$
From the time evolution of the strain, the spatial evolution of the peel stress along the adhesive layer is determined by integrating the viscoelastic-viscoplastic damage law, as described in section 3. For the simulations presented in the following sections, a minimum length of $L = 35$mm is chosen for the FPZ extension and a minimum number of $N = 3000$ nodes is chosen which is sufficient to achieve both adequate temporal and spatial integration.

## 3.2. Analysis of steady state crack propagation regime

Before discussing the effect of rate-dependent adhesive layer behavior on the crack propagation regime, we introduce a basic elastic perfectly plastic behavior that illustrates the Eulerian and finite difference implementation of the steady-state crack propagation regime analysis. The results are shown in Figure 11.

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Description of the Fracture process zone in case of elastic perfectly plastic adhesive layer with damage.</figcaption></figure>
In this situation, the critical SERR, $G_c$, can be obtained by integrating the stress versus strain evolution along the adhesive layer with the relation:

$$G_c = t_a \int_{-\infty}^{0} \sigma(x) \frac{d\epsilon(u)}{du} du \qquad (30)$$

With the current elastoplastic parameters, $G_c = 945\text{J.m}^{-2}$. Finally using relation (3) the critical bending moment which leads to extension of the interface decohesion is obtained:

$$M_c = \sqrt{wEIG_c} \qquad (31)$$

$M_c = 20,75\text{Nm}$ is found as critical bending moment value which matches with the obtained applied bending moment at crack tip position. Since, no strain rate dependence is introduced at in the material behavior, the critical SERR value is independent on the crack propagation rate.
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Evolution of the adhesive stress, (a), and strain, (b), distribution in the FPZ for various crack propagation rate. Resulting stress vs strain evolution, (c), and critical SERR evolution vs crack propagation rate, (d). Compression locked hardening case.</figcaption></figure>

Now that the eulerian description has been introduced, the method is applied to determine the effect of the crack propagation rate on the critical bending moment and the critical SERR considering the full viscoelastic-viscoplastic-damage material behavior. First, the simplest hardening condition is considered (i.e. viscoplastic flow blocked under compressive loading conditions). The results are presented in Figure 12, focusing on the sole evolutions of the adhesive layer deformation and stress along the adhesive layer and the stress versus strain evolution. These evolutions are obtained for various crack propagation rates in between limits corresponding to strain rate independent fully relaxed materials and glassy brittle behaviors. A clear peak resistance to crack propagation is observed, marking the transition between the slow crack propagation rate regime and the brittle regime. Such peak resistance has been observed in slightly different contexts such as peel tests (Derail C. A., 1997; Derail C. A.,
1998) where viscous dissipative processes predominate. Prior to the peak $G_c$ value, a regular increase in the critical SERR as a function of the crack propagation rate is observed. As observed in Figure 12, this effect is mainly controlled by the increase in the apparent plastic plateau as the crack propagation rate increases. However, saturation is observed after the peak $G_c$ value position due to the decrease in elongation at break as the maximum stress failure criterion is reached. At the same time, the apparent elastic modulus increases due to viscosity which reduces the $G_c$ value again.

# 4. Influence of bonded layer rheological parameters on crack propagation conditions

A total of 15 parameters are necessary to depict the complete viscoelastic-viscoplastic-damage mechanical behavior of the adhesive layer and then to evaluate the $G_c(du/dt)$ master curve. The relative contribution of all these effects is qualitatively demonstrated by the strain and stress distribution along the FPZ and the stress-strain evolution. However, for better understanding of the specific contribution of viscoelastic and viscoplastic processes, these situations are analyzed separately.

## 4.1 Viscoelastic contribution

In order to evaluate the sole influence of the viscoelastic behavior of the adhesive layer, an incommensurate yield stress value is chosen so that the failure criterion is reached before any viscoplastic deformation occurs. When considering both the maximum stress and progressive damage criterion, the ultimate stress associated with brittle failure is observed in all cases before progressive damage could occur, so this situation is presented first. Nevertheless, two distinct damage criteria are examined, maximum elongation at break first, then ultimate stress criteria.
A series of simulations are performed to evolution of adhesive stress and strain evolution along the FPZ as a function of crack propagation rate considering purely viscoelastic behavior and assuming 1.5% maximum peel strain is reached at the bondline edge position. The critical SERR versus crack propagation rate evolution is then determined from the critical bending moment applied at bondline edge position using relation (3), $G_c = M_c^2/wEI$. The viscoelastic relaxation spectrum is mainly governed by four parameters (i.e. glassy modulus, $E_0$, rubbery modulus, $E_\infty$, central time constant $\tau_c$, relaxation peak width $\Delta\tau$). Same glassy and rubbery modulus have been considered but distinct $\tau_c$ and $\Delta\tau$ values have been considered to evaluate the effect on the slow to fast crack propagation rate conditions transition.

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Critical SERR as a function of normalized crack propagation rate, viscoelastic adhesive bondline behaviour, maximum elongation at break criteria.</figcaption></figure>

The evolutions of critical SERR as a function of normalized crack propagation rate are presented in Figure 13. When normalized crack propagation rate, $\tau_c \cdot da/dt$, is used all evolutions obtained for a same value of $\log_{10}(\tau_c) \cdot \Delta\tau$ are found the same. A sigmoidal function
is obtained whose upper and lower $G_c$ limits are simply obtained by considering linear elastic behavior in the glassy and rubbery states respectively. The normalization indicates that crack propagation rate which governs the transition between slow fast crack propagation regimes is controlled by the characteristic time of the relaxation spectrum, $\tau_c$. Similarly, the width of the transition is controlled by $log_{10}(\tau_c) \cdot \Delta\tau$, that is the larger the peak of the relaxation spectrum, the larger the transition between slow to fast crack propagation regime. Interestingly, it should be observed that the presence of a narrow peak in the relaxation spectrum produces a peak $G_c$ value prior the fast crack propagation regime is reached.

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Stress versus strain evolution for various crack propagation rate. a) narrow viscoelastic relaxation peak, b) wide viscoelastic relaxation peak.</figcaption></figure>

The origin of the peak $G_c$ value can be interpreted by analyzing the evolution local peel stress versus peel strain during the crack propagation until failure. The consequence of the viscoelastic behavior is predominantly observed through the stiffening of the adhesive layer as crack propagation rate increased leading to higher local deformation rate. Additionally, due to alternating compression and tension loading the viscoelastic behavior produces and hysteretic response marking the presence of viscoelastic dissipation. Due no varying deformation rate along the bondline, the peak viscoelastic dissipation is not observed for intermediate crack propagation rate but prior fast crack propagation rate regime is reached.
Also the hysteretic behavior and then the viscoelastic dissipation is more pronounced for
narrow peak in the relaxation spectrum.

<figure><img src="image_17.png" /><figcaption>**Figure 15.** Critical SERR versus normalized crack propagation rate evolution. Viscoelastic behaviour, Ultimate stress failure criteria.</figcaption></figure>

To conclude with the analysis of the viscoelastic contribution, the same calculations were
used but now considering ultimate stress failure criteria σ = σmax = 30MPa. Since, linear
viscoelastic calculations are performed, a multiplicative factor is applied to convert the
ultimate strain criteria to ultimate stress criteria. Results are presented in Figure 15. Sigmoidal
evolution are found again. However, while considering same elongation at break an increase
of the critical SERR with increasing crack propagation rate is found, a decrease is found when
maximum stress criteria is considered due to loss of ductility and embrittlement. Also, is
should be noticed that monotonous evolution of Gc is found and no peak or minimum value is
found is this case. Choosing the appropriate criteria is then crucial to ensure accurate results,
as opposite outcomes can be found depending on whether the ultimate stress or ultimate
elongation at break criteria is selected. In “soft” systems, the SERR typically increases with
velocity which is observed when constant elongation at break criteria is used. However, such
criterion is not suitable for crosslinked rigid adhesives unless incommensurable strength is be
introduced in the model. These systems also exhibit embrittlement when facing rapid crack
propagation conditions. The decrease of G<sub>c</sub> would be due to a decrease of adhesive ductility,
not its strength. It is important to note that the sole viscoelastic component does not represent
all classes of behavior observed in adhesive systems. Therefore, it is necessary to investigate
viscoplastic behavior, as proposed below.

## 4.2 Viscoplastic contribution

We now consider an elastic viscoplastic damage behavior and assume no viscoelastic
contribution. Below the yield linear elastic behavior is considering and the rubbery modulus
E∞ = 500MPa is used. At first viscoplastic behavior with compression blocked viscoplastic
deformation is introduced. Kinematic and isotropic hardening will be introduced latter.

<figure><img src="image_18.png" /><figcaption>**Figure 16.** Evolution of the critical SERR as a function of the normalized crack propagation rate for various hardening exponent values.</figcaption></figure>
Numerous simulations are first performed to evaluate the evolution of the critical bending moment as a function of the crack propagation rate and considering distinct values of $K_{vp}$ modulus and $n_{vp}$ hardening exponent. From the critical bending moment, the critical SERR is determined using relation (3) again. Add to the elongation at break criteria, the damage rule introduced in section 3.3 is used here. The evolutions of the critical SERR as a function of normalized crack propagation rate, $(da/dt) \cdot (K_{vp}/\sigma_{ref})^{n_{vp}}$, are finally reported in Figure 16.

A series of sigmoidal master curves are found which are distinct for distinct $n_{vp}$ values. This indicates that the transition between slow to fast crack propagation is controlled by the characteristic time $(K_{vp}/\sigma_{ref})^{n_{vp}}$. The width of the transition increases with increasing hardening exponent value. In Figure 16, it should be noticed that all master curves intercept the same point. However, the location of the interception point depends on the value chosen for $\sigma_{ref}$. Here $\sigma_{ref} = 40MPa$ is chosen as an intermediate reference stress value between yield and ultimate stress value.
<figure><img src="image_19.png" /><figcaption>**Figure 17.** Stress versus strain evolution for various hardening exponent, $n_{vp}$, value a) for same critical SERR, $G_c$, value (Fig 16 - white square markers), b) at a same normalized crack propagation rate value (Fig 16 - grey square markers).</figcaption></figure>

To investigate the role of the Norton's law exponent on slow rate crack propagation condition the stress versus strain evolution obtained for different $n_{vp}$ values are compared for both same $G_c$ value (2000 J.m⁻²) and same normalized crack propagation rate value as indicated in Figure 16 with the two dotted lines as seen in Figure 17. For a same critical SERR, $G_c$, value (see white square markers in Figure 16) the average stress in the viscoplastic regime is then same but the slope of the stress is increasing with decreasing hardening exponent value. The interpretation is not straightforward, since this results from the coupling phenomena with adherend deformation. The comparison of stress versus strain evolution from a same crack propagation value (see grey square markers in Figure 16) is more straightforward. Since very similar evolution is observed here but an increase of the plastic stress plateau is observed for higher $n_{vp}$ value. This is a direct consequence of the Norton's viscoplastic flow rules. Indeed, for a given crack propagation rate, the deformation rate of the adhesive bondline edge is
mainly governed by the adhesive average bondline compliance in the fracture process zone as discuss in section 2. Then the deformation rate being essentially the same, an increase of $n_{vp}$ leads to increasing applied stress.

<figure><img src="image_20.png" /><figcaption>**Figure 18.** Influence of the viscoplastic strain hardening exponent, $n_{vp}$, on the peel stress evolution along the fracture process zone, a) for same critical SERR, $G_c$, value (see Fig 16 - white markers), b) at a same normalized crack propagation rate value (see Fig 16 - grey markers).</figcaption></figure>

Similar conclusions can be drawn by comparing the spatial evolution of the peel stress along the adhesive layer for different values of the Norton's law exponent and comparing the conditions that leads to the same $G_c$ and same $da/dt$ values respectively. These evolutions are plotted in Figure 18. They show that large variation of $da/dt$ are required to obtain the same $G_c$ value but varying $n_{vp}$ values, but the spatial evolution of the peel stress is left almost unchanged except near the crack tip position where the peel stress gradient increases with $n_{vp}$. For a given normalised $da/dt$ value, $G_c$ increases with $n_{vp}$. At the bondline scale, this is due to an overall increase in the peel stress variation amplitude with decreasing $n_{vp}$ value.
Meanwhile, the length of the plastic plateau ahead of the crack tip decreases. These descriptions are valid for the curves plotted here and for normalized crack propagation rate value smaller that observed for before the normalized characteristics speed smaller than the one of intersection point of all master curves. At this point, a more detailed discussion on the role of each of the model parameters would not lead to any firm conclusion at this stage. It is important to remember that the changes observed macroscopically on the master curves reflect reconfigurations of the spatial and temporal changes in the peel stresses and strains in the adhesive.

## 4.2 *Effect of hardening law.*

Finally, three different hardening behaviors have been introduced but not yet discussed. The contribution of the viscoplastic dissipation is mainly observed under slow crack propagation rate condition. As a consequence, the effect of hardening will be studied for normalized crack propagation rate ranging in the interval $[10^{-5} \ 10^{-1}]m$ as seen in Figure 19.
<figure><img src="image_21.png" /><figcaption>**Figure 19.** Critical SERR versus normalized crack propagation rate evolution. Combined Ultimate stress maximum elongation failure criteria. Effect of hardening rule.</figcaption></figure>

Considering combined maximum elongation at break and maximum stress failure criteria, similar master's curve evolutions are found for the three considered hardening rules (isotropic, kinematic, compression blocked) which are described in section 3. Only a slight offset is found for the three considered case before the peak $G_c$ value occurs. As expected, when viscoplastic flow is activated under compression loading condition, a higher $G_c$ value is found. Also a slightly higher value is found for the isotropic hardening than for the kinematic one.
<figure><img src="image_22.png" /><figcaption>**Figure 20.** a) stress and b) strain distribution at the crack tip vicinity, c) stress versus strain evolution. Effect of crack propagation rate and hardening law</figcaption></figure>

In Figure 20 are represented the peel stress and strain evolution ahead the crack tip for various crack propagation rate and considering the three different hardening rules. The resulting stress versus strain evolutions are also represented. When the viscoplastic flow is activated under compression loading conditions, an hysteretic loop is observed on the stress versus strain evolution which leads to higher dissipation and higher $G_c$ value. In case of kinematic hardening the dissipation effectively appears smaller than for the isotropic hardening due to decrease of the yield stress under tension loading condition after a compressive load is applied. The main consequence of the additional viscoplastic flow is seen on the stress distribution ahead the crack tip position. Indeed, an additional and rather large plastically deforming region is seen along the fracture process zone sustaining compressive peel stresses. This modification of the stress distribution slightly increases the overall process zone size length but also significantly increases the amplitude of the compressive peel deformation. For crack propagation rates larger than the one leading to peak $G_c$ value a rapid decrease of $G_c$ is observed for all three hardening laws due to activation of maximum stress criteria and decrease of elongation at break value. No significant effect of hardening law is then observed since the total viscoplastic flow in the compression region decreases as the crack propagation
rate decreases and the overall response tends to be elastic. The trends identified here are valid for the parameters of the viscoplastic model which used in this study. A more pronounced influence of the type of hardening may certainly be observed, for example, by modifying the fracture criteria parameters. More extensive parametric analysis could be undertaken in the future to investigate in details the role of the other parameters of the rheological model on “master curves” describing the effect of crack propagation rate on G_c.

## 5. Conclusion

In the present paper, the steady state propagation of a decohesion along a viscoelastic-viscoplastic-damaged adhesive layer is studied from a theoretical point of view, considering a Eulerian description of the crack propagation regime. The aim of the model is to establish a correspondence between the rheological behavior of the adhesive layer and the evolution of the critical SERR as a function of the crack propagation rate. After introducing the adhesive layer rheological model and the numerical resolution procedure simple rate-independent behavior situations are presented to introduce a pure bending moment loading configuration which is representative of the steady state crack propagation regime and facilitate the evaluation of the SERR. Then, by introducing a viscoelastic-viscoplastic-damage model with dual failure criterion, three different crack propagation regimes are observed. Under the slow crack propagation rate condition, the adhesive layer exhibits rubbery and viscoplastic behavior, and the critical SERR increases slowly with the crack growth rate due to increasing viscous dissipation. A critical SERR peak is observed when maximum stress is reached at the bond scale. The critical SERR then decreases with increasing crack propagation rate due to the embrittlement of the adhesive and loss of ductility. In order to clarify the role of the numerous adhesive rheological parameters, a parametric analysis is carried out, focusing first
on the viscoelastic and then on the viscoplastic behavior. It is observed that a one-to-one
correspondence can be found between the position and width of the slow and fast crack
propagation rate regime and the viscoelastic spectrum in the viscoelastic situation and
similarly with the Norton's law parameters. Finally, the influence of the hardening behavior
on the FPZ configuration is discussed. While, the energy balance is not significantly
modified, a large increase in the FPZ length is observed when both kinematic and isotropic
hardening are considered. Indeed, the potential plastic deformation in the compressive region
is rarely discussed, whereas it is clearly predominant here.

The proposed model and approach are very versatile and the parametric analysis could be
extended to evaluate a greater variety of situations and/or introduce additional refinements to
the adhesive layer behavior and/or to study specific situations. However, now that the
dependence of the crack propagation rate on G<sub>c</sub> related has been established in relation to the
strain rate dependence of the adhesive layer behavior, it is important to generalize the
evaluation of the SERR to other specimen loading conditions that are more representative of
DCB tests (cf. application of transverse force). Crack length correction or alternative SERR
calculation methods capable of accounting for the rate dependent nature of the system should
then be introduced. Similarly, possible rate-dependent behavior of the adherend could be
introduced. The steady-state analysis proposed above is a very versatile tool for analyzing
these different situations. At a more macroscopic scale, the crack propagation rate dependence
of the adhesive layer critical SERR is also an important effect to considered when analyzing
crack propagation tests that exhibit stick slip behavior but also to assess the resistance of
bonded joints and/or laminates under dynamics or impact loading conditions.
6. Références

Abdel-Monsef, S. T. (2023). Accurate simulation of delamination under mixed-mode loading using a multilinear cohesive law. *Engineering Fracture Mechanics*, 284, 109233.

Adams, R. D. (2009). The relative merits of the Boeing wedge test and the double cantilever beam test for assessing the durability of adhesively bonded joints, with particular reference to the use of fracture mechanics. *International Journal of Adhesion and Adhesives*, 29(6), 609-620.

Armstrong, K. B. (1997). Long-term durability in water of aluminium alloy adhesive joints bonded with epoxy adhesives. *International Journal of Adhesion and Adhesives*, 17(2), 89-105.

Azab, M. E. (2022). Influence of the DCB/Wedge beam kinematics on the identification of the cohesive parameters of interfaces. *International Journal of Solids and Structures*, 252, 111810.

Broughton, W. R. (1999). Creep testing of adhesive joints analysis of creep rupture data.

Budzik, M. J. (2009). Accurate and continuous adhesive fracture energy determination using an instrumented wedge test. *International Journal of Adhesion and Adhesives*, 29(7), 694-701.

Budzik, M. K. (2011). An in situ technique for the assessment of adhesive properties of a joint under load. *International Journal of Fracture*, 171, 111-124.

Chen, Y. S. (2021). Ratcheting and recovery of adhesively bonded joints under tensile cyclic loading. *Mechanics of Time-Dependent Materials*, 1-20.

Chen, Y. S. (2023). A viscoelastic-viscoplastic model for adhesives subjected to reversed cyclic load. *Mechanics of Time-Dependent Materials*, 1-17.

Costa, J. M. (2023). Influence of adhesive rheological behavior on stable mode II crack propagation in bonded interface subjected to shear loading condition. *Mechanics of Advanced Materials and Structures*, 1-15.

Derail, C. A. (1997). Relationship between viscoelastic and peeling properties of model adhesives. Part 1. Cohesive fracture. *The Journal of Adhesion*, 61(1-4), 123-157.

Derail, C. A. (1998). Relationship between viscoelastic and peeling properties of model adhesives. Part 2. The interfacial fracture domains. *The Journal of Adhesion*, 68(3-4), 203-228.

Djeumen, E. C. (2020). Creep investigations on adhesively bonded fasteners developed for offshore steel structures. *Marine Structures*. 69, 102660.

Elmukashfi, E. C. (2017). A theoretical and computational framework for studying creep crack growth. *International Journal of Fracture*, 208(1-2), 145-170.

Ferracin, T. L. (2003). On the determination of the cohesive zone properties of an adhesive layer from the analysis of the wedge-peel test. *International Journal of Solids and Structures*, 40(11), 2889-2904.

Groth, H. L. (1990). Viscoelastic and viscoplastic stress analysis of adhesive joints. *International Journal of Adhesion and Adhesives*, 10(3), 207-213.

Hui, C. Y. (2022). Steady state crack growth in viscoelastic solids: A comparative study. *Journal of the Mechanics and Physics of Solids*, 159, 104748.
Ilioni, A. B. (2018). A viscoelastic-viscoplastic model to describe creep and strain rate effects on the mechanical behaviour of adhesively-bonded assemblies. *International Journal of Adhesion and Adhesives*, 82, 184-195.

Jhin, G. A. (2013). Crack growth rate and crack path in adhesively bonded joints: Comparison of creep, fatigue and fracture. *International Journal of Adhesion and Adhesive, 46*, 74-84.

Ji, S. S. (2004). Polymer/metal interfacial crack growth characterized by C*. *International journal of fracture*, 129, 63-73.

Jirásek, M. A. (2019). Analysis of self-similar rate-dependent interfacial crack propagation in mode II. *International Journal of Fracture*, 220, 45-83.

Jones, G. (1997). *Analysis of beams on elastic foundations: using finite difference theory*. Thomas Telford.

Jumel, J. (2017). Crack propagation along interface having randomly fluctuating mechanical properties during DCB test finite difference implementation–Evaluation of Gc distribution with effective crack length technique. *Composites Part B: Engineering*, 116.

Juul, K. J.-P. (2019). Steady-state fracture toughness of elastic-plastic solids: Isotropic versus kinematic hardening. *Engineering Fracture Mechanics*, 207, 254-268.

Li, X. F. (2016). Fracture of a thin power-law nonlinear material with a crack using the DCB model. *International Journal of Fracture*, 201, 119-125.

Loy, R. J. (2015). Interconversion of Prony series for relaxation and creep. *Journal of Rheology*. 59(5), 1261-1270.

Majda, P. S. (2009). A modified creep model of epoxy adhesive at ambient temperature. *International Journal of Adhesion and Adhesives*, 29(4), 396-404.

McNichol, D. (2008). " The Big Dig" Unglued: NTSB Report on Ceiling Collapse Spreads the Blame. *Roads & Bridges*, 46(2).

Mizah, B. R. (2016). Novel method to measure the creep strength of adhesively bonded butt joints subjected to constant loading using a hydro-pneumatic testing machine. *The Journal of Adhesion*, 92(1), 65-79.

Mizuno, M. &. (2005). Simplified analysis of steady-state crack growth of piezoelectric ceramics based on the continuum damage mechanics. *Acta mechanica*, 179(3-4), 157-168.

Mohapatra, P. C. (2021). Adhesive hardening and plasticity in bonded joints. *International Journal of Adhesion and Adhesives*, 106, 102821.

Moreira, R. D. (2022). Mixed-mode I+ II fatigue/fracture characterisation of bi-material Aluminium/CFRP bonded joints. *Composites Part B: Engineering*, 246, 110240.

Nakamura, K. S. (2023). Creep Crack Growth Behavior during Hot Water Immersion of an Epoxy Adhesive Using a Spring-Loaded Double Cantilever Beam Test Method. *Materials*, 16(2), 607.

Nuwayer, H. M. (2018). Flexural creep behavior of adhesively bonded metal and composite laminates. *International Journal of Adhesion and Adhesives*, 84, 220-226.

Ogawa, Y. N. (2022). Evaluation of crack growth behaviors under Mode I static loading for two-part polyurethane adhesives. *International Journal of Adhesion and Adhesives*, 117, 103172.
Park, S. L. (2006). A nonlinear viscoelastic fracture analysis of concrete/FRP delamination in aggressive environment. *International journal of fracture*, **142**, 9-27.

Peng, Z. W. (2021). Static and sustained loading behavior of a basalt FRP shell-concrete composite bridge deck: An experimental and numerical study. *Engineering Structures*, **230**, 111689.

Plausinis, D. S. (1995). Application of a new constant G load-jig to creep crack growth in adhesive joints. *International journal of adhesion and adhesives*, **15**(4), 225-232.

Riccio, A. S. (2023). Numerical-experimental study on the compressive behaviour of repaired composite panels. *Composites Part B: Engineering*, **251**, 110466.

Rocha, I. B. (2019). Numerical/experimental study of the monotonic and cyclic viscoelastic/viscoplastic/fracture behavior of an epoxy resin. *International Journal of Solids and Structures*, **168**, 153-165.

Sahellie, S. P. (2015). Expectancy of the lifetime of bonded steel joints due to long-term shear loading. *Archives of civil and mechanical engineering*, **15**(4), 1061-1069.

Schapery, R. A. (2022). A theory of viscoelastic crack growth: revisited. *International Journal of Fracture*, **233**(1), 1-16.

Shitikova, M. V. (2022). Models of viscoelastic materials: a review on historical development and formulation. In I. P. Giorgio, *Theoretical Analyses, Computations, and Experiments of Multiscale Materials: A Tribute to Francesco dell'Isola*. (pp. 285-326). Springer Nature.

Siddika, A. A. (2020). Performances, challenges and opportunities in strengthening reinforced concrete structures by using FRPs-A state-of-the-art review. *Engineering Failure Analysis*. **111**, 104480.

Škec, L. &. (2023). Experimental and numerical study of rate-dependent mode-I failure of a structural adhesive. *The Journal of Adhesion*, **99**(8), 1323-1355.

Szekrényes, A. (2023). Fracture mechanical stability of DCB and 4ENF tests complemented with linear springs. *Theoretical and Applied Fracture Mechanics*, **125**, 103900.

Tanulia, V. W. (2022). Experimental and computational assessment of disbond growth and fatigue life of bonded joints and patch repairs for primary airframe structures. *International Journal of Fatigue*, **159**, 106776.

Wang, H. L. (2016). The roles of cohesive strength and toughness for crack growth in visco-elastic and creeping materials. *Engineering Fracture Mechanics*, **160**, 226-237.

Yang, T. G. (2022). Rate dependent fracture along a silicon/epoxy interface under mixed-mode loading conditions. *International Journal of Solids and Structures*, **257**, 111129.

Yusoff, N. I. (2011). Modelling the linear viscoelastic rheological properties of bituminous binders. *Construction and Building Materials*, **25**(5), 2171-2189.

Zehsaz, M. V.-T.-S. (2015). Modified creep constitutive equation for an epoxy-based adhesive with nonlinear viscoelastic behavior. *The Journal of Strain Analysis for Engineering Design*, **50**(1), 4-14.
# Appendix: Main constitutive equations - Beam on elastic foundation theory

In order to evaluate the peel stress and strain evolutions along the adhesive layer during a DCB test, the local static equilibrium the adherend can be written first, leading to the relations:

$$ \frac{dM}{dx} + T = 0 \qquad \frac{dT}{dx} - w\sigma = 0 \qquad (A1_1) $$

where $M$ and $T$ are the local bending moment and shear cohesive force along the adherend, respectively. $\sigma$ is the local peel stress in the adhesive, $w$ is the width of the adherend and $x$ is the position along the DCB specimen. The Timoshenko shear deformable beam model is then introduced to describe the adherend deflection. The corresponding constitutive equations are given by:

$$ M = EI \frac{d\varphi}{dx} \qquad T = \kappa GS \left[ \frac{dv}{dx} - \varphi \right] \qquad (A1_2) $$

where $v$ and $\varphi$ are the local adherend deflection and cross section rotation respectively. $E$ and $G$ are respectively the adherend Young's and shear moduli, $S$ and $I$ are respectively the area and second moment of inertia of the cross section and $\kappa=5/6$ is the shear correction factor. Combining equations (A1_1) and then (A1_2) we obtain two second order differential equations:
$$ \frac{d^2 M}{dx^2} + w\sigma = 0 \qquad \frac{d^2 v}{dx^2} - \frac{M}{EI} - \frac{w\sigma}{\kappa GS} = 0 \qquad (A1_3) $$

In the case of linear elastic adhesive bondline behavior, a simple linear relation is introduced between local peel stress and adherend deflection:

$$ \sigma = 2 \frac{E_a^*}{t_a} v = k_a v \qquad (A1_4) $$

where $E_a^*$ is the apparent Young's modulus of the adhesive, and $t_a$ is the thickness of the adhesive layer. Combining equations (A1_3) a characteristic equation is obtained which governs the peel stress and strain evolutions along the adhesive bondline:

$$ \frac{d^4 v}{dx^4} - \frac{w k_a}{\kappa GS} \frac{d^2 v}{dx^2} + \frac{w k_a}{EI} v = 0 \qquad (A1_5) $$

The characteristic roots of equation (A1_5) is equal to:

$$ \lambda_i = \pm \lambda_\sigma \sqrt{2 \left( \varepsilon^2 \pm i \sqrt{1 - \varepsilon^2} \right)} \qquad (A1_6) $$

with $\varepsilon = \lambda_\tau / 2\lambda_\sigma$ and:

$$ \lambda_\sigma = \frac{\sqrt{2}}{2} \left( \frac{wk_a}{EI} \right)^{1/4} \qquad \lambda_\tau = \sqrt{\frac{wk_a}{\kappa GS}} \qquad (A1_7) $$

$\lambda_\tau$ and $\lambda_\sigma$ controls the extension of the process zone ahead of the crack tip for respectively shear deformable and Euler Bernoulli beams only. For the DCB test, and considering a large
span length compared to the FPZ size, the general form of the peel stress distribution ahead of the crack tip is given by:

$$ \sigma = A\exp(\lambda_1 x) + B\exp(\lambda_2 x) \qquad (A1_8) $$

where $x = 0$ at the crack tip position, and $x < 0$ along the adhesive layer. $\lambda_1$ and $\lambda_2$ being the characteristic roots of equation (A1_5) with positive real part. Coefficients A and B are found by considering the boundary conditions corresponding to the specimen loading conditions and integrating equations (A1_1) to (A1_2). For DCB test we have:

$$ T(x = 0) = P \qquad M(x = 0) = aP \qquad (A1_9) $$

With, P, being the applied transverse force and, a, the instant crack length value.