Determining mode I cohesive law of Pinus pinaster by coupling double cantilever beam test with digital image correlation

The direct identification of the cohesive law in pure mode I of Pinus pinaster is addressed. The approach couples the double cantilever beam (DCB) test with digital image correlation (DIC). Wooden beam specimens loaded in the radial-longitudinal (RL) fracture propagation system are used. The strain energy release rate in mode I ( I G ) is uniquely determined from the load-displacement ( P   ) curve by means of the compliance-based beam method (CBBM). This method relies on the concept of equivalent elastic crack length ( eq a ) and therefore does not require the monitoring of crack propagation during test. The crack tip opening displacement in mode I   I w is determined from the displacement field at the initial crack tip. The cohesive law in mode I I I ( ) w   is then identified by numerical differentiation of the I I G w  relationship. Moreover, the proposed procedure is validated by finite element analyses including cohesive zone modelling. It is concluded that the proposed data reduction scheme is adequate for assessing the cohesive law in pure mode I of P. pinaster.

ABSTRACT.The direct identification of the cohesive law in pure mode I of Pinus pinaster is addressed.The approach couples the double cantilever beam (DCB) test with digital image correlation (DIC).Wooden beam specimens loaded in the radial-longitudinal (RL) fracture propagation system are used.The strain energy release rate in mode I ( I G ) is uniquely determined from the load-displacement ( P   ) curve by means of the compliance-based beam method (CBBM).This method relies on the concept of equivalent elastic crack length ( eq a ) and therefore does not require the monitoring of crack propagation during test.The crack tip opening displacement in mode I   I w is determined from the displacement field at the initial crack tip.The cohesive law in mode I I I ( ) w   is then identified by numerical differentiation of the I I

INTRODUCTION
ood is a hierarchical, anisotropic and heterogeneous composite material formed by trees.Recently, green composites based on lignocellulosic fibres and forest-based resources have attracted increasing interest in both research and market [1,2].Moreover, in a policy of sustainability, wood and wood products are increasingly used nowadays, for instance, in structural and semi-structural applications [3].However, for a better and efficient utilisation of wood material, several issues must be further investigated.One fundamental aspect concerns the fracture mechanical behaviour of wood.Relatively extensive fracture process zones (FPZ) are observed in wood due to fibre bridging and micro-cracking ahead of the crack tip [4].However, the microstructural mechanisms in wood fracture are usually confined to a region of reduced thickness [5].Therefore, at the macroscopic scale, the wood behaviour in the FPZ can be conveniently described through a phenomenological constitutive cohesive law [6,7].In order to obtain the cohesive law, one approach consists in minimising an objective function quantifying the difference between numerical and experimental load-displacement ( P   ) curves by inverse analysis, assuming a given shape of the softening law.This approach, however, is semi-empirical and does not guarantee the uniqueness of the solution.Notwithstanding, it has been shown that the inverse identification of cohesive laws provide good agreement between experimental and numerical finite element simulations [7,8].Instead, a direct method for evaluating the cohesive law can be proposed based on independent determination of strain energy release rate and crack tip opening displacement (CTOD) [9].The advantages of this approach are: (i) the shape of the cohesive law does not need to be assumed a priori; (ii) the cohesive law is determined based on local measurements.In this work, a direct identification of the cohesive law in mode I of P. pinaster was investigated by coupling the double cantilever beam (DCB) test with digital image correlation (DIC).Specimens oriented in the radial-longitudinal (RL) propagation system were used.The strain energy release rate in mode I ( I G ) was explicitly determined from the P   curve by means of the compliance-based beam method (CBBM).This data reduction scheme is based on the concept of equivalent elastic crack length ( eq a ) and, therefore, does not require the measurement of the crack length during test.An independent evaluation of CTOD in mode I ( I w ) was determined from displacement fields at the initial crack tip.The direct differentiation of the I I G w  curve and the reconstruction of the I I w   cohesive law, by means of least-squares regression using a continuous approximation function, were addressed.The proposed procedure was also validated by finite element simulations including cohesive zone modelling.
where  accounts for root rotation effects and can be determined from finite element analysis [4], and G LR is the shear moduli of the material.The CBBM is based on eq a , which is considered to account for the FPZ effect at the crack tip as given by: eq FPZ a a a      [7,10].Finally, the application of CBBM to the DCB test yield the following expression for the strain energy release rate in mode I (resistance or R-curve) [4] It is worth noting that this procedure is less sensitive to experimental errors.This is supported by the fact that the measurement of the crack length during the fracture test is not required.Besides, the inherent elastic properties variability among specimens is taken into account by computing a flexural modulus for each P   curve Eq. ( 2).In mode I loading, strain energy release rate ( I G ) and CTOD ( I w ) can be related by the following expression [9] The cohesive law ) can then be directly obtained by differentiating the above equation This data reduction scheme, however, requires the accurate evaluation of the relationship.Moreover, a suitable differentiation algorithm must be used to avoid noise amplification in the reconstruction of the constitutive cohesive law.In order to solve Eq. ( 5), it is proposed here to fit the I I G w  data by a continuous function described by the following expression (logistic function) where 1 A , 2 A , p and I,0 w are constants to be determined in least-square sense.In this function, the 2 A parameter must provide an estimation of the critical strain energy release rate: In the literature, several analytical expressions for the cohesive law in mode I have been proposed.Among them, there is the Xu and Needleman exponential law [11] Ic in which Iu w is the crack tip opening displacement at maximum stress ( Iu  ) (see cohesive law in Fig. 2).The logistic (Eq.6) and exponential (Eq.7) cohesive laws are then compared and discussed during the analysis.

DIGITAL IMAGE CORRELATION
ull-field optical techniques have become very important tools in experimental solid mechanics.Among them, DIC has become widely used, following the development of digital cameras and automatic digital image processing techniques [12].This computer vision technique has the advantages of a simple principle and experimental set-up, which can switch from large down to small scales of observation.In DIC-2D, a planar object is imaged by a single camera-lens optical system connected to a computer for real-time visualisation.It is assumed that the surface of interest has a local random textured pattern uniquely characterising the material surface.A matching process is then carried out between images taken before and after deformation.Typically, the reference (undeformed) image is divided into subsets, whose number of pixels defines the displacement spatial resolution (i.e., the smaller distance separating two independent displacement measurements).The selection of these measuring parameters, together with the quality of the speckle pattern, are key issues for determining the spatial resolution and accuracy associated to DIC measurements.Therefore, they should be carefully chosen in a compromise between correlation (small subsets) and interpolation (large subsets) errors.

Crack tip opening displacement
The CTOD in mode I I ( ) w was determined by processing the displacements provided by DIC.For a complex material such as wood, this technique can be advantageous with regard to standard monitoring methods.As a procedure, the initial crack length is firstly identified in the reference (undeformed) image.A suitable pair of subsets near the crack tip is then selected.During the test, the relative displacement of these points is evaluated.The I w can then be determined by [13,14] where I w  and I w  are the components of the displacement in the direction perpendicular to the crack propagation associated, respectively, to the upper and the lower cracked surface and .represent the Euclidean norm.

Mechanical Properties
Fracture properties (Mode I, trilinear cohesive law

FINITE ELEMENT SIMULATION
wo-dimensional finite element analyses (FEA) of the DCB test under plane strain state, including cohesive zone modelling, were performed in order to validate the proposed procedure using CBBM.Nominal dimensions of the DCB specimen were (Fig. 1): 2h = 20 mm, 1 L = 300 mm, L = 280 mm, B = 20 mm and 0 a = 100 mm.F T Isoparametric 8-node planar solid finite elements were used for the bulk material, whilst 6-node cohesive elements were applied along the FPZ, which in this case was confined to a line at mid-height of the specimen (Fig. 2).Wood was modelled as an orthotropic linear elastic material with properties summarised in Tab. 1 [7,[15][16][17].

EXPERIMENTAL WORK
pecimens for the DCB test were cut from a single P. pinaster tree.Specimens were selected in the mature region within the stem at a given location to minimise material variability.Twelve specimens were prepared with axes oriented along the RL propagation system.The nominal dimensions of the DCB specimens were those already used in the numerical simulation (Fig. 1).The crack was generated by using a band saw of 1 mm thickness followed by an impacted blade to guarantee a sharp initial crack surface.The fracture tests were performed in an INSTRON 1125 universal testing machine, with a controlled cross-head displacement rate of 3 mm/min.The load was measured by means of a 5 kN load cell, setting the gain at 50 N/V.Specimens were tested after stabilisation at laboratory conditions of 60-65% relative humidity and temperature of 20-25 ºC.The ARAMIS DIC-2D system was used in this work [18][19][20].The speckled pattern required in the DIC method was painted over the region of interest by means of an airbrush to guarantee suitable granulometry, contrast and isotropy at the scale of magnification.For DIC analyses, a subset size of 15×15 pixel 2 (0.270×0.270 mm 2 ) and a subset step of 13×13 pixel 2 (0.234×0.234 mm 2 ) were selected for enhancing spatial resolution [19].A resolution in displacement of 1-2×10 −2 pixel (0.18-0.36 μm) was obtained.For measuring CTOD I ( ) w , the coordinates of the initial crack tip were firstly identified in the reference image.The I w was then determined by post-processing the displacement of subsets chosen upper and lower the crack tip during the test (Eq.8).The base distance between the two subsets was 0.468 mm .The algorithm proposed for this evaluation together with a further comparison between eq a and DIC a as a function of the applied displacement in the DCB test is described in detailed in Ref. [10].Hereafter, a systematic evaluation of mode I fracture properties (R -curve and cohesive law) obtained from both CBBM and Irwin-Kies equations is addressed.To start with, the R-curves determined from both Irwin-Kies equation and CBBM (Eq. 3) were analysed.The compliance versus crack length ( DIC a ) function was fitted in the least-square sense by a cubic function of the form: The analytical differentiation of this function was then used in order to compute I G .Fig. 3a shows the P   curves obtained experimentally together with the numerical one resulting from FEA using the trilinear cohesive law (Fig. 2).This law was defined using the average values of Ic G and Iu  in order to outline the area circumscribed by the law and the peak stress value, respectively.The coordinates of the inflection point ( Ib w = 0.04 mm, Ib  = 2.0 MPa) were taken from the experimental cohesive laws, considering the issuing average values.The scatter on the initial compliance among the curves ( 0 0.072 0.0076 C   mm/N) is expected due to the inherent variability of the material.Moreover, qualitatively, the numerical prediction of the P   curve was in good agreement with the experimental ones.In Fig. 3a it is also shown a macroscopic visualisation of crack propagation.As it can be seen, microcracking and fibre bridging can be identified.This confirms the difficulties in measuring accurately the crack length using conventional monitoring techniques.Some authors (e.g., [21]) report that the main mechanism of mode I fracture is fibre bridging.However, these observations suggest that both micro-cracking and fibre bridging contribute significantly for the energy dissipation in the FPZ.The R-curves in mode I obtained from the DCB test by both CBBM and Irwin-Kies equations are shown in Fig. 3b and 3c, respectively, together with the numerical resistance curve.The wide dispersion of the experimental curves is most likely a reflection of the local variability of wood microstructure at the initial crack tip (e.g., earlywood and latewood constituents).From the R-curves, the evaluation of the strain energy release rate in mode I was carried out at two distinct stages.The first corresponds to the starting point of the non-linearity in the P   curve and therefore the initialisation of the FPZ ( Ii G ), whilst the second is defined at the maximum loading ( Im G ). Due to the fact that some of the R-curves do not reveal a clear plateau identifying the critical strain energy release rate, it was assumed that:

Specimens
. This value is then related to the complete development of the FPZ and initial steady-state crack propagation [7,10].The Ii G and Ic G values for both CBBM and Irwin-Kies equations are reported in Tab. 2, together with density and flexural modulus (Eq.2).

T These results point out that an underestimation of I
G can be obtained if the actual crack length is used (Eq. 1) because a fraction of fracture energy dissipation at the FPZ is not properly taken into account.It should be noted that the density values among specimens have a coefficient of variation (C.V.) lower than 4%.Consequently, it is not surprising that scatter of Ii G and Ic G was not statistically correlated with density.Indeed, as already pointed out, this scatter is likely due to natural variability of wood cellular structure at the crack tip among specimens.As can be concluded, the evaluation of the strain energy release rate in mode I by the Irwin-Kies equation is slightly lower than the one from CBMM.This difference is of 11.8% and 12.9% for Ii G and Ic G respectively, which is lower than the coefficients of variation among the tested specimens.The higher scatter observed in the R-curves resulting from CBBM relative to Irwin-Kies is explained by the scatter on elastic properties.In fact, since CBBM is based on specimen compliance it becomes more susceptible to material variability.From the DIC measurements both normal and transverse CTOD, with regard to the crack plane, were determined during the DCB test.As expected, CTOD in mode II ( II w ) was negligible.Characteristic R-curves in mode I I I ( ) G w  were then obtained as shown in Fig. 4a and Fig. 5a for the CBBM and Irwin-Kies equations, respectively.The numerical characteristic R-curve obtained by FEA of the DCB test was generically in relative good agreement with the experimental ones.For determining the cohesive law (Eq.5), a logistic function (Eq.6) was firstly fitted to the experimental data as shown in Figs.4b (CBBM) and 5b (Irwin-Kies).For this analysis, only the data points just before initial crack propagation, where the FPZ is assumed to be completely developed were considered.As it can be seen, a relatively good approximation was obtained in both cases.This procedure allows filtering experimental data and provides a basis for analytical differentiation, which is less prone to noise amplification.The cohesive laws in mode I obtained from the DCB test are shown in Figs.4c (CBBM) and 5c (Irwin-Kies), together with the numerical curve.In this case, a systematic lower evaluation was obtained from this fitting with regard to the independent measurements based on the R- curves (Tab.2).It is interesting to notice that Iu  (Tab.3, 4) is of the same order of magnitude as the radial tensile strength (7.93 MPa) measured from tensile tests on P. pinaster (with density of 0.7 g/cm 3 ) in [22].The mean cohesive law in mode I for P. pinaster was determined from the mean values of the parameters governing the logistic function, as shown in Figs.4d (CBBM) and 5d (Irwin-Kies).For comparison purposes, the exponential (Eq.7) cohesive law is also plotted in Figs.4d (CBBM) and 5d (Irwin-Kies), determined from mean values of the fitting parameters.Qualitatively both laws have a similar behaviour as summarised in Tab. 3 and 4.
is schematically shown in Fig.1.The specimen is a 1 2 L h B   mm 3 rectangular beam.The resistance curve (R-curve) can then be determined from the Irwin-Kies equation compliance.From the Timoshenko beam theory and Castigliano theorem an expression for the compliance of the DCB specimen can be obtained.This equation can be solved for the flexural modulus ( f E ) using an initial compliance 0 ( ) C and the corrected initial crack length 0 ( ) a   as[4]

Figure 2 :
Figure 2: Finite element model of the DCB test (cohesive law, mesh and boundary conditions).

Figure 3 :
Figure 3: DCB test: (a) curves and macroscopic visualization of the crack propagation; (b) R-curve obtained by CBBM; (c) R-curve obtained by Irwin-Kies equation

Figure 4 : 3 (
Figure 4: DCB test processed by CBBM: (a) characteristic R-curves; (b) least-square regression with the logistic function; (c) cohesive laws in mode I; (d) comparison between logistic and exponential mean cohesive laws.The parameters ( 1 A , 2 A , p and I,0 w ) of the logistic function (Eq.6) obtained from this study are reported in Tab. 3 (CBBM) and 4 (Irwin-Kies), together with the characteristic values of maximum stress ( Iu  ) and relative displacements ( Iu w , Ic w in Fig. 2) in mode I.The mean value of 2 A (Tab. 3 and 4) represents an estimation of Ic G .In this case, a systematic lower evaluation was obtained from this fitting with regard to the independent measurements based on the R-

Table 1 :
Mechanical and fracture properties of P. pinaster used in the finite element analyses of the DCB test.

Table 3 :
Parameters of the logistic function ( 1 A , 2A , p and I,0 w ), and characteristic values of maximum stress ( Iu  ) and relative displacements ( Iu w , Ic w , see Fig.2) determined by CBBM.