Complete analytical solutions for double cantilever beam specimens with bi-linear quasi-brittle and brittle interfaces

In this work we develop a complete analytical solution for a double cantilever beam (DCB) where the arms are modelled as Timoshenko beams, and a bi-linear cohesive-zone model (CZM) is embedded at the interface. The solution is given for two types of DCB; one with prescribed rotations (with steady-state crack propagation) and one with prescribed displacement (where the crack propagation is not steady state). Because the CZM is bi-linear, the analytical solutions are given separately in three phases, namely (i) linear-elastic behaviour before crack propagation, (ii) damage growth before crack propagation and (iii) crack propagation. These solutions are then used to derive the solutions for the case when the interface is linear-elastic with brittle failure (i.e. no damage growth before crack propagation) and the case with infinitely stiff interface with brittle failure (corresponding to linear-elastic fracture mechanics (LEFM) solutions). If the DCB arms are shear-deformable, our solution correctly captures the fact that they will rotate at the crack tip and in front of it even if the interface is infinitely stiff. Expressions defining the distribution of contact tractions at the interface, as well as shear forces, bending moments and cross-sectional rotations of the arms, at and in front of the crack tip, are derived for a linear-elastic interface with brittle failure and in the LEFM limit. For a DCB with prescribed displacement in the LEFM limit we also derive a closed-form expression for the critical energy release rate, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G_c$$\end{document}Gc. This formula, compared to the so-called ‘standard beam theory’ formula based on the assumptions that the DCB arms are clamped at the crack tip (and also used in standards for determining fracture toughness in mode-I delamination), has an additional term which takes into account the rotation at the crack tip. Additionally, we provide all the mentioned analytical solutions for the case when the shear stiffness of the arms is infinitely high, which corresponds to Euler–Bernoulli beam theory. In the numerical examples we compare results for Euler–Beronulli and Timoshenko beam theory and analyse the influence of the CZM parameters.

of it even if the interface is infinitely stiff. Expressions defining the distribution of contact tractions at the interface, as well as shear forces, bending moments and cross-sectional rotations of the arms, at and in front of the crack tip, are derived for a linear-elastic interface with brittle failure and in the LEFM limit. For a DCB with prescribed displacement in the LEFM limit we also derive a closed-form expression for the critical energy release rate, G c . This formula, compared to the so-called 'standard beam theory' formula based on the assumptions that the DCB arms are clamped at the crack tip (and also used in standards for determining fracture toughness in mode-I delamination), has an additional term which takes into account the rotation at the crack tip. Additionally, we provide all the mentioned analytical solutions for the case when the shear stiffness of the arms is infinitely high, which corresponds to Euler-Bernoulli beam theory. In the numerical examples we compare results for Euler-Beronulli and Timoshenko beam theory and analyse the influence of the CZM parameters.

A
Cross-sectional area of a single DCB arm a Crack length a 0 Initial crack length b Width of a DCB C i Integration constants for the undamaged part (i = 1, . . . , 6) C j i Constants depending on C 3 , C 4 , ξ 3 and ξ 4 (i = 3, 4 and j = M, T ) C i Integration constants for the undamaged part that are zero because of the boundary conditions for a semi-infinite DCB (i = 1, . . . , 6) c 2 function cos(L cz κξ 2 ) c 2 function cos(L cz κξ 2 ) ch 1 function cosh(L cz κξ 1 ) ch 1 function cosh(L cz κξ 1 ) D Roots of the characteristic equation for the undamaged part (i = 1, . . . , 4) s 2 function sin(L cz κξ 2 ) s 2 function sin(L cz κξ 2 ) sh 1 function sinh(L cz κξ 1 ) sh 1 function sinh(L cz κξ 1 ) T Shear force in the upper DCB arm T 1 Shear force in the upper DCB arm on the part where the interface is undamaged T L 1 LEFM limit value of T 1 T 2 Shear force in the upper DCB arm on the part where the interface is damaged t i Roots of the characteristic equation for the damaged part (i = 1, . . . , 4) v Transversal displacement of the upper DCB arm v 1 Transversal displacement of the upper DCB arm on the part where the interface is undamaged v L 1 LEFM limit value of v 1 v 2 Transversal displacement of the upper DCB arm on the part where the interface is damaged x Co-ordinate along the interface x 0 Co-ordinate x corresponding to zero stress at the interface for the case of a linear-elastic interface with brittle failure x min Co-ordinate x corresponding to the minimum (maximum compressive) stress at the interface for the case of a linear-elastic interface with brittle failure α A constant defined as α = δ 0 /δ c β i Constants depending on ξ j , D kl , κ and a 0 (i = 1, . . . , 8, j = 1, 2, k = 1, . . . , 4 and l = 1, 2) γ Shear strain in the upper DCB arm Prescribed crack mouth opening displacement of a DCB a Part of the crack mouth opening displacement due to bending of the arms assuming that they are clamped at the crack tip δ CT Part of the crack mouth opening displacement due to the opening at the crack tip ϕ CT Part of the crack mouth opening displacement due to the rotation at the crack tip δ Relative displacement at the interface in mode I δ 0 Linear-elastic limit value of the relative displacement at the interface in mode I δ c Relative displacement at the interface corresponding to the total loss of interconnection in mode I ε r Relative error due to using Euler-Bernoulli instead of Timoshenko beam theory ε r Relative difference between CCM and LEFM solutions ζ i Constants depending on ω (i = 1, . . . , 4) ζ i Constants depending on ω (i = 1, . . . , 4) η A constant defined as η = λ/κ η E A constant defined as η E = λ E /κ θ Prescribed rotation on a DCB arm κ A constant depending on the bending stiffness of DCB arms and the softening part of the σ − δ traction-separation law of the interface λ A constant depending on the bending stiffness of DCB arms and the linear-elastic stiffness of the interface λ E A constant defined as λ E = √ 2λ/2 μ Shear modulus of DCB arms ν Poisson's ratio of DCB arms ξ i Constants depending on ψ (i = 1, 2) ρ i Constants defined as ρ i = ζ i /λ (i = 1, 2) Since introduced in 1960s by Dugdale (1960) and Barenblatt (1959), the use of cohesive-zone models (CZM) has become one of the most popular ways of describing fracture processes within the research community (Hillerborg et al. 1976; Alfano and Crisfield 2001;Park and Paulino 2011). Nowadays CZMs are widely implemented within the framework of finiteelement analysis (FEA) and interface elements, based on CZMs, can be found in element libraries of many commercial softwares for FEA used to solve delamination/debonding problems in 2D (Alfano and Crisfield 2001) and 3D (Park and Paulino 2011). However, because the accuracy of such computations is always dependent on the size of the FE mesh, they can be computationally demanding and suffer from convergence problems.
One of the ways to reduce the computational cost of the analysis was already proposed by the first and third author and consists of using beam finite elements instead of plane solids to model the bulk material of the specimens in 2D analysis of delamination. The results were presented for single-mode (I and II) and mixedmode delamination problems in geometrically linear (Škec et al. 2015) and non-linear analysis (Škec and Jelenić 2017). Compared to models which use plane solid FEs, the computational efficiency of the beam model was improved (the reduction of total number of degrees of freedom can go up to 40 %) without any significant loss in the accuracy. However, the beam model still suffered from convergence problems and spurious oscillations around the exact solution for cases of brittle interfaces when the mesh is not sufficiently refined. For this reason, it is always very useful when an analytical or semi-analytical solution can be found for cases of engineering interest.
In this paper we focus on mode-I delamination and the double cantilever beam (DCB) test, which is the standard test for determining fracture toughness in mode I. We use quasi-static and geometrically linear analysis where the arms of the DCB are modelled as Timoshenko beams and at the interface, a bi-linear CZM is used. Since the traction-separation law at the interface is composed of two linear parts and a final part with zero tractions, the solution of the problem can be obtained analytically, separately for each part of the interface whose state falls within one of the linear parts of the traction separation law. Because all the quantities of the problem will be expressed exactly with no need for discretisation, the computational cost is negligible and the obtained solutions have no spurious oscillations, which typically occur when delamination problems are solved using FE analysis (Alfano and Crisfield 2001;Blackman et al. 2003b;Škec et al. 2015). However, the idea of using analytical solutions for a DCB is not new and many researchers have contributed to the field in the last 50 years.
The simplest way to analytically model a DCB would be to assume that the arms of the specimen act as if they were cantilever beams clamped at the crack tip. Under this assumption, Benbow and Roesler (1957) used Euler-Bernoulli beam theory to establish the Griffith's energy balance for a flat-strip specimen where the crack gradually propagates down the middle by holding the specimen in a state of lengthwise compression. During 1960s and early 1970s, Ripling et al. (1971 introduced the DCB and tapered double cantilever beam (TDCB) tests and specimens. They also provided analytical formulae based on Irwin's energy approach (Irwin 1956) and Timoshenko beam theory (assuming that the arms are clamped at the crack tip) to compute the fracture toughness of the adhesive, which in 1974 became part of the American standard for determining fracture toughness of adhesive joints in mode I. The current version of that standard, ASTM D3433-99 (reapproved in 2012) (ASTM D3433-99 2012), still exclusively uses the same formulae proposed in Ripling et al. (1971). The formula for the DCB used in ASTM D3433-99 (2012) is also used in BS ISO 25217:2009(2009, where it is called 'simple beam theory' (SBT). We will adopt this terminology and use the term 'simple beam theory' (SBT) for formulations based on simple beam theories (Euler-Bernoulli or Timoshenko) and the assumption that the DCB arms act as if they were clamped at the crack tip.
However, Ripling et al. (1971) noticed that the SBT formula gives smaller deflections than those obtained from the experiments and they attributed it to not taking into account the rotations of the arms which take place at the crack tip. They also suggested that the results could be simply corrected by increasing the measured crack length by a fixed amount. Although they did not propose a method to obtain this crack length correction, this concept was further developed by other researchers (Blackman et al. 2003a) and became the basis for a data reduction scheme in BS ISO 25217:2009(2009 called 'corrected beam theory' (CBT) based on Euler-Bernoulli beam theory. de Moura et al. (2008) developed the so-called 'compliance-based beam method' (CBBM) where the corrected-crack-length concept was used with Timoshenko beam theory. Kanninen (1973) presented an analytical model for a DCB where the upper arm was modelled as a Euler-Bernoulli beam on elastic foundation (Winkler model), allowing for relative displacements and rotations at the crack tip and ahead of it. The very next year (Kan-ninen 1974) extended his formulation to account for the shear deformability of the arms and rotational stiffness of the interface, which was accomplished by using Timoshenko beam theory and Pasternak elastic foundation. However, Gehlen et al. (1979) (with Kanninen as the third author) showed that the rotational stiffness of the foundation springs is not relevant in a symmetrical DCB configuration. Kanninen's model (Kanninen 1974) was later extended by Williams (1989) to account for orthotropic material behaviour. Shahani and Forqani (2004), Shahani and Amini Fasakhodi (2010) developed solutions for a DCB model of finite length consisting of a Timoshenko beam lying on an elastic Winkler foundation for the conditions of fixed force and fixed displacement. Most of the mentioned papers (Kanninen 1973(Kanninen , 1974Gehlen et al. 1979;Shahani and Forqani 2004;Shahani and Amini Fasakhodi 2010) also investigate the dynamic analysis of unstable crack propagation and arrest in a DCB test. We will refer to the beam-on-elastic-foundation DCB models as 'enhanced beam theory' (EBT) models.
In EBT models the interface acts elastically linear up to a certain point where brittle failure occurs. However, fracture processes can usually introduce a certain level of quasi-brittle behaviour, which cannot be captured using EBT models. It is, however, well know that the quasi-brittle behaviour of the interface in a DCB test has an important influence on the structural response before the crack starts to propagate, whereas during the crack propagation the influence of the interface ductility is negligible. One way of introducing a quasi-brittle behaviour of the interface in the model is to use CZMs which account for progressive softening/damage after a certain critical value of the traction at the interface has been reached. Stigh (1988) developed an analytical solution for a DCB where the arms are modelled as Euler-Bernoulli beams and a bi-linear CZM is embedded at the interface. This solution was revisited and extended to account for a finite length of the specimen and a trapezoidal CZM by Dimitri et al. (2017). de Morais (2013) proposed a solution for a DCB with prescribed rotations (loaded with moments) where the arms are modelled as Timoshenko beams and a bi-linear CZM is embedded at the interface. Unlike Williams (1989) who gave the complete solution for the linear-elastic phase of the interface behaviour, de Morais (2013) took into account only real roots of the characteristic equation of the differential equation of the problem. We discuss this issue in detail in Sect. 2.2. We will refer to the models with quasi-brittle crack as 'cohesive crack models' (CCM) (Dimitri et al. 2017).
To the best of authors' knowledge, a complete analytical solution for a DCB with arms modelled as Timoshenko beams and the interface modelled using a bilinear CZM, which covers any of the two cases of prescribed rotations and displacement, is not available in the literature. Therefore, one aim of this paper is to fill this gap and to show a clear relationship between SBT, EBT and CCM solutions for a DCB for both Euler-Bernoulli and Timoshenko beam theory.
Furthermore, reducing the general CCM solution down to the SBT solution shows that even in the limit case of LEFM, rotations at the crack tip and in front of it still occur when Timoshenko beam theory is used. Thus, the assumption made in SBT that the arms act as if they were clamped at the crack tip is not applicable even for an infinitely stiff perfectly brittle interface, which is assumed in LEFM. This is somehow an expected result since even an infinitely stiff interface cannot prevent the bulk material of the arms to deform (rotate) around the interface. However, the fact that we can capture this behaviour using Timoshenko beam theory, to the best of authors' knowledge, has not been addressed in the literature so far. More generally, a comprehensive investigation of the analytical solutions for the LEFM limit is lacking in the literature. Thus, we call this novel approach 'enhanced simple beam theory' (ESBT). However, we will show that when the shear deformability of the arms is excluded from the model, which corresponds to Euler-Bernoulli beam theory, ESBT is equivalent to SBT. A second aim of this paper is to determine the interface stresses and the stress resultant profiles in the LEFM cases. We will also show that a novel LEFM-based formula for the determination of the critical energy release rate, G c , can be derived. This formula takes into account the rotation at the crack tip, unlike those currently available in the standards (ASTM D3433-99 2012;BS ISO 25217:20092009.
The outline of the paper is as follows. In Sect. 2, we define the problem and derive the general solutions of differential equations of the problem. In Sects. 3 and 4, we determine the integration constants for the cases of DCB with prescribed rotations and DCB with prescribed displacement, respectively. The solutions are given in a unified and compact form, thus avoiding cumbersome expressions which can be often encountered in such analytical solutions.
The presented general CCM solutions (based on Timoshenko beam theory and a bi-linear CZM at the interface) are then used to derive solutions for different particular cases. EBT solutions presented in Sect. 5 are obtained from CCM solutions by making the interface brittle (removing the softening branch in the traction-separation law of the interface), whereas ESBT (LEFM-based) solutions presented in Sect. 6 are obtained from EBT solutions by letting the interface stiffness go to infinity. The solutions for Euler-Bernoulli beam theory (including CCM, EBT and SBT solutions) given in "Appendix B" are obtained from the Timoshenko beam theory solutions by letting the shear stiffness become infinite.
Results obtained by means of the presented analytical solutions are studied for a number of cases, including sensitivity analyses on the interface strength. In particular, results for a linear interface with brittle behaviour and in the LEFM limit are presented and discussed in Sect. 7.3.

Problem description
Consider a double cantilever beam (DCB) of length L composed of two identical arms with depth h, width b and the initial crack length a 0 , as shown in Fig. 1. Each arm is modelled as a Timoshenko beam with a linearelastic constitutive law, where material properties are defined by Young's modulus, E, and shear modulus, μ. In a general case E and μ can have independent values, while for an isotropic material μ = 0.5E/(1 + ν), where ν is Poisson's ratio. Mode-I problem is considered, so that loads are applied symmetrically with respect to the mid-plane of the interface between the two arms. Therefore, stresses and strains in a DCB are symmetrical with respect to the mid-plane of the interface, which means that, for the sake of simplicity, only one arm can be considered in the analysis. In the present paper we will consider only the upper arm and assume that the x-axis is the centroidal axis of the arm (reference axis), while y and z axes are the principal axes of the arm's cross section (see Fig. 1).
In our approach, Timoshenko beam theory is used to model the DCB arms, which implies that we assume that displacements and rotations of the arms are relatively small compared to the specimen's dimensions.
The edges of the beams (where the interface is attached) in the general case can move both in x and y directions. However, for our problem, which is symmetric with respect to the mid-plane of the interface, there is no relative displacement at the interface responsible for mode-II delamination. This is because the upper and the bottom arm of a DCB at the same co-ordinate x have the same, but opposite cross-sectional rotation. Since Timoshenko beam theory is a geometrically linear theory, all points in a single cross-section of the arm experience the same displacement in the y-direction, i.e. v(x, y) = v(x). Therefore, for a DCB with symmetrical arm deformations, opening (mode-I) relative displacement at the interface, δ(x), corresponds to the sum of transverse displacements of both arms. This can be written as where v is the displacement of the upper arm. Note that in Fig. 1 the origin of the co-ordinate system is positioned at the crack tip, which means that at the left-hand end of the specimen x = −a 0 . The two types of tests we investigate in this paper are the DCB with prescribed rotations, θ , and the DCB with prescribed displacement, , as shown in Fig. 1. In the first case the crack propagates with a constant cohesive zone length (i.e. crack propagation is steadystate) (Suo et al. 1992;Škec et al. 2018), the cohesive zone being where softening/damage of the interface takes place. In the second case the crack propagation is not steady-state, but it tends to being so in the limit of infinitely long cracks, i.e. as a → ∞, where a denotes the crack length. Complete solutions for both cases are given in Sects. 3 and 4, respectively.
For the interface we use a bi-linear CZM consisting of a linear-elastic branch and a linear softening branch, followed by zero tractions for relative displacements greater than the critical value δ c , as shown on the right-hand side in Fig. 2. Here we emphasise that CCM solutions presented in this paper are not general solutions valid for any shape of the traction-separation law of the CZM, but are only valid when the interface behaviour can be assumed as bi-linear with progressive damage. Thus, for a bi-linear CZM law the solution will be obtained for three different phases in the crack propagation process, which is explained in detail in the following subsection.

Three-phase solution
In order to solve the problem of a DCB with an interface with embedded bi-linear CZM, the usual approach is to develop solutions for three different phases (Stigh 1988;de Morais 2013;Dimitri et al. 2017) which are namely: (i) linear-elastic behaviour, (ii) damage growth before crack propagation and (iii) crack propagation. The solution for each of these phases is derived in detail in the following sections according to Fig. 2, where only the upper half of the beam is shown. The origin of the co-ordinate system is in the left-most undamaged point for every phase, and it is moving to the righthand side as the damage at the interface increases and the crack propagates. Transversal displacements of the upper arm are denoted by v 1 (x) for x > 0 and by v 2 (x) otherwise for reasons which will be explained in Sect. 2.2.

Linear-elastic behaviour
In this phase the entire DCB acts as a linear-elastic body, which means that besides the DCB arms (with a linear-elastic constitutive law), the interface behaviour is linear-elastic, too. For all points at the interface undergoing separation (including points A and B in Fig. 2a), δ(x) < δ 0 and σ (x) < σ max , where δ 0 and σ max are linear-elastic limit values of the separation δ and the traction σ on the interface, respectively. This also means that, as long as δ(x) < δ 0 and σ (x) < σ max , no energy dissipation (damage) at the interface occurs and the initial (undeformed and undamaged) configuration can be recovered if the load is removed. Once at the crack tip δ(0) = δ 0 , we enter a new phase where damage at the interface starts to develop. This phase is explained in the following section.
It is well known that a zone of compressive stresses (σ (x) < 0, resulting in δ(x) < 0) ahead of the crack tip exists in a DCB (Alfano and Crisfield 2001;de Morais 2015;Dimitri et al. 2017). However, our CZM assumes that in compression no softening (or damage) occurs, and the behaviour is still linear elastic. In Fig. 2a we can see that point C experiences negative (compressive) contact tractions and negative relative displacements at the interface. For a zero-thickness adhesive layer this results in a non-physical overlapping of the arms, which in our model is allowed and does not create any additional stresses. But because the interface thickness in our model does not influence the results (although it is well know that in reality a different thickness of the same adhesive results in a different structural behaviour) and we can assume that the interface is sufficiently thick to prevent the arms coming in direct contact.

Damage growth before crack propagation
In this phase, the interface can be divided in two zones. In the first zone, just ahead of the crack tip, damage is developing (although the crack is still not propagating) and δ(x) ∈ [δ 0 , δ c ), where x ∈ (−L cz , 0]. In the second zone interface behaviour is again linear-elastic, i.e. δ(x) < δ 0 , where x > 0. Note that according to Fig. 2b, x = −L cz corresponds to the initial crack tip and x = 0 corresponds to the point where the first and the second zone meet. In this phase, the damage builds up with loading and L cz increases from 0 (corresponding to δ(−L cz ) = δ(0) = δ 0 ) to a certain limit value (corresponding to δ(−L cz ) = δ c ) at which the crack begins to propagate. Note that, according to Fig. 2b, the σ − δ relationship in the first zone is also linear, but with softening (damage).

Crack propagation
This phase is similar to the previous phase (Sect. 2.1.2) with the difference that here the relative displacement at the current crack tip, δ(−L cz ) = δ c , remains unchanged during the whole phase. For a DCB with prescribed rotations the cohesive zone length remains constant for any position of the crack (steady-state crack propagation), while for a DCB with prescribed displacement the cohesive zone length will change (non-steady-state crack propagation). Thus, when the crack propagation is steady state the deformed shape of the interface shown in Fig. 2c, and the contact traction distribution over the interface, σ (x), simply translate to the right-hand side. This is why point A in Fig.  2c is the point at the interface where the crack tip is currently located. When the crack propagation is not steady-state, we still have δ(−L cz ) = δ c and δ(0) = δ 0 for any position of the crack tip, but the deformed shape and the contact traction distribution at the interface (including L cz ) change as the crack propagates. The part of the interface which has been completely damaged (σ (x) = 0) is excluded from the domain of the solution for v 2 (x) and becomes a part of the DCB arm.
2.2 Solutions of differential equations of the problem

The differential equation of the Timoshenko beam reads
where E I is the bending stiffness (with second moment of area I = b h 3 /12), μA s is the shear stiffness (with corrected shear area A s = b h k s , where k s is the shear correction coefficient) and q(x) is distributed transverse load along the beam axis, assumed to be positive when pointing upwards. Furthermore, in accordance with the co-ordinate system from Fig. 1, we have where T (x), γ (x) and ϕ(x) are the shear force, shear strains and cross-sectional rotation, respectively. We assume that the self-weight of the DCB arms is negligible compared to the magnitude of the external forces or moments acting on the DCB. Therefore, only contact tractions will act on the arms as a distributed load, and thus where the negative sign appears because positive (tensile) contact tractions at the upper arm are pointed downwards. The part of the DCB arms separated by the initial crack, or where the interface has been completely damaged, is excluded from the domain of the solution v(x). However, the moment or the force applied at the point of the prescribed rotation or displacement, respectively, (which is outside of the mentioned domain) are taken into account via the boundary conditions at the crack tip. This will be explained in detail in the following sections. We define σ (x) according to the traction-separation law of the CZM. Thus, in our case we will define σ (x) separately for the linear-elastic (δ(x) ≤ δ 0 ) and for the linear softening part (δ 0 < δ(x) ≤ δ c ) as Substituting (6) in (5) and then in (2) and taking into account (1) we obtain two differential equations where Equation (7) is used in all phases on the undamaged part of the interface (x ≥ 0), while Eq. (8) is used only in phases 2 and 3 on the damaged part of the interface (x ∈ [−L cz , 0)). We will denote the solutions of Eqs. (7) and (8) by v 1 (x) and v 2 (x), respectively, and derive them in the following sections.

Solution on the undamaged part of the interface
Assuming the solution of Eq. (7) in a form v 1 (x) = e r x , where r is a constant, results in a characteristic equation with four roots, namely where Since ω ≥ 0, we will have all real roots for ω > 1, all complex roots for ω < 1 and multiple real roots for ω = 1. For each of these cases, we can now give: 1. The solution of Eq. (7) for ω > 1: where C 1 , C 2 , C 1 and C 2 are integration constants.
2. The solution of Eq. (7) for ω < 1: where and C 3 , C 4 , C 3 and C 4 are integration constants. Note that because of we have 3. The solution of Eq. (7) for ω = 1: where C 5 , C 6 , C 5 and C 6 are integration constants.
In our approach we will assume that during crack propagation the crack tip is always sufficiently distant from the right-hand end of the DCB, in a way that any boundary conditions at the right-hand end (free or clamped) do not influence the results. This is equivalent to assuming a semi-infinite DCB. Therefore the domain of the undamaged part of the interface will have boundary conditions at x = 0 and x = ∞, where we can write: where ϕ 1 (x) is the cross-sectional rotation on the undamaged part of the interface. According to Timoshenko beam theory is the cross-sectional shear force on the undamaged part of the interface. Because T 1 (∞) = 0, applying boundary conditions (19) to solutions (13), (14) and (18) gives C i = 0, where i = 1, . . . , 6. Thus, the solution of Eq. (7) for a semi-infinite DCB can be written in a general form as where x ≥ 0.
Remark 2.1 For the sake of simplicity and because of the extreme unlikelihood that the value ω = 1 occurs in a real case, the solutions in the following sections are given only for the cases when ω > 1 and ω < 1. However, the results for ω = 1 are given separately in "Appendix A" for completeness. In the numerical examples presented in Sect. 7 we will show that, unlike stated in de Morais (2015), all solutions from Eq. (20) are possible for realistic values of geometrical and material properties of a DCB.

Solution on the damaged part of the interface
Assuming the solution of Eq.
where t is a constant, results in a characteristic equation with four roots, namely where Since ξ 1 and ξ 2 are real for any value of ψ, roots t 1 and t 2 are always real, whereas t 3 and t 4 are always imaginary. Thus, the solution of Eq. (8) can be written as In the following sections, the problem is solved and the integration constants are determined for each phase for a DCB with either prescribed rotations or prescribed displacement.
opposite rotation θ is prescribed, causing the opening of the DCB along the interface due to equal, but opposite concentrated moments, M acting at the point of the prescribed rotation. This implies that the disconnected parts of the DCB arms are in pure bending and the shear force at the crack tip is zero during all phases. Each of the three solution phases is derived in detail in the following sections.

Linear-elastic phase
Using (4), (5) and (6), we can express the crosssectional bending moment, M 1 (x) and shear force, T 1 (x), in the upper DCB arm as: where index 1 refers to the solution of Eq. (7), i.e. linear elastic behaviour. Using solutions (20), Eqs. (24) and (25) can be written as where Please recall that the solutions for the case when ω = 1 are given in "Appendix A". Using expressions (26) and (27), from the following boundary conditions at the crack tip: we obtain for ω > 1, and for ω < 1. The transition from the linear-elastic phase to the phase of damage growth before crack propagation, is defined by the condition for both ω > 1 and ω < 1 one has Hence, the moment M L , leading to the transition from the first phase (linear-elastic behaviour) to the second phase (damage growth before crack propagation) reads for any value of ω (including ω = 1, as reported in "Appendix A").
The relative displacement of the DCB arms at the point of application of the moment (crack mouth opening displacement) is computed as where are the crack mouth opening displacements due to bending of the arms behind the crack tip, rotation of the crack tip and opening at the crack tip, respectively. Note that displacement is the total opening of the crack mouth, and it takes into account both arms, which is why factor 2 is used in expressions (37)-(39). Equation (36) can be now rewritten as where v 1 (0) is defined in (34) and v 1 (0) can be obtained from (20) as in which the final expression is valid for any value of ω. It follows that is obviously linear for any value of ω (including ω = 1, as reported in "Appendix A"). We will use Eq. (36) as a general solution for the crack mouth opening displacement of a DCB, which is valid for all three solutions phases, not only for a DCB with prescribed rotations. However, in general a , ϕ CT and δ CT are computed differently each time.
and σ (x) for this phase are given in "Appendix E.1".

Phase of damage growth before crack propagation
In this phase, a cohesive (or damage-process) zone is developing in front of the crack tip. As already mentioned, for the cohesive zone (x ∈ [−L cz , 0]) we use the solution (23), whereas for the zone of linear-elastic behaviour (x ≥ 0) we use the solution (20).
There are six constants to determine, two for the undamaged part (C 1 and C 2 if ω > 1, or C 3 and C 4 if ω < 1) and four for the damaged part (D 1 , . . . , D 4 ) in order to obtain the complete solution. First, from we obtain We then impose the continuity conditions at the origin of the co-ordinate system: where ϕ 2 (x) is the cross-sectional rotation on the damaged part of the interface, M 1 (x) and T 1 (x) are defined according to (26) and (27), respectively, and are the cross-sectional bending moment and the shear force on the damaged part of the interface, respectively. Using the solution (23) these expressions can be written as where the property ξ 1 ξ 2 = 1 is used to simplify the expressions. Note also that, because (45) we can then express constants D i (i = 1, . . . , 4), in terms of C 2 for ω > 1 or C 3 for ω < 1 in the following general form where j = 2 for ω > 1 and j = 3 for ω < 1, D i1 and D i2 are constants depending on the value of ω, as explained below, and D 0 = 2(ξ 2 1 + ξ 2 2 ) ≡ 4(ψ 2 + 1).
As indicated in (50) and explained later, D i and C j are not true constants, but parameters depending on L cz . If we set the values of the constants D i1 and D i2 (i = 1, . . . , 4) are the following: with η = λ/κ. Parameter C j (L cz ) is determined from the last (sixth) remaining boundary condition as where j = 2, 3 and To summarise, parameters D i (L cz ) needed to define displacement field (23) within the damaged zone are determined from (50) for any value of ω, where C j (L cz ) is computed from (55). For ω > 1, C 2 (L cz ) is defined according to (55) ( j = 2) and C 1 (L cz ) is obtained from (44) (55) ( j = 3) and, according to (44), C 4 = δ 0 /2 for any value of L cz .
The value of the applied moment M is increasing during this phase and L cz is strictly increasing with increasing M. Hence, M will also depend on L cz . This function, M(L cz ), can be obtained from the following condition with M 2 (−L cz ) computed from (48). Obviously, L cz and M(L cz ) can increase only up to a certain limit after which the crack begins to propagate and we enter the third phase of the solution. The crack starts to propagate as soon as the relative opening at the crack tip reaches the critical value δ c or This condition represents a highly non-linear equation in terms of L cz , which contains trigonometric and hyperbolic functions, and cannot be solved in a closed form. Thus, a numerical solver is needed and in the present work we use a simple Newton-Raphson iterative procedure (more information regarding the numerical solver is given in Sect. 7). We will denote the solution of (58) as L cz and the maximum applied moment by M max = M(L cz ). According to (36), where a can be still defined according to (37), but now the crack mouth opening can be computed as where v 2 (−L cz ) and v 2 (−L cz ) are evaluated from (23) at the co-ordinate x = −L cz .

Crack propagation phase
In the previous phase (damage growth before crack propagation) the applied moment is increasing from M L to M(L cz ) = M max . Since for a DCB with prescribed rotations there is no shear force at the crack tip during all phases, it means that only the applied moment M (same at the crack tip as at the point of application) is responsible for crack propagation. Obviously, the crack will propagate when the applied moment reaches the value M max and this value will not change as the crack propagates. A constant value of M max during crack propagation implies that the boundary conditions at the crack tip remain constant and that L cz = L cz , during crack propagation. This kind of behaviour is known as 'steady-state crack propagation'. Thus, unlike in the previous phase, in this phase L cz is not a variable.
The interface is again divided in two domains, the undamaged one (x ≥ 0), where function v 1 (x) is defined according to (20), and the damaged one (x ∈ [−L cz , 0]), where function v 2 (x) is defined according to (23). Continuity conditions (45) still apply and constants D i are now obtained as which is similar, but not equivalent to (50), because D i and C j are now true constants and not functions of L cz . On the other hand, constants D 0 , D i1 and D i2 are still defined according to (51) and (53). Constant C j is now determined from the condition where j = 2, 3 and To summarise, constants D i (i = 1, . . . , 4) are determined from (62) for any value of ω. For ω > 1, (64) and where C j is defined in (64) and C 4 = δ 0 /2, according to (44).
Since in the crack propagation phase the applied moment remains constant, it can be computed from (57) as M max = M(L cz ). Crack mouth opening can be then computed as where, theoretically, the value of a can go to infinity. Note also that v 2 (−L cz ) and v 2 (−L cz ) are constants and thus the function (a) in this phase is quadratic. Because in this phase M does not change, does not depend on M and we cannot define (M).

DCB with prescribed displacement
In this section we consider a DCB with prescribed displacement where, according to Fig. 1b, at the left-hand end the bottom arm is pinned, whereas the upper arm is pulled upwards. In order to prescribe a displacement at the left-hand side of the upper arm, a vertical force F must be applied at the same place and in the same direction. Thus, unlike in the case of a DCB with prescribed rotations, at the cracked portion of a DCB with prescribed displacement there is bending and shear in the arms, which will make the problem slightly more complex. Furthermore, because the crack propagation in the case of a DCB with prescribed displacement is not steady-state (L cz changes during crack propagation), the solution for the third phase (crack propagation) will be also more complex, compared to the case of DCB with a prescribed rotations where L cz = L cz is constant during crack propagation. Each phase of the solution is explained in detail in following sections.

Linear-elastic phase
In this phase the boundary conditions at the crack tip read from which, using (26) and (27), constants for ω > 1 and for ω < 1 are determined. The limit value of the force, F L , at which damage starts to develop at the interface is obtained from condition (33), which gives so that Note that a real value of F L is obtained for any value of ω (including ω = 1, as reported in "Appendix A"). The same applies to v 1 (0) in (70). The prescribed displacement is in this case equal to the crack mouth opening , which we define according to (36), where now we have The crack mouth opening displacement can be then written as where v 1 (0) is defined in (70) and v 1 (0) = − F E I λ 2 2 ω + 1 + a 0 λ 2(1 + ω) . (76) Note that although the terms responsible for shear deformations in (72) and (73) cancel out in (75), shear deformability is taken into account through ω in (70) and (76). Equation (76) is valid for any value of ω (including ω = 1, as reported in "Appendix A"). Equation (75) can be finally written as The final forms of functions v 1 (x), v 1 (x), ϕ 1 (x), M 1 (x), T 1 (x) and σ (x) for this phase are given in "Appendix E.2".
4.2 Phase of damage growth before crack propagation In this phase, we divide the interface into the undamaged domain (x ≥ 0) and the damaged domain (x ∈ [−L cz , 0]), and continuity conditions (45) between the two still apply. Using these conditions and condition (43) we can again define parameters D i (L cz ), where i = 1, . . . , 4, using (50) and constants C 1 and C 4 using (44). However, solution (55) for parameters C j (L cz ), where j = 2, 3, is no longer valid because the boundary condition (54) in the case of prescribed displacement becomes which cannot give us the solution for C j (L cz ) because the function F(L cz ) is yet unknown. An additional boundary condition gives where M(L cz ) is defined in (57). Now, using (80), from (78) it follows that where j = 2, 3 and Let us recall that, because of definitions (52) and (53), solution (81) for C j (L cz ) automatically returns the value of C 2 (L cz ) for ω > 1 and C 3 (L cz ) for ω < 1. Thus, in the case when ω > 1 we compute C 2 (L cz ) from (81) and then C 1 (L cz ) follows from (44) as C 1 (L cz ) = δ 0 /2 − C 2 (L cz ). For ω < 1 we compute (81), whereas C 4 is defined in (44) as C 4 = δ 0 /2. Using solution (81), from (50) we can then obtain parameters D i (L cz ), where i = 1, . . . , 4, which are needed to compute M(L cz ) (see (57)) and finally F(L cz ) according to (80). As explained in Sect. 3.2, in this phase L cz grows from 0 to a value corresponding to the initiation of crack propagation, i.e. transition to the third phase. This is a maximum value for L cz , because during crack propagation L cz decreases and asymptotically tends to a minimum value when a → ∞, as is discussed in the next section. Therefore, this initial maximum value of L cz at the initiation of crack propagation will be denoted by L max cz . In order to obtain this value, the same approach as for a DCB with prescribed rotations is followed, i.e. condition (58) is imposed. This is again a highly nonlinear equation in terms of L cz , which is solved numerically (in our approach Newton-Raphson procedure is used).
The prescribed displacement in the second phase can be computed as where v 2 (−L cz ) and v 2 (−L cz ) are evaluated using (23) at the co-ordinate x = −L cz .

Crack propagation phase
As previously mentioned, in the case of a DCB with prescribed displacement, the cohesive zone length decreases during crack propagation from L max cz asymptotically approaching a lower limit value L min cz . This means that the crack propagation is not steady state, but it approaches steady state for infinitely long cracks (Dimitri et al. 2017). Because L min cz corresponds to a steady-state crack propagation, it must have the same value as L cz found in the identical DCB loaded with prescribed rotations, where crack propagation is always steady state, i.e. L min cz = L cz . Note that both L max cz and L min cz are obtained by numerically solving Eq. (58), where constants D i (i = 1, . . . , 4) in (23) in the former case are computed using (81) in (50), whereas in the latter case they are computed using (55) in (50).
Continuity conditions (45) and condition (43) are still valid in the third phase, which gives us solution (50). From the condition we obtain where j = 2, 3. Note that this expression is different from (64) because here L cz is a variable. Function F(L cz ) is obtained from conditions (78) and (49). Note that a is also a function of L cz , i.e. a = a(L cz ), and can be determined from the condition which gives where M(L cz ) and F(L cz ) are defined according (57) and (78), respectively. Prescribed displacement can be now expressed as a function of L cz as where v 2 (−L cz ) and v 2 (−L cz ) are evaluated using (23) at the co-ordinate x = −L cz .
Remark 4.1 Solutions developed in Sects. 3 and 4 (based on Timoshenko beam theory) represent general solutions from which we can easily derive three other particular cases, for DCBs with either rotations or displacement prescribed: 1. Solutions for Euler-Bernoulli beam theory. These solutions are obtained by simply letting the shear modulus μ → ∞ and are presented in "Appendix B". 2. Solutions for a linear-elastic interface with brittle failure. These solutions, presented in Sect. 5 for Timoshenko beam theory and in "Appendix C" for Euler-Bernoulli beam theory, are obtained by letting δ c → δ 0 , which means that there is no damage before crack propagation, i.e. only the first and the third solution phases remain. 3. LEFM solutions. These solutions, given both for Timoshenko beam theory in Sect. 6 and Euler-Bernoulli beam theory in "Appendix C", are obtained from the solutions for a linear-elastic CZM with brittle failure by letting δ 0 → 0.

Analytical solutions for a DCB with linear-elastic interface with brittle failure (EBT)
Here we assume that the interface behaves as linearelastic up to a certain point after which brittle failure (instantaneous loss of cohesion) occurs. In our model, this is achieved by setting δ 0 = δ c , which removes the softening branch in σ −δ diagram shown in Fig. 2. This also means that the second part of the solution (damage growth before crack propagation) does not exist and that L cz = 0 during crack propagation. Referring to Sect. 2.2, we can now rewrite Eq. (6) as Solutions (20) for v 1 (x), as well as definitions (9), (12) and (15) for λ, ω and ζ i (i = 1, . . . , 4) are still valid. Solution (23) for v 2 (x) is no longer applicable because L cz = 0. We will now derive the solution for a DCB with a brittle interface, first for a DCB with prescribed rotations and then for a DCB with prescribed displacement.

DCB with prescribed rotations
The first (linear-elastic) phase of the solution presented in Sect. 3.1 applies completely to the case of linearelastic interface with brittle failure. Thus, no modifications in the expressions presented in Sect. 3.1 are needed. Moreover, this phase can be now called 'linearelastic behaviour before crack propagation.' The crack mouth opening displacement in the phase of crack propagation now reads where M max is the value of the applied moment when the crack starts to propagate corresponding to M L defined in (35).
Since for a linear-elastic interface with brittle failure the critical energy release rate, G c , can be written as we can rewrite Eq. (9) 1 as and from (35) obtain which is equivalent to the well-known formula used to compute G c (or the critical value of the J integral, J c ) for a DCB with prescribed rotations (Rice 1968;Freiman et al. 1973;Suo et al. 1992;Sørensen et al. 1996). Note that M max and G c are independent of the beam theory used, i.e. the shear deformability of the arms does not influence their values.

DCB with prescribed displacement
For a DCB with prescribed displacement the first part of the solution presented in Sect. 4.1 is entirely valid in the case when the interface is linear-elastic with brittle failure. The peak load at the point when the crack starts to propagate, F max , corresponds to F L defined in (71), from where by substituting (9) and (92) it follows that For the crack propagation phase we have where Remark 5.1 Note that presented solutions in terms of the applied load and the crack mouth opening displacement of a DCB with either prescribed rotations or prescribed displacement for the case of a linear-elastic interface with a brittle failure (EBT) are valid for all values of ω (ω < 1, ω = 1 and ω > 1).

LEFM solutions for a DCB with prescribed rotations and a DCB with prescribed displacement: enhanced simple beam theory (ESBT)
The presented model for a linear-elastic interface with a brittle crack still allows some opening at the interface in the linear-elastic range before crack starts to propagate (δ(x) < δ 0 ). If this initial linear-elastic behaviour is excluded from the model by letting δ 0 → 0, while keeping G c constant (which also results in σ max → ∞), we obtain solutions equivalent to those given by linear-elastic fracture mechanics (LEFM). However, we will not a-priori assume that the DCB arms act as if they were clamped at the crack tip, which is usually done in the SBT approach. In this way we will show that in the limit case of LEFM the arms rotate at the crack tip and in front of it (even though their centre-lines there remain straight), which means that the clamped conditions at the crack tip cannot be obtained even for an infinitely stiff perfectly brittle interface. This is due exclusively to the shear deformability of the arms, which is accounted for in Timoshenko beam theory. For Euler-Bernoulli beam theory, the limit case of LEFM indeed corresponds to SBT and the arms do not rotate at the crack tip and in front of it. Thus, for our LEFM-limit solution for Timoshenko beam theory we will adopt the term 'enhanced simple beam theory' (ESBT). In the following subsections we will present only the final results for a DCB with prescribed rotations and a DCB with prescribed displacement, respectively, while the complete derivation is given in "Appendix E".

DCB with prescribed rotations
As shown in "Appendix E", in the limit case of LEFM for any where D(x) is the Dirac distribution centred at zero (which means that, rigorously speaking, the stress is not a 'proper' function but a generalised one), and H(x) is the Heaviside function defined by: From the above expressions we can see that, although there are no relative displacements at the crack tip and in front of it (v 1 (x) = 0 for any x ≥ 0), rotations at the crack tip and in front of it are still allowed to occur when Timoshenko beam theory is used. On the one hand, this result is expected, because it is intuitive that the independence of rotation and deflection in Timoshenko beam theory allows this theory to capture the deformation in front of the crack tip, which occurs also in the LEFM limit. This is in contrast with Euler-Bernoulli theory (see "Appendix E.3"), for which the absence of displacements also means a zero rotation and, ultimately, no deformation in front of the crack tip, but also with the widely used assumption, made in the SBT, that the arms of a DCB act as if they were clamped at the crack tip. In other words, to the best of the authors' knowledge, the ability of Timoshenko's beam theory to capture the crack tip rotation also in the LEFM limit has not been explored so far, although something similar was done in EBT for a linear elastic interface with finite stiffness and brittle failure (Kanninen 1973;Williams 1989). This is why we call our approach ESBT.
Moreover, because the DCB arms deform (rotate) in front of crack tip, contact tractions at the interface, as well as the shear forces and bending moments in the arms, appear in front of the crack tip, with an exponential decay as x → ∞. We can notice that at the crack tip there is a jump in the shear force (from 0 to −M √ μA s /E I ) which corresponds to a transition in the bending moment diagram from the constant value M in the cracked portion of the arms to the function (103) in front of the crack tip. This implies that in the limit case of LEFM there is a concentrated transversal cohesive force exchanged at the crack tip, so that the interface stress is the sum of a compressive smooth part and the Dirac distribution centred at zero. Because at the crack tip the cross-sectional rotations of the arms must be continuous and there is a jump in the shear force in the arms, the function v 1 (x) is also discontinuous due to ϕ 1 (x) = v 1 (x) + T 1 (x)/μA s . Expressions (98)-(103) are valid only for the phase of linear-elastic behaviour before crack propagation. However, analogous expressions for the phase of crack propagation can be obtained simply by substituting M with M max . Finally, according to (40), we can express the crack mouth opening displacement as where M max is defined in (93) and the second term in the parentheses in both expressions represents the rotation of the arms at the crack tip.

DCB with prescribed displacement
As shown in "Appendix E", in the limit case of LEFM for any where again we see that even in the limit case of LEFM the arms rotate at and in front of the crack tip when Timoshenko beam theory is used. The discussion regarding Eqs. (98)-(103) also applies here, with the only difference that for a DCB with prescribed displacement in the cracked portion of the arms we have a constant shear force and a linear distribution of bending moments. Note that expressions (107)-(112) are valid only for the phase of linear-elastic behaviour before crack propagation. However, analogous expressions for the crack propagation phase can be obtained by substituting F with F(a) (defined in (116)) and a 0 with a.
The crack mouth opening displacement before crack propagation follows from (75) as where from (95) we have During crack propagation the crack mouth opening displacement is given by where from (97) we have Note that in (113) and (115) the term outside the parentheses represents the arm deflection due to bending according to Euler-Bernoulli beam theory, the second term in the parentheses is due to shear deformability of the arm, while the third term in the parentheses represents the influence of the rotation of the arms at the crack tip. From (116) we can obtain the critical energy release rate of a Timoshenko DCB with prescribed displacement as where, for the sake of simplicity, we denote F(a) simply by F. Usually in the SBT solution (Ripling et al. 1971;ASTM D3433-99 2012;BS ISO 25217:20092009, only the first two terms in expression (117) are taken into account because it is assumed that the DCB arms are clamped at the crack tip. The third therm in (117), which takes into account the rotation at the crack tip, to the best of authors' knowledge, has not been recognised so far for the limit case of LEFM. Equation (117) represents the ESBT expression for G c .
Remark 6.1 In "Appendix E.3" we show that if Euler-Bernoulli beam theory is used, there are no crosssectional rotations, shear forces or bending moments of the arms in front of the crack tip. However, singularity of contact tractions at the interface, as well as discontinuity of the shear stresses and bending moments in the arms, take place at the crack tip. These conditions are equivalent to clamping DCB arms at the crack tip and explain why the formulae obtained for the limit case of LEFM in "Appendix D" indeed correspond to widely used formulae in SBT.

Numerical examples
In this section, for a DCB with a bi-linear CZM at the interface, the analytical solutions derived in this paper using Timoshenko beam theory will be compared to the numerical results obtained with an equivalent finite-element (FE) model, in which the same beam theory and CZM are used (Škec et al. 2015), and to the Euler-Bernoulli beam theory analytical solutions derived in "Appendix B". The latter will allow us to investigate the influence of shear deformability of DCB arms on the results. LEFM solutions obtained in Sect. 6 will also be presented as limit cases for a brittle interface. Referring to Fig. 1, we consider a DCB with dimensions h = 6 mm, b = 25 mm and a 0 = 30 mm. In the present analytical solution, the length of the specimen is assumed to be infinite (L = ∞), i.e. the crack is always sufficiently distant from the right-hand (non-loaded) end of the DCB. Material data for the DCB arms and the interface used in numerical examples is presented in Table 1, where it can be noted that the maximum contact traction σ max is varied between 7.5 and 120 MPa, while keeping the area under the traction-separation law, , constant. This gives us 5 cases of different brittleness of the interface, where σ max = 7.5 MPa represents an extremely ductile case and σ max = 120 MPa an extremely brittle case. Ratio α = δ 0 /δ c is kept constant in the first set of examples (α = 0.01 for all cases), whereas in the last example is varied. Values of E and ν for the DCB arms presented in Table 1 correspond to aluminum, the shear modulus μ is calculated as for an isotropic material, i.e. μ = 0.5E/(1 + ν), and for the rectangular crosssection considered, k s = 5/6. The numerical model used is the multi-layer beam model presented in Škec et al. (2015) where we assume a total length of the specimen L = 200 mm. A total number of 2000 2-node Timoshenko beam elements are distributed evenly over the upper half of the DCB, meaning that the element length is 0.1 mm. Such a fine mesh is used to eliminate or at least minimise the influence of discretisation-caused spurious oscillations on the results (Alfano and Crisfield 2001;Škec et al. 2015). A 4-node interface element is attached to every beam FE from x = a 0 to x = L making a total of 1700 interface elements. The solution is obtained using displacement control and Newton-Raphson iterative procedure. Because our numerical model has 4002 degrees of freedom (one transverse displacement and one crosssectional rotation per node), in each iteration of each increment, 4002 linear equations are solved in order to obtain the cross-head displacement. In our analytical solution, we obtain the cross-head displacement from a single closed-form solution. The same applies to any other quantity we want to obtain. Furthermore, all analytical solutions, unlike the numerical ones, are perfectly smooth.
In the following sections we will present the results for the DCB first with prescribed rotations and then with prescribed displacement.  Fig. 3, the reaction moment, M, is plotted against the crack mouth opening displacement, . From Fig. 3b, which is the zoom of Fig. 3a, we can first observe that the results for Timoshenko beam theory from the analytical and the FE model perfectly match, as expected.
We can also note that the results approach the LEFM solution, which is obtained from our model by letting δ 0 → 0 and σ max → ∞ (as shown in Sect. 5.1), as σ max increases. The same behaviour can be observed when Euler-Bernoulli beam theory is used. Differences between Timoshenko and Euler-Bernoulli beam theories are more pronounced for higher values of σ max and especially for the LEFM solution. From Fig. 4a it can be clearly seen that the crack will start to propagate sooner (i.e. for smaller crack mouth opening displacements) when the interface is more brittle. The numerical model again agrees perfectly with the analytical solution, and there are some differences between Euler-Bernoulli and Timoshenko beam theory, which again become more significant as σ max increases.
In Fig. 4b we can finally compare the cohesive zone lengths, L cz , for Timoshenko and Euler-Bernoulli beam theories. As expected, L cz , which is highly influenced by the value of σ max , remains constant during crack propagation. Differences between Timoshenko and Euler-Bernoulli beam theory solutions are now even more pronounced, especially for more brittle cases. Again, the numerical results match perfectly with those obtained from the analytical solution for Timoshenko beam theory.

DCB with prescribed displacement
In this section, the same geometrical and material data as in Sect. 7.1 is used for the case of a DCB with prescribed displacement. We present a comparison of the solution from Sect. 4 for Timoshenko beam theory with the analogous solution for Euler-Bernoulli beam theory (presented in "Appendix B.2") and the LEFM solu- tion. This is done by comparing the values of the applied force, F, computed for each model at the same value of the crack mouth opening displacement, . Recall that in Sects. 4.2 and 4.3, in the solution phases which include damage, we have defined all quantities as functions of L cz . Therefore, once the limit values L max cz and L min cz have been obtained numerically, a series of suitable values of L cz can be chosen for the second and the third phase and the solution of the problem is obtained analytically in a closed form. On the other hand, for the purpose of the present rigorous numerical comparison, it is necessary to express these quantities in terms of . However, deriving functions F( ), a( ) and L cz ( ) from expressions given Sects. 4.2 and 4.3 is impossible, because the functional dependence on L cz is highly non-linear. Thus, values of F, a and L cz for a certain value of are determined numerically using Newton-Raphson procedure. This computational effort is still negligible in comparison with the standard FEA because we only need to solve a single non-linear equation for each increment of .
In Fig. 5, the reaction force, F, is plotted against the prescribed crack mouth opening displacement, , for different values of σ max . Both Timoshenko and Euler-Bernoulli beam theories are used. We can see that σ max has a considerable influence on the results in the first two phases (before crack starts to propagate), especially in the second phase, when the stiffness of the DCB progressively decreases as damage is developing in front of the crack tip. This also results in a reduction of the peak force. In Fig. 5b it can again be seen that differences between Timoshenko and Euler-Bernoulli beam theory become more pronounced as σ max increases, i.e. the interface becomes more brittle. Results from the FE model based on the Timoshenko beam theory perfectly match the results obtained from the analytical solution. In the third phase (crack propagation) all curves are extremely close, but they do not coincide perfectly. This is better shown in Fig. 6a, where it can be noted that the differences of the results in the crack propagation phase indeed exist, but they are too small to be appreciated on a normal scale.
In fact, a simple argument explains why the curve corresponding to a finite value of σ max must lie above the LEFM solution after a point which coincides with the start of crack propagation. For finite values of σ max , the gradual development of the cohesive zone ahead of the initial crack tip, before the crack starts propagating, is responsible for the nonlinear deviation of the loaddisplacement curve from the initial straight line of the LEFM solution, culminating in the rounded part so that, the lower σ max , the lesser the peak load with respect to  Comparison of a −F data during crack propagation (a closer look) and b relative errors when using Euler-Bernoulli instead of Timoshenko beam theory for different values of σ max the theoretical peak load predicted by LEFM. In terms of energy, this means that, in the case of a finite value of σ max , before the crack starts propagating, less external work is performed by the external force, F, than for the LEFM case. However, once the DCB is completely separated, the total amount of external work must be equal to the interface area times the area under the traction-separation law, that is the work of separation, . Therefore, for a DCB with an infinite length, the only way that the total external work spent is the same for the two cases is that, for a finite value of σ max , the curve lies above the case for σ max → ∞ during crack propagation, so that the increase in external work in this part of the curve compensates the lower amount of external work before crack propagation.
We can define a relative error due to using Euler-Bernoulli instead of Timoshenko beam theory as where, for the same value of σ max , F T and F E are the values of the applied force F at the same computed using Timoshenko and Euler-Bernoulli beam theory, respectively. By comparing Fig. 6b with Fig. 5b, we can see that this error has the highest values in the first (linear-elastic) phase. During the second phase the relative error decreases, and the higher the strength σ max the more rapid the reduction in the error. In the third phase, the error is negligible. We can also compare the relative differences in the results with respect to LEFM solutions by introducing factor ε r defined by: where F σ max and F ∞ are the values of the applied force F computed for the same crack mouth opening displacement for different values of σ max and for σ max = ∞, respectively. We can compute ε r separately for Timoshenko and Euler-Bernoulli beam theory. These results are presented in Fig. 7 for Timoshenko beam theory and Fig. 8 for Euler-Bernoulli beam theory. As expected, the highest values of the relative differences are obtained for more ductile interfaces (with smaller σ max ). However, the most important result is that these relative differences are extremely small during the crack propagation phase for both beam theories. This means that the most important parameter in the CZM, , can be computed accurately enough using F − data in the crack propagation phase and the simple LEFM formula (D.7), derived in "Appendix D" and based on Euler-Bernoulli beam theory assuming that = G c . This is because we have already demonstrated that the differences in the F − curves computed with Timoshenko and Euler-Bernoulli beam theories are negligible in the crack propagation phase (see Fig. 6b) for the case under examination. However, we have to keep in mind that for anisotropic materials (such as composites), where the shear modulus can be significantly smaller than Young's modulus, the differences between the results obtained with Timoshenko and Euler-Bernoulli beam theory could be more pronounced.
In Fig. 5 we noted that σ max significantly influences the peak load reached before the crack starts to propagate and dictates how far from brittle behaviour (LEFM) the considered CZM solution is. However, changing the ratio, α, between δ 0 and δ c has a noticeable influence on the results, too. We will assume that α can vary between a value very close to 0 (meaning that δ 0 is almost negligible compared to δ c ) and 1 (meaning that δ 0 = δ c ). The latter case, which is covered in Sect. 5.2, implies that we have a linear-elastic behaviour at the interface up to δ 0 (or δ c ) followed by brittle failure (leading to L cz = 0). In Figs. 9 and 10 we can see that increasing α reduces the stiffness of the interface. For α = 1 the phase of damage growth before the crack tip does not exist and the behaviour before crack propagation is linear-elastic. We can see that the influence of α on the results is more pronounced for σ max = 7.5 MPa than for σ max = 120 MPa, because, for the same , smaller σ max results in larger δ c . It can be noted, however, that α does not significantly influence the value of the peak force, especially for the brittle case.

Behaviour in front of the crack tip for a DCB with
linear-elastic interface with brittle crack and in the limit case of LEFM In Sect. 6 (see also "Appendix E") we showed that in the limit case of LEFM stresses and strains are found in front of the crack tip when Timoshenko beam theory is used to model the arms. In "Appendix E.3" we show that this is not the case when Euler-Bernoulli beam theory is used to model the arms. In this section, using the same geometrical and material data for the bulk material as in the previous examples, we will show that the behaviour of a DCB with a linear-elastic interface with brittle crack approaches the behaviour described in Sect. 6 for the limit case of LEFM as the stiffness of the interface increases. We will consider only the case of a DCB with prescribed displacement and investigate the case when the crack starts to propagate, which means that in Eqs. (107)-(112) we use F(a) (defined in (116)) and a instead of F and a 0 . However, in this example we will assume that a = a 0 = 30 mm, which means that we will investigate the case when the crack starts to propagate from its initial position. Note that, because the crack propagation for a DCB with prescribed displacement is not steady state, the presented results would change for a > a 0 , eventually approaching the steady-state solutions for a → ∞. These solutions are given by Eqs. In Fig. 11a we see that the cross-sectional rotations of the upper arm in front of the crack tip, ϕ 1 (x), reduce as the interface becomes stiffer, but for the limit case of LEFM, instead of approaching zero, they approach the limit values given by function (109). According to this formula, the rotation at the crack tip for the limit case of LEFM is ϕ(0) = − 0.0025 rad. In Fig. 11b we see how the contact tractions at the crack tip (σ (0) = σ max ) increase for higher values of σ max and tends to ∞ in the LEFM limit. In front of the crack tip, there is a distribution of compressive stresses which in the limit case of LEFM converges to the curve given by Eq. (110) according to which lim x→0 + σ (x) = − 105.98 MPa. This shows that the area of tensile stresses reduces to zero, while the stress at x = 0 tends to infinity. It is shown in "Appendix E" that the resultant of these tensile stresses tends to a finite value. In other words, in the LEFM limit, a Dirac distribution centred at the crack tip is found.   Figure 12 shows the distribution of shear forces in the upper arm in front of the crack tip, T 1 (x). We can see that for higher values of the stiffness the force rapidly drops from the positive value T 1 (0) = F(a) to a negative value. In the limit case of LEFM, according to (111), a discontinuity is found between lim x→0 − T 1 (x) = F(a) and lim x→0 + T 1 (x) = −8209.26 N. In Fig. 12b, the distribution of bending moments in the upper arm in front of the crack tip, M 1 (x), is shown. At the crack tip M 1 (x) = F(0)a, which for the limit case of LEFM, according to (112), is M 1 (0) = 25435.47 Nmm.
It is worth noting that using Euler-Bernoulli beam theory in the limit case of LEFM we have ϕ and M 1 (0) = F(a)a at the crack tip (x = 0).

Conclusions
In this paper we have derived complete analytical solutions for DCB specimens where the arms are modelled using simple beam theories (Timoshenko or Euler-Bernoulli). At the interface, three different models have been assumed: (i) a quasi-brittle bi-linear CZM, (ii) a liner-elastic CZM with brittle failure and (iii) a perfectly brittle and infinitely stiff interface (corresponding to LEFM solutions). The models obtained for each mentioned type of interface are called 'cohesive crack model' (CCM), 'enhanced beam theory' (EBT) and 'enhanced simple beam theory' (ESBT), respectively. In our approach EBT solutions are obtained from CCM solutions by removing the softening branch (responsible for progressive damage) from the CZM. From there, ESBT solutions can be obtained by letting the interface stiffness go to infinity. We have introduced a new term ESBT because we show that in the limit case of LEFM, EBT model does not correspond to 'standard beam theory' (SBT) model where the DCB arms act as if they were clamped at the crack tip. In ESBT, the arms are allowed to rotate at and in front of the crack tip, which is due exclusively to the shear deformability of the arms. This is why for the case when the arms are not sheardeformable (Euler-Bernoulli beam theory), ESBT corresponds to SBT. We have also derived the CCM, EBT and SBT solutions for the case when Euler-Bernoulli beam theory is used to model the arms by simply letting the shear stiffness of the arms go to infinity. All the mentioned solutions are derived for two different types of DCB: one with prescribed rotations and one with prescribed displacement. The presented solutions allow us to easily compute the crack mouth opening displacement, applied load (moment or force), contact tractions at the interface, displacement and rotations of the arms ahead of the crack tip, and shear forces and bending moments in the arms for different DCB models.
The main novel contributions of the paper are as follows: 1. The complete analytical solutions given in a unified and compact form for DCBs (either with prescribed rotations or prescribed displacement) with a bi-linear CZM at the interface. 2. The complete analytical solutions given in a unified and compact form for DCBs (either with prescribed rotations or prescribed displacement) with linearelastic interface with brittle failure. It is shown that such EBT solutions, compared to CMM solutions, are very accurate in the phase of crack propagation. However, they are not able to capture the quasibrittle behaviour before crack propagation. 3. The solutions for the limit case of LEFM which, in the case of Timoshenko beam theory (ESBT), show that at and in front of the crack tip the arms are allowed to rotate even if the interface is infinitely stiff. This implies that Timoshenko beam theory is, in a way, capable of capturing the realistic behaviour of the DCB in front of the crack tip and that the widely used assumption that in LEFM the arms act as if they were clamped at the crack is not necessary. However, we show that this assumption is still valid when the arms are modelled as Euler-Bernoulli beams. 4. An expression for the critical energy release rate, G c , for a DCB with prescribed displacement is pro-posed, which, to the best of authors' knowledge, is original. This expression, compared to the formula derived under the assumption that the arms are clamped at the crack tip, has an additional term which is dependent on the shear deformability of the arms and accounts for the rotations of the arms at the crack tip. If the arms are non-deformable in shear (Euler-Bernoulli beam theory), the expression for G c corresponds to the one obtained under the assumptions that the arms are clamped at the crack tip.
Future work will include the assessment of the accuracy of the formula for G c which in LEFM limit case takes into account the rotations of the arms at the crack tip. Similar analyses have been already done by Biel and Stigh (2008) and the authors of the present work (Škec et al. 2018). Based on the data obtained from the experiments or more sophisticated numerical models, our intention is to compare the accuracy of this formula with formulae from BS ISO 25217:2009(2009.
It is also worth noting that, using the approach presented in this paper, obtaining the analytical solution for a DCB with trapezoidal CZM at the interface using the Timoshenko beam theory to model the arms should be straight-forward and will also be covered in future work.

Free software made available
All the results are implemented in a software application with a user-friendly graphic interface where Euler-Bernoulli or Timoshenko beam theory for the arms and CCM, EBT or ESBT models for the interface can be selected. Results of the analysis can be plotted and exported. Because the computations are based on the presented analytical solutions, the results in our software are obtained instantaneously, even on a regular laptop computer. The software is free and can be downloaded at http://dx.doi.org/10.17633/rd. brunel.7223795.

Supplementary data
Supplementary material related to this article can be found online at http://dx.doi.org/10.17633/rd.brunel.

7212218.
Acknowledgements 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. The third author wishes to acknowledge the financial support of the Croatian Science Foundation (Research Project IP-2013-11-1631. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. In the phase of linear-elastic behaviour, Eqs. (24) and (25) become which, after substituting (20) for the case when ω = 1, gives From here, using the boundary conditions (30) we obtain which is the same result obtained in (34) for ω > 1 and ω < 1. This means that the moment M L can be expressed using (35) for any value of ω (including ω = 1). Because is equivalent to setting ω = 1 in (41), Eq. (42) is valid for any value of ω, and for ω = 1 it can be written as In the phase of damage growth before crack propagation, from boundary condition (43) we obtain From the continuity conditions (45) we obtain where the constants D i1 and D i2 read Parameter C 6 (L cz ) is defined as for C j (L cz ) in (55).
In the phase of crack propagation, constant C 5 is still obtained from (A.9) and constants D i from the expression analogous to (62), i.e. . . . , 4, (A.12) where C 6 is defined as for C j in (64). The crack mouth opening displacement (a) is computed according to (66).
Appendix A.2: DCB with prescribed displacement Using expressions (A.1) and (A.2), from boundary conditions (67) we obtain which gives as the limit value of the force .
The latter result implies that Eq. (71) is valid also for the case when ω = 1. Furthermore, (A.14) and show that expressions (70) and (76) are valid also for ω = 1. Thus, the crack mouth opening can be computed according to (77) as For the second phase of the solution we again obtain (A.9) and (A.10) from condition (43) and the continuity conditions (45), respectively. Constant C 6 (L cz ) is defined as for C j (L cz ) in (81), where the definitions (A.11) for constants D i1 and D i2 (i = 1, . . . , 4) must be used. We can express F(L cz ) using (80) and (L cz ) using (83).

Appendix B: Solutions for Euler-Bernoulli beam theory
In this section, we will assume that DCB arms are shear-undeformable by letting μ → ∞. This will simply transform solutions for Timoshenko beam theory presented in Sects. 3 and 4 into solutions for Euler-Bernoulli beam theory.
According to (9) 2 and (10) 2 , in the limit of μ → ∞ both ω and ψ tend to 0. This means that on the undamaged part of the interface the solution will be defined as for the case when ω < 1 in Sect. 2.2.1. Setting ω = ψ = 0, from (15) and (22), respectively, we now obtain The displacements of the upper arm on the undamaged and damaged part of the interface read The three-phase solution procedure reported in Sect. 3 is followed. For the liner-elastic phase we have The crack mouth opening displacement in the linearelastic phase can be expressed using (42) as For the phase of damage growth before crack propagation, we use (50) to determine constants D i (i = 1, . . . , 4), where now C j (L cz ) = C 3 (L cz ). From (44) we know that By taking into account (B.1), from (51) we get D 0 = 4, (B.8) and from (53) (55)  Equations (57) and (48) give the solution for M(L cz ) which is then used to numerically compute the value of L cz from Eq. (58) so that M max = M(L cz ). The crack mouth opening displacement in the second phase is computed from (61), where v 2 (L cz ) and v 2 (L cz ) are obtained from (B.3).
In the third phase we obtain constants D i (i = 1, . . . , 4) from (62) by using expressions (B.8) and (B.9) for D 0 , D i1 and D i2 and assuming that now C j = C 3 . Equation (64) now becomes The crack mouth opening displacement in the phase of crack propagation is computed according to (66).
In the third phase we use (85) to express C 3 (L cz ) as where D i1 and D i2 (i = 1, . . . , 4) are defined in (B.9). Function F(L cz ) is obtained from condition (78) and (49), whereas a(L cz ) is expressed using (87). Finally, the crack mouth opening displacement for the third phase can be computed according to (88).
where M max is defined in (93).

Appendix C.2: DCB with prescribed displacement
The crack mouth opening displacement in the linearelastic phase before crack propagation is computed according to (B.14). The peak load, as defined in (95), now becomes In the crack propagation phase, the crack mouth opening displacement is computed as where according to (97) we now have (20)) is considered in the limit case of LEFM, i.e.
v 1 (x) = e −λζ 1 x C 1 + e −λζ 2 x C 2 . (E.1) Constants C 1 and C 2 are defined according to (31) for a DCB with prescribed rotations, or according to (68) for a DCB with prescribed displacement. In the following subsections we will present the limit (LEFM) values of functions v 1 , v 1 , ϕ 1 , σ , T 1 and M 1 for both cases.
Appendix E.1: DCB with prescribed rotations Substituting (31) for C 1 and C 2 in (E.1) results in where the relation ζ 1 ζ 2 = 1 is used to simplify the expression. If we introduce the following notation where χ = E I /(2μA s ), we can rewrite (E.2) as v 1 (x) = M E I λ 2 (ρ 1 − ρ 2 ) ρ 1 e −λ 2 ρ1 x − ρ 2 e −λ 2 ρ2 x . (E.5) Multiplying (E.3) and (E.4) and because ζ 1 ζ 2 = 1, one then has that λ 2 ρ 2 = 1/ρ 1 , so that we can write: Note that in the exponent of the second term in the parenthesis we avoid λ 2 ρ 2 and use 1/ρ 1 instead because lim λ→∞ ρ 1 = 2χ = E I μA s , and lim λ→∞ ρ 2 = 0, (E.7) show that the limit value of λ 2 ρ 2 , unlike 1/ρ 1 , is less convenient for computer implementation. From (E.6) it then follows that From (E.10) it can be shown that the interface stress profile for finite values of λ starts from a positive value at x = 0, decreases to zero at a coordinate x = x 0 = ln(ρ 1 /ρ 2 )/ λ 2 (ρ 1 − ρ 2 ) , continues decreasing to a minimum, negative value, at a coordinate x min = 2x 0 , and then remains negative, but also increases and approaches zero asymptotically for x → ∞. For λ → ∞, x 0 → 0 and x min → 0. Stress profiles for different values of λ are reported in Fig. 11b. Note that the integral of the stress profile σ L (x) in (E.12) is finite: μA s E I (E.14) and therefore, for the vertical equilibrium to be satisfied the resultant F 0 = b lim c→0 c 0 σ L (x)dx of the infinite stress σ L (0) in (E.11) must be finite: where F 0 is a tensile concentrated force at the crack tip pointing downwards. Hence, the LEFM-limit stress profile, σ L , is not a proper function but a generalised one, which is the sum of (E.12) and of F 0 /b multiplied by a Dirac distribution D centred at zero: Using (E.9) we can obtain The interface stress profile is now very similar to the one defined in (E.12) for the case of a DCB with prescribed rotations. It can be shown that for (E.32) zero stress and minimum stress correspond to x 0 = 1 λ 2 (ρ 1 − ρ 2 ) ln ρ 1 (a 0 + ρ 1 ) ρ 2 (a 0 + ρ 2 ) , and x min = 1 λ 2 (ρ 1 − ρ 2 ) ln ρ 2 1 (a 0 + ρ 1 ) ρ 2 2 (a 0 + ρ 2 ) , (E.34) respectively, where as λ → ∞, both x 0 → 0 and x min → 0. Considering that Fa 0 in (E.33) and M in (E.13) both come from the same moment boundary condition -(30) 1 for the case of prescribed rotations and (67) 1 for the prescribed displacement -these two results are obviously the same, i.e. σ L (0 + ) = −μA s M 1 (0)/(bE I ). Again, in the limit case where λ → ∞, the resultant of the tensile stresses at the interface tends to a finite limit value F 0 . The integral of the limit stress profile σ L (x) in (E.32) has a finite value which reads where F 0 is a tensile concentrated force at the crack tip pointing downwards. Note that using (30) and (67) These expressions confirm that the DCB arms do not deform at and in front of the crack tip when they are modelled as Euler-Bernoulli beams. Analogous expressions for the phase of crack propagation could be obtained by substituting M with M max , F with F(a) and a 0 with a.