Measuring the Cohesive Law in Mode I Loading of Eucalyptus globulus

Assessing wood fracture behavior is essential in the design of structural timber elements and connections. This is particularly the case for connections with the possibility of brittle splitting failure. The numerical cohesive zone models that are used to simulate the fracture behavior of wood make it necessary to assume a cohesive law of the material that relates cohesive tractions and crack opening displacements ahead of the crack tip. This work addresses the determination of the fracture cohesive laws of Eucalyptus globulus, a hardwood species with great potential in timber engineering. This study centres on Mode I fracture loading for RL and TL crack propagation systems using Double Cantilever Beam tests. The Compliance-Based Beam Method is applied as the data reduction scheme in order to obtain the strain energy release rate from the load-displacement curves. The cohesive laws are determined by differentiating the relationship between strain energy release rate and crack tip opening displacement. The latter is measured by the digital image correlation technique. High strain energy release rates were obtained for this species, with no big differences between crack propagation systems. The difference between the crack systems is somewhat more pronounced in terms of maximum stress that determines the respective cohesive laws.


Introduction
Hardwood species are increasingly used for structural purposes, and this is shown, for instance, by in the development of new products with great impact in the European market. In particular, Eucalyptus globulus Labill is seen as a hardwood species with major potential in timber engineering, because of its high mechanical performance and durability, its aesthetic qualities, and the large stock of eucalyptus resources. This situation requires continuous research in different fields, where improvement of the drying process and the development of laminated products are two of the main ongoing research focuses [1][2][3][4].
Fracture characterization of wood is of essential importance, especially in the design of timber elements and connections in engineering structures [5,6]. Connections are most particularly some of the most critical parts, as they may lead to a dangerous situation in cases of brittle splitting failure produced by tension perpendicular to the grain. The orthotropic average values for the radial modulus of elasticity E R = 1820 MPa, the tangential modulus of elasticity E T = 821 MPa, the shear modulus of elasticity in the LR plane G LR = 1926 MPa, and the shear modulus of elasticity in the LT plane G LT = 969 MPa, are taken from [24] using Galician Eucalyptus globulus with a similar density to the boards used in this study. These parameters were obtained by compression tests coupled with a stereovision system (DIC 3D). DCB specimens were prepared from these boards according to the specifications shown in Section 2.3.

Compliance-Based Beam Method (CBBM)
The procedure applied to determine the cohesive law corresponds to a direct method that requires establishing the relationship between the strain energy release rate in mode I loading (G I ), the crack tip opening displacement (w I ), and the traction tension (σ I ), according to Equation (1).
The cohesive law in mode I, defined as σ I = f(w I ), can be then determined by differentiating Equation (1), as follows: This requires the accurate measurement of G I evolution as a function of w I in the course of an experimental fracture test (in this case a DCB test, see details in Section 2.3). The classical data reduction schemes used for this purpose are based on beam theory or compliance calibration and require crack length (a) measuring during testing [25]. However, the fracture process zone (FPZ) ahead of the crack tip in wood involves toughening mechanisms, such as microcracking, crack-branching, or fiber-bridging, hindering the identification of the crack tip and therefore also the a-measurement.
To overcome this problem, the Compliance Based Beam method (CBBM) [20,26] is shown to be a suitable alternative. It is based on Timoshenko beam theory and it introduces the concept of an equivalent crack length (a eq ), accounting for the FPZ effect given by a eq = a + ∆ + ∆a FPZ . Accordingly, compliance for a DCB specimen during crack propagation can be written as where G LR is the shear modulus in the LR plane; B and h the specimen dimensions; and, E f the corrected flexural modulus (instead of E L ) to take into account the cross-section rotation effects at the crack tip during testing and local stress concentrations. E f can be estimated from Equation (4) when considering the initial compliance (C 0 ) and a corrected initial crack length (a 0 + ∆) where ∆ represents the Williams correction term given by [27] in the form: An iterative process can be used to solve Equations (4)-(6) until a converged value of E f is reached. It must be noted that E T and G LT values should be used instead of E R and G LR in Equations (3)-(7) when the TL crack propagation system is considered.
The equivalent crack length, a eq , that meets the specimen compliance recorded during propagation is evaluated from a polynomial function solved with Matlab ® (Mathworks, Madrid, Spain), according to [20].
Let us consider the Irwin-Kies equation [28] G I = P 2 2B dC da (7) the strain energy release rate in mode I (G I ) is obtained by combining Equations (3) and (7). It represents the resistance curve (R-curve) of the material to the crack growth.
The CBBM method has the definitive advantage of only requiring the experimental load-displacement (P-δ) curve to derive the evolution of G I without crack length monitoring, making it less sensitive to experimental errors. The G I is then correlated with crack tip opening displacement in mode I (w I ), measured by the digital image correlation (DIC) technique during the test (see details in Section 2.3) and its derivative yields in the cohesive law expressed in Equation (2). It is therefore important to accurately evaluate the G I = f(w I ) relationship. This was performed in two ways for further comparison and discussion: (a) a smoothing spline using Matlab ® was adjusted to the experimental curve in order to soften the noise before differentiation; (b) the G I -w I data were fitted, in the least-square sense, by a continuous approximation function (logistic function), as follows, where A 1 , A 2 , p, and w I,0 are constants determined by regression analysis. Although this function has no particular physical meaning, it is simply a tool for the analytical differentiation that is required to obtain the cohesive law. The A 2 parameter must provide an estimation of the critical strain release, as The direct approach presented in this data reduction scheme can be potentially extended to other fracture modes, including mode II by means of the end notched flexure (ENF) test [29,30].

Double Cantilever Beam (DCB) Test Coupled with Digital Image Correlation
Thirteen DCB specimens that were oriented along the RL crack propagation system and fourteen oriented along the TL system were prepared for fracture tests. The first letter indicates the loading direction (Radial and Tangential, respectively) and the second letter refers to the crack propagation direction (Longitudinal). The DCB specimens consist of a rectangular beam with L 1 × 2h × B mm (250 mm × 20 mm × 20 mm) nominal dimension, as schematically shown in Figure 1. A mid-height pre-cracked surface of 100 mm in length and 1 mm thickness was initially performed. This initial notch was then lengthened a few millimeters with a band saw in order to guarantee a sharp initial crack. The actual a 0 value for each specimen was measured after testing. A symmetrical pair of 3 mm diameter holes were drilled at 10 mm from the specimen end, where the load (P) perpendicular to the pre-cracked surface was applied. The applied load was transferred to the specimen by means of two 3 mm diameter steel pins that were inserted into the holes.

Double Cantilever beam (DCB) test coupled with digital image correlation
Thirteen DCB specimens that were oriented along the RL crack propagation system and fourteen oriented along the TL system were prepared for fracture tests. The first letter indicates the loading direction (Radial and Tangential, respectively) and the second letter refers to the crack propagation direction (Longitudinal). The DCB specimens consist of a rectangular beam with L1 × 2h × B mm (250 mm × 20 mm × 20 mm) nominal dimension, as schematically shown in Figure 1. A mid-height precracked surface of 100 mm in length and 1 mm thickness was initially performed. This initial notch was then lengthened a few millimeters with a band saw in order to guarantee a sharp initial crack. The actual a0 value for each specimen was measured after testing. A symmetrical pair of 3 mm diameter holes were drilled at 10 mm from the specimen end, where the load (P) perpendicular to the pre-cracked surface was applied. The applied load was transferred to the specimen by means of two 3 mm diameter steel pins that were inserted into the holes. Prior to testing, the specimens were conditioned at 20 °C and 65% relative humidity until equilibrium moisture content was reached. The mean value of moisture content was approximately 11%.
The fracture tests were carried out using an INSTRON 1125 universal testing machine (Instron, Barcelona, Spain) with a load cell having a maximum capacity of 5 kN and 50 N/V gain. Specimens were loaded under 3 mm/min displacement control.
Crack mouth opening displacement was recorded using the optical system ARAMIS DIC-2D (GOM mbH, Braunschweig, Germany) [31,32] (Figure 2). This is a non-contact system that applies the principles of digital image correlation (DIC). This technique shows clear advantages in comparison with traditional measurement methods since it makes it possible to measure the deformation field of a whole specimen area, providing more robust results. It is composed of an eightbit charge-coupled device (CCD) camera with a telecentric lens that was mounted on a translation bar for fine aligning of the optical axis with regard to the planar specimen surface. The specimens were illuminated by two cold light sources incorporated in the measuring device. A speckled pattern with black ink on a white matte surface is applied to the specimen by an airbrush IWATA, model CM-B (Anesta Iwata Iberica SL, Barcelona, Spain), so that proper granulometry contrast and isotropy at the magnification scale is ensured. The region of interest is focused on the area just in front of the crack tip, where the crack starts to propagate. The different components of the optical system and the measuring parameters selected for this work are compiled in Table 2. Complete P-δ curves were obtained in all tests with an acquisition rate of 5 Hz, while the acquisition of images from DIC was made with 1 Hz frequency. Prior to testing, the specimens were conditioned at 20 • C and 65% relative humidity until equilibrium moisture content was reached. The mean value of moisture content was approximately 11%.
The fracture tests were carried out using an INSTRON 1125 universal testing machine (Instron, Barcelona, Spain) with a load cell having a maximum capacity of 5 kN and 50 N/V gain. Specimens were loaded under 3 mm/min displacement control.
Crack mouth opening displacement was recorded using the optical system ARAMIS DIC-2D (GOM mbH, Braunschweig, Germany) [31,32] (Figure 2). This is a non-contact system that applies the principles of digital image correlation (DIC). This technique shows clear advantages in comparison with traditional measurement methods since it makes it possible to measure the deformation field of a whole specimen area, providing more robust results. It is composed of an eight-bit charge-coupled device (CCD) camera with a telecentric lens that was mounted on a translation bar for fine aligning of the optical axis with regard to the planar specimen surface. The specimens were illuminated by two cold light sources incorporated in the measuring device. A speckled pattern with black ink on a white matte surface is applied to the specimen by an airbrush IWATA, model CM-B (Anesta Iwata Iberica SL, Barcelona, Spain), so that proper granulometry contrast and isotropy at the magnification scale is ensured. The region of interest is focused on the area just in front of the crack tip, where the crack starts to propagate. The different components of the optical system and the measuring parameters selected for this work are compiled in Table 2. Complete P-δ curves were obtained in all tests with an acquisition rate of 5 Hz, while the acquisition of images from DIC was made with 1 Hz frequency.
with black ink on a white matte surface is applied to the specimen by an airbrush IWATA, model CM-B (Anesta Iwata Iberica SL, Barcelona, Spain), so that proper granulometry contrast and isotropy at the magnification scale is ensured. The region of interest is focused on the area just in front of the crack tip, where the crack starts to propagate. The different components of the optical system and the measuring parameters selected for this work are compiled in Table 2. Complete P-δ curves were obtained in all tests with an acquisition rate of 5 Hz, while the acquisition of images from DIC was made with 1 Hz frequency.   The crack tip opening displacement in mode I (w I ) was obtained by post-processing the displacements monitored by DIC. The initial crack length was firstly identified in the undeformed image. The relative displacement between a pair of subsets selected close to the crack tip is evaluated afterwards. The value of w I is calculated as the Eucledian norm, as shown in Equation (9) [33,34].
where w + I and w − I are the displacement components in the direction perpendicular to the crack propagation associated to the upper and the lower cracked surface, respectively. This approach has the limitation of defining the cohesive law based on surface measurements, which may not be fully representative of the crack front over the volume of the FPZ. Current research interests have been expanded to the experimental observation of the volumetric crack front ahead of the crack tip by X-ray computed tomography [35,36] and exploring paths, which includes digital volume correlation, for the quantitative extraction of relevant mechanical [37,38] and fracture parameters [39].

Resistance Curve from CBBM
The load-displacement curves obtained from the DCB specimens for RL ant TL crack propagation systems are shown in Figure 3.

Resistance curve from CBBM
The load-displacement curves obtained from the DCB specimens for RL ant TL crack propagation systems are shown in Figure 3.  Both groups of P-δ curves show quite consistent results with variation typical of wood. The non-linear behaviour that was observed before the curves peak reveals that a non-negligible FPZ develops ahead of the crack tip. This phenomenon is characteristic of quasi-brittle materials, like wood and results in micro-cracking and fiber bridging, as confirmed by the macroscopic image in Figure 4. These observations are in agreement with other authors' research in wood [19,25]. This fact supports the difficulties in measuring the crack length during testing using conventional techniques and the convenience of applying the alternative CBBM method, when considering an equivalent crack length. Both groups of P-δ curves show quite consistent results with variation typical of wood. The nonlinear behaviour that was observed before the curves peak reveals that a non-negligible FPZ develops ahead of the crack tip. This phenomenon is characteristic of quasi-brittle materials, like wood and results in micro-cracking and fiber bridging, as confirmed by the macroscopic image in Figure 4. These observations are in agreement with other authors´ research in wood [19,25]. This fact supports the difficulties in measuring the crack length during testing using conventional techniques and the convenience of applying the alternative CBBM method, when considering an equivalent crack length.
The maximum loads that were attained in the tests are shown in Tables 3 and 4 for RL and TL crack propagations systems, respectively. The mean value from TL specimens is slightly lower than the one obtained from RL, but all of the values vary within a similar range. This means that the high proportion of radially oriented rays acting as reinforcement in this direction typical of hardwoods, like beech and ash [40], is not shown so pronouncedly in eucalyptus. However, microstructural constrictions may act in both directions. In most specimens, the decrease in load did not occur suddenly after reaching the peak load, as it arose gradually after a considerable increase in displacements, displaying an uneven behavior, due, once again, to the complex structure of wood [26,40]. The initial compliance C0 is calculated using Matlab ® , as this is the result that provides the maximum R 2 in every P-δ curve. A representative example of a DCB specimen is shown in Figure 5.
The C0 values resulting from all the tests are also included in Tables 3 and 4 for both crack propagation systems. There are minor quantitative differences between both orientations, while scatter is within the expected range for wood. The R-curves that were obtained from DCB tests in RL and TL crack propagation systems were evaluated by applying the CBBM and they are shown in Figure 6. From them, following an initial rising domain characterized by the development of the FPZ (corresponding to the non-linearity The maximum loads that were attained in the tests are shown in Tables 3 and 4 for RL and TL crack propagations systems, respectively. The mean value from TL specimens is slightly lower than the one obtained from RL, but all of the values vary within a similar range. This means that the high proportion of radially oriented rays acting as reinforcement in this direction typical of hardwoods, like beech and ash [40], is not shown so pronouncedly in eucalyptus. However, microstructural constrictions may act in both directions. In most specimens, the decrease in load did not occur suddenly after reaching the peak load, as it arose gradually after a considerable increase in displacements, displaying an uneven behavior, due, once again, to the complex structure of wood [26,40].
The initial compliance C 0 is calculated using Matlab ® , as this is the result that provides the maximum R 2 in every P-δ curve. A representative example of a DCB specimen is shown in Figure 5. The C 0 values resulting from all the tests are also included in Tables 3 and 4 for both crack propagation systems. There are minor quantitative differences between both orientations, while scatter is within the expected range for wood.  Both groups of P-δ curves show quite consistent results with variation typical of wood. The nonlinear behaviour that was observed before the curves peak reveals that a non-negligible FPZ develops ahead of the crack tip. This phenomenon is characteristic of quasi-brittle materials, like wood and results in micro-cracking and fiber bridging, as confirmed by the macroscopic image in Figure 4. These observations are in agreement with other authors´ research in wood [19,25]. This fact supports the difficulties in measuring the crack length during testing using conventional techniques and the convenience of applying the alternative CBBM method, when considering an equivalent crack length.
The maximum loads that were attained in the tests are shown in Tables 3 and 4 for RL and TL crack propagations systems, respectively. The mean value from TL specimens is slightly lower than the one obtained from RL, but all of the values vary within a similar range. This means that the high proportion of radially oriented rays acting as reinforcement in this direction typical of hardwoods, like beech and ash [40], is not shown so pronouncedly in eucalyptus. However, microstructural constrictions may act in both directions. In most specimens, the decrease in load did not occur suddenly after reaching the peak load, as it arose gradually after a considerable increase in displacements, displaying an uneven behavior, due, once again, to the complex structure of wood [26,40]. The initial compliance C0 is calculated using Matlab ® , as this is the result that provides the maximum R 2 in every P-δ curve. A representative example of a DCB specimen is shown in Figure 5.
The C0 values resulting from all the tests are also included in Tables 3 and 4 for both crack propagation systems. There are minor quantitative differences between both orientations, while scatter is within the expected range for wood. The R-curves that were obtained from DCB tests in RL and TL crack propagation systems were evaluated by applying the CBBM and they are shown in Figure 6. From them, following an initial rising domain characterized by the development of the FPZ (corresponding to the non-linearity beginning in P-δ curve), resistance to crack growth tends to a horizontal asymptote, despite the noise in the measurements, which defines the value of critical strain energy release rate (GIc) and it represents the material´s toughness to crack-growth. The R-curves that were obtained from DCB tests in RL and TL crack propagation systems were evaluated by applying the CBBM and they are shown in Figure 6. From them, following an initial rising domain characterized by the development of the FPZ (corresponding to the non-linearity beginning in P-δ curve), resistance to crack growth tends to a horizontal asymptote, despite the noise in the measurements, which defines the value of critical strain energy release rate (G Ic ) and it represents the material´s toughness to crack-growth.    In this research, most of the specimens had plateaus for a given crack extent, which means that the FPZ has been completely developed. Therefore, the fracture energy could be determined as a mean value over the horizontal domain of the curves. In cases where the R-curve did not show a clear plateau, the strain energy release rate corresponding to the maximum load (G I,Pmax ) could be assumed as critical strain energy release rate. Both values are reported in Tables 3 and 4, together with the maximum load reached in the tests, the corrected flexural modulus of elasticity from every specimen and the initial compliance. The last three values are input parameters in the CBBM formulation detailed in Section 2.2. The wide dispersion of R-curves may be due to local variability of wood microstructure at the crack tip among the specimens, e.g., earlywood and latewood [26].
The mean G Ic values that were obtained from both crack propagation systems were found to be the same: 0.77 N/mm. This value is considerably higher than that for other species. In particular, DCB specimens of Pinus pinaster gave a mean G Ic value of 0.31 N/mm when applying the same data reduction method in [25]. Eucalyptus globulus also displays higher fracture energy than other hardwood species. For instance, average G f values of 0.48 and 0.40 N/mm were obtained for beech and ash, respectively, in previous work by the author [40].
As the differences between G Ic and G I,Pmax were minimal (see Tables 3 and 4), the latter can be taken as a practical measure of mode I critical strain energy release rate.

Cohesive Law
G I evolution was correlated with the crack tip opening displacement (CTOD) values that were provided by DIC during testing in order to obtain the cohesive law in mode I. Normal and transverse CTOD with respect to the crack plane were determined. Representative normal and transverse CTOD-δ curves are shown in Figure 7. As can be seen, CTOD in mode II (w II ) was found to be negligible in DCB mode I tests.
GI evolution was correlated with the crack tip opening displacement (CTOD) values that were provided by DIC during testing in order to obtain the cohesive law in mode I. Normal and transverse CTOD with respect to the crack plane were determined. Representative normal and transverse CTOD-δ curves are shown in Figure 7. As can be seen, CTOD in mode II (wII) was found to be negligible in DCB mode I tests.
Characteristic GI-wI curves were then obtained, as shown in Figures 8 and 9.   Characteristic G I -w I curves were then obtained, as shown in Figures 8 and 9.

Cohesive law
GI evolution was correlated with the crack tip opening displacement (CTOD) values that were provided by DIC during testing in order to obtain the cohesive law in mode I. Normal and transverse CTOD with respect to the crack plane were determined. Representative normal and transverse CTOD-δ curves are shown in Figure 7. As can be seen, CTOD in mode II (wII) was found to be negligible in DCB mode I tests.
Characteristic GI-wI curves were then obtained, as shown in Figures 8 and 9.   The cohesive law for each specimen was finally determined by fitting a logistic function to the experimental data, as shown in Figure 10. The parameters corresponding to the logistic function expressed in Eq. 9 (A1, A2, p, and wI0), the area circumscribed by the cohesive laws (Glaw,I), the maximum stress (σI,u), and the relative displacements (wIu and wIc) are all included in Tables 5 and 6 for RL and TL crack propagation systems, respectively. The cohesive law for each specimen was finally determined by fitting a logistic function to the experimental data, as shown in Figure 10. The parameters corresponding to the logistic function expressed in Equation (9) (A 1 , A 2 , p, and w I0 ), the area circumscribed by the cohesive laws (G law,I ), the maximum stress (σ I,u ), and the relative displacements (w Iu and w Ic ) are all included in Tables 5  and 6 for RL and TL crack propagation systems, respectively.
least-square regression with the logistic function (right).
The cohesive law for each specimen was finally determined by fitting a logistic function to the experimental data, as shown in Figure 10. The parameters corresponding to the logistic function expressed in Eq. 9 (A1, A2, p, and wI0), the area circumscribed by the cohesive laws (Glaw,I), the maximum stress (σI,u), and the relative displacements (wIu and wIc) are all included in Tables 5 and 6 for RL and TL crack propagation systems, respectively.     35 32 In these tables, the mean value of A 2 depicts an estimation of the critical strain energy release rate, G Ic , and it acquires the values of 0.78 and 0.76 N/mm in the RL and TL crack propagation systems, respectively. These values are quite close to the mean ones obtained from the horizontal asymptote of the R-curves, that is 0.77 N/mm in both propagation systems (see Tables 3 and 4).
The mean parameters of Tables 5 and 6 are used to build the mean experimental cohesive law in mode I that can be considered for Eucalyptus globulus in RL and TL crack systems. It is highlighted by a bold curve in Figure 10. As can be seen, the difference between both propagation systems is more pronounced in terms of cohesive law than it is for the other previously derived fracture parameters. Mean maximum stress in RL is approximately 50% higher than in the TL system. In general, the maximum stresses that were displayed by the RL specimens are reached with lower crack tip opening displacements than in TL.
The mean cohesive laws in mode I obtained for eucalyptus could be implemented in numerical cohesive zone models to simulate the development of the FPZ and crack growth, and thereby analyze the actual fracture behavior of a timber structure.

Conclusions
The cohesive laws in mode I of Eucalyptus globulus for RL and TL crack propagation systems were determined by means of DCB fracture tests. The CBBM data reduction method was applied to derive the strain energy release rate (G I ) from the load-displacement curves when considering an equivalent crack length (a eq ) instead of the actual one, which would be difficult to measure. The G I was correlated with the crack tip opening displacements measured by digital image correlation technique to obtain the cohesive laws.
The same G Ic value of 0.77 N/mm was obtained for RL and TL crack propagation systems from the horizontal asymptote of the R-curves. The estimation of G Ic from the logistic function parameters that was used to attain the cohesive law was also within the same range (0.78 and 0.76 N/mm in RL and TL, respectively).
The behavioral difference between the two crack propagation systems was more significantly displayed in the cohesive laws. In this sense, RL laws showed higher values of mean maximum cohesive stress at lower crack tip opening displacements in comparison with the TL system. The cohesive laws definition in both crack propagation systems makes it possible to implement them in numerical cohesive zone models to accurately simulate crack growth along the FPZ and quantify the actual fracture behavior of the timber structure in question, especially in elements and connections involving the possibility of brittle splitting failures.
The excellent mechanical properties added to the high fracture toughness shown in this study underline that Eucalyptus globulus L. is a hardwood species of great interest for structural design.