In this work, an analytical model consisting of Timoshenko beams coupled with a cohesive zone is used in order to analyze the extent of the cohesive fracture process zone within mechanical systems displaying beam-like features. This is the case, for example, for the double cantilever beam (DCB) specimen or for the wedge test specimen. The predictions of the model are displayed under the form of diagrams and formulas involving non-dimensional parameters, that can be readily used for example in the context of the identification of cohesive parameters from experimental data. The predictions are also compared to those of previous papers using beam models. It appears that in some configurations, for example soft materials with moderate toughness, the kinematic hypothesis of the Timoshenko beam is necessary in order to get accurate estimate of process zone lengths. The influence of the shape of the cohesive law is also discussed.

Keywords: Fracture process zone; Cohesive zone model; Adhesion; DCB test; Wedge test; Process zone length

The description of fracture mechanisms by cohesive zone models is widely used by the micromechanics community. In particular, it enables a local description of the fracture process, that can incorporate features of the physical mechanisms at work at the micro or nano scale during material failure. In the cohesive zone model (CZM) description, a crack is considered a gradual process extending along a finite distance ahead of the crack tip, called the cohesive process zone (or cohesive zone). Some of the first works incorporating the application of CZM to model fracture is found in Dugdale ([15]), Barenblatt ([2]) and Rice ([35]). Further on, Needleman ([25], [26], [27]) used numerical implementations of the cohesive zone model in the framework of the finite elements method for modelling interface debonding while at the same time Carpinteri ([9], [10]) introduced the 'Cohesive Crack Model' terminology, and studied the softening and snap-back behaviour of cohesive solids. The applicability of CZM for dynamic fracture through the insertion of cohesive elements a priori into the FEM mesh was demonstrated in Xu and Needleman ([51], [52]). In contrast, Camacho and Ortiz ([7]) and Ortiz and Pandolfi ([29]) used an approach with cohesive elements added during the simulations for modelling dynamic fracture. In recent days, CZM has been largely applied to debonding and delamination problems where the crack path is generally known beforehand. Different shapes of CZMs have been proposed, namely the linear softening shape (Carpinteri et al. [11]; Williams and Hadavinia [48]), exponential forms Tvergaard ([45]) and Xu and Needleman ([50], [51]), trapezoidal forms Tvergaard and Hutchinson ([46]) among others. There have been recently some attempts to develop traction-separation relationship using molecular dynamics (MD) Spearot et al. ([41]), Yamakov et al. ([53]) and Zhou et al. ([55]), an area of research still ongoing. There have been several comprehensive reviews over the years among which de Borst et al. ([13]) and Needleman ([28]). Despite the decades of development in CZM, most often the input parameters are largely chosen for numerical reasons rather than physical insights.

CZM are mainly characterized by three physical parameters, that are the strength (critical traction at which the crack nucleates in the material) that will be denoted as $T_{\mathit{max}}$ in the following, the toughness (work of separation per unit surface in order to reach complete separation of the two interacting surfaces) that will be denoted $G_c$ and the critical opening displacement at full separation of the interface that will be denoted as ${\delta }_{\mathit{max}}$ . If the cohesive traction versus separation law shape is known a priori, then the three parameters could be reduced only to $T_{\mathit{max}}$ and $G_c$ .[1] These models can hence incorporate both nucleation and propagation of the cracks. They also obviously involve an internal length scale (proportional to $\frac{\mathit{EGc}}{T_{\mathit{max}}^2}$ ). It is important to note that, if the length of the process zone is necessarily related to this internal length, it also depends on the overall geometry of the specimen, and will be different if the system is an infinite body or a slender specimen like DCB.

CZM approach is used for a wide variety of materials, ranging from concrete Hillerborg et al. ([18]) to composites Schellekens and de Borst ([36]), Camanho et al. ([8]) and Turon et al. ([43]). Whereas in metals and ceramics, the extent of the fracture process zone is on the order of a few tens of nanometers Needleman ([25]), Spearot et al. ([41]) and Yamakov et al. ([53]), it is much larger in the case of polymers, typically a few tens of micrometers and it is then possible to make indirect measurements of the process zone length using optical techniques like digital image correlation (DIC) to track displacements fields in the vicinity of the crack front Döll and Könczöl ([14]) and Réthoré and Estevez ([34]). This allows for experimental identifications of the cohesive zone parameters. However, the use of a good mechanical model is necessary in order to accurately extract those parameters from experimental measurements. Along another line, a few studies provide a comparison between analytical expressions and numerical results for the estimation of the length of the process zone, namely Smith ([38]), Bao and Suo ([1]), Harper and Hallett ([17]), Turon et al. ([44]) and Soto et al. ([40]). Other prominent works address estimates of the length under different types of loading like pure mode I in orthotropic materials Sih et al. ([37]) and Yang et al. ([54]), or infinite media Irwin ([20]), Palmer and Rice ([30]) and Planas and Elices ([33]).

While most of the above discussed references deal with the calculation of a process zone size for infinite bodies, the estimate of the process zone size in thin structures, which are for example particularly common in automotive and aeronautics industries, is also an important issue. Double Cantilever Beam (DCB) and wedge tests are commonly used for experimentation, and we will focus on these tests. The work of Smith ([38]) is one of the first to report an estimate of a process zone length in a DCB test, using two different types of traction versus separation law shapes, i.e. (i) a Dugale type (with a uniform traction value over the whole process zone) (ii) a linear softening cohesive law. He used the Euler–Bernoulli beam theory to estimate the length of the cohesive zone before the onset of debonding. Following a similar trend, Suo and Bao ([42]) considered a numerical solution to estimate the process zone length for the same beam problem and related the process zone length with the energy release rate. Coincidentally, both of their results can be combined into one analytical relation for the process zone length:

Graph

where E is the Young's Modulus, h the beam thickness, $T_{\mathit{max}}$ the cohesive strength, $G_c$ the toughness and $M_I$ a parameter dependent on the traction versus separation law shape. For a Dugdale type cohesive law, both found $M_I=0.76$ and for the linear softening law, Smith ([38]) found $M_I=0.84$ . It is interesting to see that, alongside with the characteristic length of the material mentioned before, the thickness of the beam also plays a part in the value of the cohesive zone length.

The work of Massabò and Cox ([24]) considered mode-II failure modes in delamination for the end-notched flexural specimen. Soto et al. ([40]) proposed expressions for determining the length in a homogeneous orthotropic material with a crack growing under pure mode I or pure mode II. However, in all the above referenced works on slender structures, the cohesive zone length is generally over-predicted with errors of up to 1200%. Harper and Hallett ([17]) discussed the limitations of existing expressions and propose a factor of 0.5 to obtained improved accuracy. In contrast, Soto et al. ([40]) considered a least square fit to obtain a statistical solution and reduce the overall error.

Besides, most of these models neglect the shear effects on the beam kinematics by using simple bending-based or Euler–Bernoulli beam theories. However, unless the process zone length $L_{\mathit{cz}}$ is way larger than the beam's thickness, the classical beam theory is not valid anymore, and a Timoshenko beam kinematics should be used instead. This aspect is further elaborated in the upcoming section.

This paper is organized as follow: in Sect. 2, we introduce an new analytical model for DCB/Wedge test using a Timoshenko beam kinematic for the description of the mechanical behaviour of the two bonded beams coupled by a cohesive zone and including a description of the compression state existing between the two beams in the unbroken area just in front of the crack tip (see Fig. 2). This new model can be used to determine the cohesive zone length during loading, and particularly at the onset of the debonding. In Sect. 3, we use our analytical model to carry out an analysis of the process zone length, and diagrams are built using non-dimensional parameters. The cohesive zone length is expressed as a function of Young's modulus, strength, fracture toughness and height of the DCB beam. Two types of traction versus separation laws are considered: a constant law (i.e. Dugdale type) and a linear softening law (see Fig. 1). In Sect. 4, we provide a new formula for computing the cohesive zone length $L_{\mathit{cz}}$ similar to the one given by Smith ([38]) and Suo and Bao ([42]) for slender beams. In Sect. 5, the prediction of our model is compared with former predictions from the literature for both constant and linear softening law. Then the influence of the shape of the cohesive law is studied in order to determine to what extent it influences the identification of the cohesive zone parameters. Finally we conclude that neglecting the shear effect can have a strong impact on estimating the fracture process zone and the use of a Timoshenko kinematic to describe the mechanical behavior of the DCB/Wedge specimen is hence recommended.

Graph: Fig. 1 Traction versus separation constitutive laws used in our study: a constant (Dugdale type) cohesive law, b linear softening law

Graph: Fig. 2 Timoshenko Model to describe DCB Test: Region 1 is the free zone with load or displacement imposed to its end; region 2 is interaction zone where tx is the cohesive forces along the interface; region 3 represents the contact between the two beams

The reasons why a new analytical model is needed to describe mode I testing (as DCB or wedge test), despite several models already available in the literature, is now exposed. There are actually two families of model for DCB/wedge testing. The first one relies on using an elastic foundation near the root of the crack tip as an alternative to cohesive zone, as was considered starting in the works of Kanninen ([23]) and Williams ([47]) and then by Cottrell et al. Williams et al. ([49]) and Cotterell et al. ([12]) and more recently by Jumel and co-workers Budzik et al. ([4], [5], [6]) and Jumel et al. ([21], [22]). More relevant to our approach is the second family of models, based on beams coupled with a cohesive zone for the interface description. Euler–Bernoulli beam theory is based on hypothesis of a dominating pure bending of the beam, assuming that: (a) plane sections remain plane after deflection (b) plane sections that are initially perpendicular to neutral axis remain perpendicular after beam deflection. However, hypothesis (b) holds only in the cases of pure bending, when shear forces are of smaller magnitude than the bending forces. When the structures are slender, i.e., the ratio of height-to-length is high, then the Euler–Bernoulli theory, can be used. However when the structure is not slender (ratio beam height to beam length are comparable), it is necessary to consider a Timoshenko beam theory instead. Timoshenko beam theory assumes that plane sections do not remain perpendicular to the neutral axis after deflection, due to an additional rotation resulting from shear, that has to be considered in the equilibrium equation of the beam.

Some earlier works have considered the use of beam theory along with a cohesive zone model to describe the behaviour of the DCB specimens. Some of the prominent works involving beam theory along with a CZM to describe fracture of DCB specimens include Williams and Hadavinia ([48]) who explore different shape factors for the cohesive law, using Euler–Bernoulli beam theory with an elastic foundation for the contact zone. The authors conclude on the negligible effect of the cohesive zone shape on the energy release rate. Another work of interest is the study of the traction-separation relation for the silicon/epoxy interface considered by Gowrishankar et al. ([16]). The authors used interferometric measurements to estimate the normal opening displacement near a crack front, and computed a J-integral based analytical solution. They also used a direct approach of fitting with data from experiments to obtain the parameters for the traction-separation curve to conclude that both the approaches compared reasonably well. Similar to the work of Gowrishankar et al. ([16]) and Blaysat et al. ([3]) also considered an optical technique to extract the parameters of the traction-separation relations. Like many earlier works, Blaysat et al. ([3]) use a Winkler foundation for the undamaged region alongside a three-parameter cohesive zone model. Overall, the arms of the DCB are considered as two beams whose kinematics comply with the classical Bernoulli theory. However, it is important to note that several works, as early as Kanninen ([23]) and Williams ([47]), had concluded that the arms of the beams considered in the DCB specimens as Euler–Bernoulli beams could lead to significant errors. Furthermore, they concluded that the shear in the bonded region is comparable to the bending moment and thus a Timoshenko beam theory was the most appropriate option for modelling the arms of the DBC specimen.

In this work, we have considered a Timoshenko beam theory to model the two arms of the DCB/wedge specimen as shown in Fig. 2. Only half of the model is represented, for symmetry reasons. The plane of equation $y=0$ is therefore a symmetry plane for the complete specimen. The overall structure is considered to be composed of three regions:

• Region 1 ( $L_{\mathit{cz}}\le x\le L_{\mathit{cz}}+a$ ) represents surfaces that are fully separated and do not exhibit any further bonding or interfacial strength. It is composed of a beam of length a.
• Region 2 ( $0\le x\le L_{\mathit{cz}}$ ) is where the actual damage occurs, and the length of this region is the total length of the cohesive zone $L_{\mathit{cz}}$ . Over this region, the Timoshenko beam is submitted to a normal traction exerted by the cohesive zone.
• Region 3 ( $-\infty \le x\le 0$ ) is ahead of the crack tip. There is a compression along the symmetry plane ( $y=0$ ). Indeed, because of the symmetry, the displacement component $w_3(x)$ as well as all its derivatives have to vanish (the beam has zero normal opening). For the beam, that represents half of the real system, this is equivalent to a rigid contact condition along the symmetry plane. We introduce a pressure $P_c$ that represents the normal component of the traction vector along the plane of equation ( $y=0$ ), i.e. $P_c=-\overrightarrow{n}.(\overline{\overline{\sigma }}\overrightarrow{n})$ , with $\overline{\overline{\sigma }}$ the stress tensor in the full specimen and $\overrightarrow{n}$ the normal to the symmetry plane.

The beam lies in the x direction, and it is deflected in the y direction during deformation. The vertical displacement is denoted w(x). The origin $x=0$ is taken at the end of the cohesive zone, as indicated in Fig. 2. It is worth noting that the length of the cohesive zone is generally smaller than the height of the beam. It is an additional argument to justify the use a Timoshenko beam kinematics for the portion of beam in region 2, as it is fit for short beams.

On the left end of the beam, a loading is applied (either a prescribed force P in the case of the DCB or a prescribed displacement U in the case of a wedge test). The complete system of equations associated with the mechanical model is provided in Appendix A. We can mention here in brief that, in addition to the classical boundary and continuity conditions applied to the beam parts, an additional condition is used in order to determine the cohesive zone length, which is that the resultant of the loading applied to the whole beam (i.e. for $-\infty \le x\le a+L_{\mathit{cz}}$ ) in the vertical (y) direction has to vanish.

The complete derivation of the model including detailed comparisons with finite elements calculation and with other models will be presented in another paper. The main objective of the present paper is to provide practical formulas and data to be used for computing the process zone length, and to help cohesive zone parameters identification using experimental data.

In the previous section and in Appendix A, we present the analytical model that we used for the determination of the process zone length $L_{\mathit{cz}}$ in the DCB/wedge test specimen. The analytical solution that we obtained was validated using finite element simulations. They have been used in order to determine the length of the cohesive zone $L_{\mathit{cz}}$ depending on the loading (applied displacement U or applied load P), with a particular focus on the value of $L_{\mathit{cz}}$ at crack initiation. The crack initiation is defined at the loading at which $w(L_{\mathit{cz}})={\delta }_{\mathit{max}}/2$ , where ${\delta }_{\mathit{max}}$ is the opening of the cohesive zone at full separation.[2] For every fixed value of a, there is always a critical value of the loading such that the above relation is fulfilled. Those critical loading values can be denoted as follows:

2a $$\begin{array}{c}\hfill \begin{array}{cc}\hfill U_{\mathit{crit}}& =U\left(T_{\mathit{max}},,,G_c,,,E,,,\nu ,,,h,,,a\right)\text{for applied displacement}U\hfill \end{array}\end{array}$$

Graph

2b $$\begin{array}{c}\hfill \begin{array}{cc}\hfill P_{\mathit{crit}}& =P\left(T_{\mathit{max}},,,G_c,,,E,,,\nu ,,,h,,,a\right)\text{for applied load}P\hfill \end{array}\end{array}$$

Graph

The relevant parameters for the problem, both material and geometry related, are the Young's modulus E, the interface toughness $G_c$ , the interface strength $T_{\mathit{max}}$ and the height h of the beam. The free beam length a has an influence on $U_{\mathit{crit}}$ (or $P_{\mathit{crit}}$ ), but was shown to not affect the value of $L_{\mathit{cz}}$ .

Material and interface properties used for the non dimensional analysis

Graph: Fig. 3 Isovalues for Lcz/h considering a constant law cohesive zone domain for: a small values of X and Y, b large values of X and Y

Graph: Fig. 4 Isovalues for Lcz/h considering a linear damage cohesive zone domain: a small values of X and Y, b large values of X and Y

Assuming the Poisson's ratio $\nu$ to be a fixed parameter, we can propose the following non-dimensional expression for the process zone length:

3 $$\begin{array}{c}\hfill \frac{L_{\mathit{cz}}}{h}=f\left(\frac{E}{T_{\mathit{max}}},,,\frac{G_c}{T_{\mathit{max}}h}\right)\end{array}$$

Graph

We denote $X=\frac{E}{T_{\mathit{max}}}$ and $Y=\frac{G_c}{T_{\mathit{max}}h}$ the two non dimensional quantities chosen for the analysis.

Graph: Fig. 5 Variation of the parameters: aα and bβ for a constant (Dugdale type) cohesive zone

Graph: Fig. 6 Variation of the parameters: aα and bβ for a linear softening cohesive zone

Equation 3 can be easily plotted in a 2-D plane. About 1.5 million configurations are considered for different sets of Young's modulus E, interfacial strength $T_{\mathit{max}}$ , toughness $G_c$ and beam height h. As shown in Table 1, $T_{\mathit{max}}$ is varied between 10 and 1000 MPa, by increments of 10 MPa; E between 2 and 300 GPa by increments of 0.15 GPa; $G_c$ from 0.1 to 1 with increments of 0.1 N/mm; and h from 1 to 10 mm with increments of 1. Thus, the sample size consists of 100 different values for $T_{\mathit{max}}$ , 150 of E, 10 of $G_c$ and 10 of h and a total of 1.5 million sets of data. The Poisson's ratio is chosen equal to $\nu =0.3$ , while the length of the pre-existing crack is fixed at $a=50$ mm.[3] In order to evaluate the fracture process zone length $L_{\mathit{cz}}$ , Eq. 21 in Appendix A is used to determine the critical force $P_c$ which represents the force at the onset of the debonding. Then, the value of $L_{\mathit{cz}}$ is obtained after solving Eq. 20 in Appendix A with $P=P_c$ . This process is repeated for all the given configurations, after which the non-dimensional parameter $L_{\mathit{cz}}/h$ can be evaluated as a function of $X=\frac{E}{T_{\mathit{max}}}$ and $Y=\frac{G_c}{T_{\mathit{max}}\times h}$ . However, the configurations where the critical opening or the displacement at the end of the beam $U_c$ is bigger than beam's height are discarded, in order to respect the small displacement hypothesis assumed by our analytical model. Finally, 1.23 million remaining sets of data are considered for discussion.

The iso-values of the process zone length obtained from the analytical expression are plotted as a function of our non-dimensional material parameters X and Y. They are plotted for a constant cohesive law (Dugdale type) in Fig. 3 and for a the linear softening law in Fig. 4.

Figures 3 and 4 provide a practical way to estimate graphically the process zone length for a wide variety of parameters. Knowing one $T_{\mathit{max}}$ (resp. $G_c$ ) and measuring $L_{\mathit{cz}}$ , it is also possible to deduce $G_c$ (resp. $T_{\mathit{max}}$ ). The possibility of using an analytical formula to compute explicitly $L_{\mathit{cz}}$ from the values of X and Y, as it is proposed in other papers, is discussed in next section.

In the study considered in the previous section, the non-dimensional process zone length is expressed as a function of two non-dimensional parameters $\left(X,,,Y\right)$ under the form of graphical representation. While this is an interesting and accurate method to calculate the process zone length $L_{\mathit{cz}}$ directly from the model, it is somehow more appealing to have an explicit expression to use, for example similar to what is proposed in Suo and Bao ([42]) and Smith ([38]). Using exponent parameters $\alpha$ and $\beta$ , the process zone length can be expressed as :

4 $$\begin{array}{c}\hfill L_{\mathit{cz}}=h{\left(\frac{E}{T_{\mathit{max}}}\right)}^{\alpha }{\left(\frac{G_c}{T_{\mathit{max}}\times h}\right)}^{\beta }\end{array}$$

Graph

where parameters $\alpha$ and $\beta$ can be obtained by least-squares fit between our analytical solution and the one given by the formula, and for all the configurations. It is interesting to note that for $\alpha =\beta =1/4$ and for linear softening cohesive zone, we obtain the particular solution given by Smith ([38]) and Bao and Suo ([1]).

Figures 5 and 6 show the values of the parameters $\alpha$ and $\beta$ for constant and linear softening cohesive zone law respectively. As anticipated, the values of the parameters vary and thus it is not possible to predict the process zone length for all the configurations with a unique formula. The figures show the variation of the parameters $\alpha$ and $\beta$ for parameters X varying from 0 to 10,000 and Y from 0 to 0.004. The minimum, maximum and the average values are shown in Table. 2.

$\alpha$ and $\beta$ variations as function of X and Y

It is interesting to estimate where typical materials would fit within this diagram. If we take the example of PMMA with $E=2$ GPa, $T_{\mathit{max}}=100$ MPa, $G_c=350{\text{J/m}}^2$ and $h=50$ mm, we find $X=20$ and $Y=7\times {10}^{-5}$ , which corresponds to rather high values of $\alpha$ and low values of $\beta$ in Fig. 5. In contrast, if cast iron is considered, with $E=100$ GPa, $T_{\mathit{max}}=1000$ MPa, $G_c=20,000{\text{J/m}}^2$ and $h=50$ mm, we find $X=100$ and $Y=4\times {10}^{-4}$ . This corresponds to rather low values of $\alpha$ and high values of $\beta$ .

In the next section, predictions obtained with the above formula are compared to predictions from other publications.

In this section, both the results of the full analytical solution and of the dimensionless formula are compared with the existing results in the literature. For comparison, two types of cohesive zone laws, namely constant and linear softening laws, are considered. As shown in both these cases, the dimensionless results discussed in this work are, within reasonably acceptable bounds of error, able to predict the results obtained from our analytical solutions.

Graph: Fig. 7 Evaluation of the process zone length Lcz for a constant (Dugdale type) cohesive law, as a function of: aTmax, bE, ch and dGc. Comparison provided between the proposed analytical solution, dimensionless forumla with with α=0.4 and β=0.5, Bao and Suo ([1]) and Harper and Hallett ([17])

There are four primary parameters, namely maximum traction $T_{\mathit{max}}$ , Young's modulus E, fracture toughness $G_c$ and height of beam h, which govern the length of the process zone. Four different cases are considered to evaluate the sensitivity of each of these parameters. The values considered are $E=70,000$ MPa, $G_c=0.1$ N/mm, $h=4$ mm, $T_{\mathit{max}}=100$ MPa. When varied, the variation ranges of these parameters are:

• $T_{\mathit{max}}$ : 10 to 970 MPa with an increment of 10 MPa.
• E: 2 to 300 GPa with an increment of 2 GPa.
• h: 2 to 10 mm with an increment of 1 mm.
• $G_c$ : 0.1 to 1 N/mm with an increment of 0.1 N/mm.

For all these comparisons, the pre-cracked length $a=50$ mm and Poisson ration $\nu =0.3$ .

Graph: Fig. 8 Comparison of relative error between the proposed analytical solution using Timoshenko model, the expression of Smith ([38]), Bao and Suo ([1]) and Harper and Hallett ([17]) when varying a, TmaxbE.

The resulting variation of the process zone length as a function of the above four parameters are displayed in Fig. 7. The solutions that have used classical beam theory to predict the process zone length, like Suo and Bao ([42]) and Smith ([38]) significantly over-estimate $L_{\mathit{cz}}$ for all cases. Even when multiplying their predictions by a factor of 0.5 as Harper and Hallett ([17]) did, the solution is still over-estimated. While the solutions obtained are close, Fig. 8 shows that the error is still significantly large. In comparison, both the analytical solution and the approximate non-dimensional-based solution provide much better estimates. The non-dimensional-based solution with $\alpha =0.4$ and $\beta =0.5$ (mean values over all the sampled cases).

We can define the relative error as the discrepancy between the exact solution, given by our analytical model, and the estimation given by our formula (Eq. 4), the estimation of Bao and Suo ([1]), Smith ([38]) and Harper and Hallett ([17]). Hence, the expression of the relative error (RE) can be written as :

5 $$\begin{array}{c}\hfill RE[\%]=\left|\frac{LczExact-LczEstimated}{LczExact}\right|\times 100\end{array}$$

Graph

From the error plot in Fig. 8, the relative error for small process zone length, in the earlier proposed solutions of Bao and Suo ([1]), Smith ([38]) and Harper and Hallett ([17]), can be as large as 500%, when varying $T_{\mathit{max}}$ , and as large as 400%, when varying E. In contrast, the non dimensional solution proposed here, mostly has about 5% error except for very small values of E and $T_{\mathit{max}}$ where the error can go up to as much as 41%. In such cases, the full analytical solution of our model should be used (Figs. 3 or 5). If we take again the examples of PMMA and cast iron introduced in Sect. 4, with the same data values, we can evaluate the process zone length $L_{\mathit{cz}}$ from our model. Using the diagram of Fig. 3, we find $L_{\mathit{cz}}=1.1$ mm for PMMA and $L_{\mathit{cz}}=5.7$ mm and for cast iron. Putting the same material data and beam thickness in Eq. 1 given by Smith ([39]) and Bao and Suo ([1]), we find $L_{\mathit{cz}}=7.4$ mm for PMMA and 17.0 mm for cast iron, which represents approximately an error of 670 $\%$ and 297 $\%$ respectively when compared to our results.

The source of these discrepancies thus mainly lies in the fact that most of the earlier works in literature use a Euler–Bernoulli beam theory in contrast to the Timoshenko beam theory used here.

From the above discussions and comparison, it has been demonstrated that the proposed non-dimensional approximate Eq. 4 provides, with the $\alpha$ and $\beta$ values computed in previous section, improved results for most configurations of interest to engineers when compared with existing solutions in the case of slender beams (e.g. Smith [38]; Bao and Suo [1]).

A comparative study between our model and the different solutions available in the literature is discussed here, in the context of a linear softening law. For such a cohesive zone model, as shown in the earlier Table 2, the $\alpha$ and $\beta$ parameters considered are 0.425 and 0.475 respectively.

Graph: Fig. 9 Evaluation of the process zone length Lcz, using a linear softening law, as a function of: aTmax, bE, ch and dGc. Comparison provided between the proposed analytical solution, non dimensional formula with α=0.425 and β=0.475, Bao and Suo ([1]) and Harper and Hallett ([17])

The observations made in the previous section are also valid for the linear softening law. The process zone length stimated in Smith ([38]) and Harper and Hallett ([17]) exceeds the value computed with our model by 400 $\%$ . The discprency is clearly visible in Fig. 9.

Another important observation is that the process zone length calculated for the linear softening law is larger than the one obtained for the constant law, for the same parameters $T_{\mathit{max}}$ , E, $G_c$ and h. This illustrates how the shape of the cohesive law can have a strong impact on the extent of the cohesive zone length $L_{\mathit{cz}}$ . As $T_{\mathit{max}}$ and $G_c$ has been chosen to be the same as in the case of the constant law, the opening at full separation ${\delta }_{\mathit{max}}$ is twice the previous one. However, it is not obvious whether the extent the process zone size follows the same trend or not. This question is investigated in the next section.

Over the years, several works have addressed the issue of the prominent parameters related to the cohesive zone law. Some of the works, like those of Tvergaard and Hutchinson ([46]), Peter Feraren ([32]) and Parmigiania ([31]), discussed that the interfacial strength $T_{\mathit{max}}$ and critical energy release rate $G_c$ are critical parameters while the actual shape of the cohesive zone law itself is secondary. In contrast, Gowrishankar et al. ([16]) and many others, emphasized that the shape of the cohesive zone law is also of importance when studying the normal crack-tip displacement and propagation. They further discussed that optimal solutions are only possible when the shape is known a priori.

In order to compare the influence of the shape of the cohesive zone law on the variation of the process zone length, calculations are carried out for both the constant law and the linear softening law, with the same set of parameters E, $G_c$ , $T_{\mathit{max}}$ and h. In the last section, it was shown that the process zone length is larger when the cohesive law considered is linear, in comparison to the constant (Dugdale type) cohesive law. For quantitative comparison of the difference between the two process zone lengths obtained from the two types of cohesive zone laws, a new parameter $\eta$ , representing the difference in the process zone lengths between the linear softening and the constant laws, is defined such as:

6 $$\begin{array}{c}\hfill \eta [\%]=\frac{Lc{z}^{\mathit{Linear}}-Lc{z}^{\mathit{Constant}}}{Lc{z}^{\mathit{Constant}}}\times 100\end{array}$$

Graph

To this end, more than 1.23 million calculations using the analytical model are carried out, with sets of data covering a wide range of cases. Each configuration has a unique set of parameters $T_{\mathit{max}}$ , E, $G_c$ and h. For every unique configuration, the process zone length is evaluated twice, firstly with the constant law, secondly with the linear softening law, then the parameter $\eta$ is calculated for each configuration using Eq. 6. A histogram for the value of $\eta$ is shown in Fig. 10 for the 1.23 million configurations.

Graph: Fig. 10 Comparison between the process zone lengths obtained with linear softening and constant cohesive law

Statistics regarding the variation of the parameter $\eta$

The statistics related to the parameter $\eta$ are given in Tables 3 and 4 with the average or mean $\mu =\frac{{\sum }_{j=1}^N{\eta }_j}{N}$ , N is the number of studied configurations equal to 1.23 million.

Statistics regarding the parameter $\eta$

It seems clear from Table 3, and for all the studied configurations in terms of sets of data, that the process zone length obtained with the linear softening law is, in average, 51.94 $\%$ greater than the one obtained with the constant law, with a standard deviation of 4.66%, indicating a small discrepancy of the data.

After analyzing our statistical data, we found that $82.04\%$ of the configurations considered are within the first standard deviation while $90.75\%$ of the configurations are within the second standard deviation from the mean. This means that the process zone length is about 42–64 $\%$ greater when the cohesive zone model considered is linear rather than constant. This is not in the same proportion that the increase of separation at total separation (a ratio of about 1.5 of $L_{\mathit{cz}}$ compared to a ratio of exactly 2. for ${\delta }_{\mathit{max}}$ ). Thus, it can be concluded that the shape of the cohesive law is of significance when considering the fracture process. The shape of the law is critical in the identification of the cohesive zone parameters like interface strength $T_{\mathit{max}}$ , maximum opening ${\delta }_{\mathit{max}}$ and process zone length $L_{\mathit{cz}}$ .

Finally, from this section, we have learned that estimating the process zone length depends mainly on the materials properties, the geometry of the test and the shape of the cohesive zone. Our predictions using the iso-values and our formula cover effectively configurations of interest, that have not been previously reached by models discussed in literature. In addition, the shape of the cohesive zones should not be neglected when dealing with fracture tests. Our results have shown that the cohesive zone identification will depend on the shape of the traction versus separation law. Going from a constant law to a linear softening law, for the same $T_{\mathit{max}}$ , double the normal opening ${\delta }_{\mathit{max}}$ and increase the process zone length by 51.94 $\%$ in average.

DCB specimens are commonly used for the identification of CZM parameters and calibration of the CZM laws. Among others, the length of the process zone is of critical importance. In addition, it is often believed that the shape of the traction versus separation cohesive zone law does not play, by itself, an important part in the fracture process.

There are been several works related to the determination and calibration of the parameters of CZM. Many of these works are related and use the DCB specimen or Wedge tests. In this work, a brief literature review on the works has been carried out. One of the common approaches in the literature consists of modelling the DCB arms as Euler–Bernoulli beams. However, this assumption does not account for the shear in the fracture process region. In contrast, we have considered the arms of the specimen as Timoshenko beams to derive an analytical expression for the length of the process zone. The analysis allows for the determination of the process zone length in the case of DCB or wedge tests. The main aim of this work is to provide practical and easy-to-use data, under graphical form and also with a formula. A wide parameter space is covered, in terms of Young's modulus, fracture toughness, thickness of the beam and interface strength.

Finally, since the fracture process zone is often smaller than the beams' height in the case of DCB or Wedge tests, Euler–Bernoulli kinematics can lead to substantial errors while estimating the fracture process zone length or calibrating the cohesive zone parameters of the interface. In this case, a Timoshenko beam kinematics is more suitable to describe the mechanical behaviour of the DCB or wedge specimen, especially around the interface. In addition, we found that the cohesive law shape can play a key role in the identification of the cohesive parameters and we have showed that going from a constant to a linear softening law could increase the fracture process zone length by 52 $\%$ in average while doubling the critical opening displacement ${\delta }_{\mathit{max}}$ . This aspect is often disregarded in literature.

Future work will focus on providing a new analytical formulation, based on Timoshenko beam theory, for both pure mode II and mixed-mode. These formulations will be helpful in the identification of cohesive zone parameters in case of delamination problems.

The Timoshenko beam theory accounts for shear stress in the beam. The bending moment M and shear force T in the beam are expressed as:

7a $$\begin{array}{cc}& \begin{array}{cc}\hfill M& =EI\frac{\partial \varphi }{\partial x},\hfill \end{array}\hfill \end{array}$$

Graph

7b $$\begin{array}{cc}& \begin{array}{c}\hfill T=\kappa GS\left(\frac{\partial w}{\partial x},-,\varphi \right),\end{array}\hfill \end{array}$$

Graph

where w is the deflection, S is the area of the cross section, $\varphi$ is the rotation of the cross section, G is the shear modulus and $\kappa$ is a shear coefficient equal to 5/6 for a rectangular cross section (see Fig. 11).

With a distributed shear load q(x) along the beam, the equilibrium equations write:

8a $$\begin{array}{cc}\hfill \frac{\partial M}{\partial x}+T& =0,\hfill \end{array}$$

Graph

8b $$\begin{array}{cc}\hfill \frac{\partial T}{\partial x}+q& =0.\hfill \end{array}$$

Graph

Graph: Fig. 11 Timoshenko Beam kinematics

In region 1 ( $L_{\mathit{cz}}\le x\le a+L_{\mathit{cz}}$ ) of the DCB/Wedge test, only a load P or a displacement U is applied at the end of the specimen, with no distributed load applied to the beam ( $q(x)=0$ ). In this case, combining Eqs 7a and 8a simply leads to:

9 $$\begin{array}{c}\hfill EI\frac{{\partial }^4{w}_1(x)}{\partial x^4}=0,\end{array}$$

Graph

where $w_1(x)$ is the deflection of region 1 in z. Integration of this relation leads to an expression for $w_1$ :

10 $$\begin{array}{c}\hfill w_1(x)=A_1{x}^3+A_2{x}^2+A_{3}x+A_4,\end{array}$$

Graph

where $A_1,A_2,A_3\text{and}{A}_4$ are constants to be determined.

The cross section rotation along region 1 is also obtained:

11 $$\begin{array}{c}\hfill {\varphi }_1(x)=\frac{\partial w_1(x)}{\partial x}+C_1.\end{array}$$

Graph

Along Region 2 ( $0\le x\le L_{\mathit{cz}}$ ) a cohesive zone is inserted between the two DCB/Wedge specimen. Hence, traction forces will develop along the interface in response to the beam displacement, creating a distributed shear load q(x) along the beam. Combining Eqs 7a and 8a this time leads to:

12 $$\begin{array}{c}\hfill EI\frac{{\partial }^4{w}_2(x)}{\partial x^4}=q(x)-\frac{\mathit{EI}}{\kappa GS}\frac{{\partial }^{2}q(x)}{\partial x^2},\end{array}$$

Graph

and

13 $$\begin{array}{c}\hfill {\varphi }_2(x)=\frac{\partial w_2(x)}{\partial x}-\frac{T(x)}{\kappa GS}=\frac{\partial w_2(x)}{\partial x}+\frac{{\int }_0^{x}q(\xi )d\xi }{\kappa GS}.\end{array}$$

Graph

The traction-separation law governing the displacement of the beam along this region for a Dugdale type cohesive law, Dugdale ([15]) is characterized by $q(x)=-T_{\mathit{max}}$ . Replacing q(x) in Eqs. 12 and 13 respectively leads to:

14 $$\begin{array}{c}\hfill EI\frac{{\partial }^4{w}_2(x)}{\partial x^4}=-T_{\mathit{max}}-\left[\frac{\mathit{EI}}{\kappa GS},\frac{{\partial }^2(-T_{\mathit{max}})}{\partial x^2},=,0\right],\end{array}$$

Graph

and

15 $$\begin{array}{c}\hfill {\varphi }_2(x)=\frac{\partial w_2(x)}{\partial x}-\frac{T_{\mathit{max}}x}{\kappa GS}+C_2,\end{array}$$

Graph

where $C_2$ is a constant to be determined.

Integration of Eq. 14 yields:

16 $$\begin{array}{c}\hfill w_2(x)=\frac{-T_{\mathit{max}}{x}^4}{24D}+B_1{x}^3+B_2{x}^2+B_{3}x+B_4,\end{array}$$

Graph

where $B_1,B_2,B_3\text{and}{B}_4$ are constants to be determined.

As previously explained, in region 3 ( $0\le x\le -\infty$ ), there is a compression along the plane ( $y=0$ ), which is a symmetry plane for the complete specimen, the half of which is represented in the model. The expression of the contact pressure $P_c(x)$ that has to be applied in this region can be obtained by using a modified version of Eq. 12, replacing q(x) by $P_c(x)$ and $w_2$ by $w_3\equiv 0$ :

17 $$\begin{array}{c}\hfill P_c(x)-\frac{\mathit{EI}}{\kappa GS}\frac{{\partial }^2{P}_c(x)}{\partial x^2}=0\mathrm{for}x\in [0,-\infty ].\end{array}$$

Graph

Solving Eq. 17 yields:

18 $$\begin{array}{c}\hfill P_c(x)=C_3{e}^{\zeta x},\end{array}$$

Graph

where $\zeta =\sqrt{\frac{\kappa GS}{\mathit{EI}}}$ and $C_3$ is the last constant of the problem to be determined.

Finally, the complete solutions of the displacement and rotation fields in the systems reads:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{w}_1(x)=A_1{x}^3+A_2{x}^2+A_{3}x+A_4\hfill \\ {\varphi }_1(x)=\frac{\partial w_1(x)}{\partial x}+C_1\hfill \\ w_2(x)=\frac{-T_{\mathit{max}}}{24D}+B_1{x}^3+B_2{x}^2+B_{3}x+B_4\hfill \\ {\varphi }_2(x)=\frac{\partial w_2(x)}{\partial x}-\frac{T_{\mathit{max}}x}{\kappa GS}+C_2\hfill \\ P_c(x)=C_3{e}^{\sqrt{\frac{\kappa GS}{\mathit{EI}}}x}.\hfill \end{array}\right)\end{array}$$

Graph

The 11 coefficients $A_i,B_i(i=1,2,3,4)$ and $C_i(i=1,2,3)$ can be found from boundary conditions and continuity of displacement w(x), neutral axis rotation $\frac{\partial w(x)}{\partial x}$ , cross section rotation $\varphi (x)$ , bending moment M(x) and shear force T(x).

• At $x=a+L_{\mathit{cz}}$ :
• The bending moment is zero:
• $$\begin{array}{c}\hfill M_1(a+L_{\mathit{cz}})=\frac{EI\partial {\varphi }_1(a+L_{\mathit{cz}})}{\partial x}=0.\end{array}$$

Graph

• In the case of the wedge test, the displacement at the end of the beam is constant and equal to:
• $$\begin{array}{c}\hfill w_1(a+L_{\mathit{cz}})=U,\end{array}$$

Graph

• where U is half the thickness of the inserted wedge.
• Alternatively, in the case of the DCB test, the shear force at the end of the beam is equal to :
• $$\begin{array}{c}\hfill \frac{EI\partial {\varphi }_1^2(a+L_{\mathit{cz}})}{\partial x^2}=P,\end{array}$$

Graph

• where P is load applied to the end of the beam.
• At $x=L_{\mathit{cz}}$ :
• the following continuity conditions are used:
• Displacement: ${\displaystyle w_1(L_{\mathit{cz}})=w_2(L_{\mathit{cz}}).}$
• Cross section rotation: ${\displaystyle {\varphi }_1(Lcz)={\varphi }_2(Lcz).}$
• Bending moment: ${\displaystyle \frac{\partial {\varphi }_1(L_{\mathit{cz}})}{\partial x}=\frac{\partial {\varphi }_2(L_{\mathit{cz}})}{\partial x.}}$
• Shear force : ${\displaystyle \frac{\partial {\varphi }_1^2(L_{\mathit{cz}})}{\partial x^2}=\frac{\partial {\varphi }_2^2(L_{\mathit{cz}})}{\partial x^2}}$
• Beam profile slope: ${\displaystyle \frac{\partial w_1(Lcz)}{\partial x}=\frac{\partial w_2(Lcz)}{\partial x}}$ .
• At $x=0$ :
• The displacement and its first derivatives are zero:
• $$\begin{array}{c}\hfill w_2(0)=\frac{\partial w_2(0)}{\partial x}=0.\end{array}$$

Graph

• In addition, using Eqs. 7a and 7b, the cross section rotation can be expressed as:
• $$\begin{array}{c}\hfill {\varphi }_2(0)=\frac{\partial w_2(0)}{\partial x}+\frac{\mathit{EI}}{\kappa GS}\frac{\partial {\varphi }_2^2(0)}{\partial x^2},\end{array}$$

Graph

• which provides a way of getting the constant $C_2$ .
• Finally, using the continuity of the bending moment between region 2 and 3 provides the equation used to determine $C_3$ :
• $$\begin{array}{cc}\hfill \frac{\partial {\varphi }_2(0)}{\partial x}=& \frac{\partial {\varphi }_3(0)}{\partial x}=\left(\frac{{\partial }^2{w}_3(0)}{\partial x^2},-,\frac{q(x)}{\kappa GS}\right)\hfill \\ \hfill =& \frac{P_c(0)}{\kappa GS}.\hfill \end{array}$$

Graph

This makes 11 relations to determine the 11 coefficients $A_i$ , $B_i$ and $C_i$ .

In order to is estimate the length of the process zone $L_{\mathit{cz}}$ , one last condition is needed. Writing the equilibrium of the external loads applied to the beam in the y direction writes:

19 $$\begin{array}{c}\hfill \underset{-\infty }{\overset{0}{\int }}{P}_c(x)dx+\underset{0}{\overset{L_{\mathit{cz}}}{\int }}q(x)dx+P=0,\end{array}$$

Graph

where P is the load applied to the extremity of the beam. Taking for example $q(x)=-T_{\mathit{max}}$ in the case of a Dugdale type cohesive zone then leads to:

20 $$\begin{array}{c}\hfill \frac{C_3}{\zeta }-T_{\mathit{max}}{L}_{\mathit{cz}}+P=0.\end{array}$$

Graph

Equation 20 thus provides a relationship to determine $L_{\mathit{cz}}$ , for a given value of the loading parameter (U or P).

Finally, in order to find the critical value of the loading parameter at the onset of propagation ( $U_{\mathit{crit}}$ or $P_{\mathit{crit}}$ ), one last relation is used:

21 $$\begin{array}{c}\hfill w_2(L_{\mathit{cz}})=\frac{{\delta }_{\mathit{max}}}{2}.\end{array}$$

Graph

Finding reliable data on experimental cohesive zone lengths represents a challenge for the fracture mechanics community. The process zone length is usually very small compared to the test geometry, and it is difficult to obtain it by direct observation. However, long process zone lengths can be found in the case of fibrous composite materials where large scale bridging can be observed in front of the crack-tip during delamination. Thus, to validate our cohesive zone length estimation, we chose to apply our methodology, developed in Sect. 3, to Huang et al. ([19]) who performed experimental DCB tests on a Bamboo-based composite, known as parallel strand bamboo or PSB. The material properties and the test geometry can be found in Table 5. A fracture toughness $G_c=145{\text{J/m}}^2(0.145\text{N/mm})$ was calculated by the compliance method. A critical opening displacement ${\delta }_{\mathit{max}}=0.038$ mm was measured using a microscopic camera. The process zone length $L_{\mathit{cz}}$ was found to be equal to 12.1 mm.

We estimate the cohesive strength assuming a linear softening law: $T_{\mathit{max}}=2\times G_c/{\delta }_{\mathit{max}}=7.63$ MPa. Injecting these values into our non-dimensional parameters, we find $X=9000/7.63\approx 1179$ and $Y=Gc/T_{\mathit{max}}\times h\approx 8.{10}^{-4}$ . Using the iso-values of $L_{\mathit{cz}}/h$ plotted in Fig.4b for the linear softening law with the values X and Y, we find $L_{\mathit{cz}}/h=0.59$ which gives $L_{\mathit{cz}}=14,16$ mm for $h=24$ mm (17 $\%$ higher than the values provided by the authors). Using Eq. 1 with the linear softening law, we find $L_{\mathit{cz}}=$ 19,81 mm, (64% higher than the provided value). This illustrates how the Euler–Bernoulli based models can provide estimates that are noticeably different from our Timoshenko based model.

Materials properties and geometry of the PSB specimen used by Huang et al. ([19])

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.