INTRODUCTION

Nowadays, wooden structural elements are increasingly often made of GLULAM – a load-bearing element consisting of several to a dozen or so beams joined together with a thin layer of glue (Fig. 1). Glued beams have different layers of early and late wood (known as growth rings).

It should be noted that GLULAM structural elements contain cracks associated with material discontinuity in wood. They can occur naturally as a result of damage to the internal structure of the tree and the appearance of intercellular cracks, and can also arise during production processes. Such defects are stress concentrators, weakening the structure and increasing the possibility of its destruction. Therefore, damage to a single layer of GLULAM caused by defects weakens the entire structure. Therefore, developing a methodology for determining the critical parameters of the cracking process is a current scientific problem.

The use of linear elastic fracture mechanics (LEFM) for wood materials based on the area method which allows the strain energy release rate to be determined using only the P-δ graph was proposed in [1,2]. It presents the results of testing three-point bending specimens with the crack propagation direction along the force application axis [1] and specimens subjected to MODE I loading [2].The results obtained showed that the field method can be used to determine the strain energy release rate in wood.

Critical values of the strain energy release rate G_Ic,G_IIcof wood samples using acoustic equipment were determined in [3]. The tests were conducted for MODE I and MODE II, as well as for crack propagation systems (RL and TL). Critical values were determined at the point of crack propagation initiation (loss of linear force increment). It has been established that the critical values of the strain energy release rate in both loading modes found for the RL system are similar to the corresponding values for the TL system. In each of the systems (RL or TL), the critical value of G_IIc is twice as high as G_Ic. It should be noted that the signal spikes obtained during the test from the sound sensors located after the sample reached its maximum force P_max may suggest that the fracture process has stabilized at this point.

The current state of knowledge on the cracking of samples made of pine (Pinus sylvestris) is presented in [4-7]. In [3], the crack was cut parallel (RL) or perpendicular (TL) to the annual rings. The study was conducted using the CBBM method based on beam theory. The strain energy release rate was determined using the compliance function C(a), where a is the length of the crack. The results were obtained for the tensile test (MODE I) using beam specimens fixed with Double Cantilever Beam (DCB) pins. It has been established that the ratio of half the height h of the sample to the initial crack length a₀ has a significant effect on the critical value of the Released-Energy Rate G_Ic, the maximum normal stress σ_n^max and the length l_cz of the fracture process zone before the crack tip.

The impact of two load tests (MODE I and MODE II) on the strain energy release rate G_I i G_II, maximum normal stress σ_nand tangential stress τ_t [4] was compared [5]. As in [4], were performed for RL and TL systems using the beam method. DCB samples (for MODE I tests) and three-point bending samples (ENF) (for MODE II tests) were used. It was determined that the maximum stress values for both types of loads were higher for the RL system. However, the critical value of the strain energy release rate and the length of the crack propagation zone were higher for the TL system. It should be noted that the values of all parameters tested for samples subjected to MODE II tests were at least twice as high as the corresponding values obtained in MODE I tests.

Unlike in [4, 5], where the crack length a was determined analytically based on beam theory, in [6] this parameter was determined experimentally. DCB wood samples in an RL configuration were subjected to four loading cycles (MODE I) [6].The strain energy release rate was determined using three methods: area, compliance calibration (CC), and modified beam theory (MBT) [8-10]. The average critical values of strain energy release rate are identical when using each of these three methods. However, the spread of the results obtained using these methods was significant – the greatest in the case of the Area method, while for the CC and MTB methods they were almost identical.

It should be emphasized that in works [1,2], [4-7] and [10,11], the critical value G_crit of the strain energy release rate is the value at which stable fracture occurs. In contrast, in work [3], the value G_crit was obtained for crack initiation.

In summary, it can be concluded that no information has been found in scientific literature to date relating to the effect of the direction of early wood and latewood layers on the critical value of the energy release rate and crack tip opening displacement (CTOD) of DCB samples. However, the effect of input parameters on the characteristics of the strain energy release rate was not presented. Scientific studies most often present the results of experimental wood fracture tests for cracks running tangentially (RL) or perpendicularly (TL) to the wood annual rings. The aim of this study is to develop a model of pine wood cracking in the event of GLULAM layer damage when the crack plane is inclined at an angle to the wood annual rings.

RESEARCH METHODOLOGY

Description of samples and test rig

The experiment used 57 DCB samples measuring 20 mm×20 mm×300 mm in which the orientation of the early and late wood layers varied between 0° (RL) and 90° (TL) (Fig. 2). The sample material was pine (Pinus sylvestris) with an initial moisture content of 65% at a temperature of 20°C. Before testing, these samples were conditioned to achieve a relative humidity of 12%. A crack tip with a length of a₀ = 105 mm was made in the samples from the point where the pins were attached. A preliminary cut 100 mm long was made using a 1.2 mm. thick saw. To obtain a crack tip, the last 5 mm of the cut was made with a thin blade 0.5 mm thick.

Fig. 2. DCB sample with crack tip

For a sample with a total height h, measurements were taken on the R-T plane of angles γ_i , i = 1,2,3, of the arrangement of early and late wood layers corresponding to heights h₁ = 0.25h, h₂ = 0.5h and h₃ = 0.75h (Fig. 3).

Next, the average value of angle was determined:

Fig. 3. Diagram for determining the angle of the sample layers

Experimental tests were performed using an MTS Insight testing machine with a non-invasive ARAMIS 3D 4M measurement system (Fig. 4).The DCB samples were installed in the machine using pins and special mounting brackets. The machine was programmed to apply a constant displacement of 5 mm/min with frequency of 10 Hz.

Fig. 4. MTS Insight measurement station with a DCB sample attached, together with the ARAMIS 3D 4M measurement system

The following subsections describe how to determine the strain energy release rate, the impact and selection of input parameters, and how to determine the opening length of the crack tip using Digital Image Correlation (DIC).

Strain energy release rate

The strain energy release rate G for a crack of length a is defined as follows [12]:

where the surface energy increase δE and the area δA the newly formed propagation surface was determined using the following formulas:

The increase in external forces δW corresponding to the increase in crack length δa in a sample with thicknessband the change in strain energy δU in formulas (3) and (4) were determined as follows:

where P(a_i) and δ(a_i) are the force and displacement values of the machine for two crack lengths a_i , i = 1,2 (Fig. 5). In this figure, the blue line in the P-δ diagram shows the curve obtained from a tensile test using a testing machine for a fixed value of crack length increase δa.

Fig. 5. Determination of the Determination of the strain energy release rate

Taking into account the relationship between equations (3)-(6) and equation (2), we obtain:

where G(a)bδa is equal to the area of the shaded region in Fig. 5.

Compliance function

In order to describe the relationship between the compliance function C(a) on the angle γ (10) DCB samples were modeled in ANSYS APDL using SOLID 186 finite elements (Fig. 6). Concentration elements were created around the top of the crack.

At the upper (I) and lower (II) sample mounting points, the relevant boundary conditions are as follows:

where u_x,u_y,u_z are the components of the displacement vector in the Cartesian coordinate system Oxyz. In order to assign material properties, a local coordinate system Ox^,y^,z was introduced, which was used to simulate the change in the arrangement of early and late wood layers in the actual sample. The origin of the global Oxyz and local Ox^,y^,z coordinate systems is located inside the sample at the center of the crack tip.

Fig. 6. Three-dimensional (3D) model of a DCB sample

The displacement value δ = 5 mm.was examined. The values of the elastic module E , Kirchhoff's modulus G and Poisson's ratio v in individual directions and planes in the Oxyz coordinate system are presented in Table 1.

Tab. 1. Material properties for Pinus sylvestris [13]

Ex	Ey	Ez	Gxy
11.7 GPa	0.94 GPa	0.41 GPa	0.95 GPa
Gxz	Gyz	vxy	vxz
0.78 GPa	0.105 GPa	0.43	0.46
vyz	vyx	vzx	vzy
0.65	0.041	0.018	0.32

The calculations were performed for angles γ ranging from 0° to 90° in increments of 5°.The length of the crack a was increased from 105 mm in 145 mm in 5 mm increments. As a result, a data matrix with dimensions (19×9) was obtained, on the basis of which the following form of the compliance function was determined using the power regression method:

where the angle γ-dependent values of the coefficients b_i, i = 0,1,2,3,4 are given in Table A.1 from Appendix A. The C_num(a)(10) function allows you to obtain points on the P - δ corresponding to the assumed crack increase δa. The compliance function C_num(a) (10) was calibrated by its relationship to the corresponding experimentally obtained value C_exp(a):

The proportionality coefficient α in formula (11) is defined as follows:

where

this is the compliance of the sample at the initial crack length a₀.

An analysis was conducted of the impact of mesh density on the change in the value of the force P read at the point of application of displacement δ, perpendicular to the plane of crack propagation. The reading was performed for various n_el numbers of finite elements converging vertically to the top of the sample crack at a height of h/2 (Fig.3). It has been established that the density of the mesh around the top of the crack does not significantly affect the measured force value. At n_el = 6, 8, 10 the force values were P = 68.045 N,  68.042 N,  68.041 N, respectively, and the total number of finite elements was n^total = 29100,  32950,  36080. respectively. Therefore, to obtain the function C_num(a) (10)-(13), a density of n_el = 8 was adopted.

Critical values of fracture parameters

It has been established that changing the δa parameter affects the characteristics of the relationship between force Pand displacement δ: smaller δa corresponds to more experimental points (Fig. 7).

Fig. 7. Matching the value of force P from displacement δ for three values of crack length increase δa at γ = 41.09°

Fig. 8. Strain energy release rate G(a) (7) for three values of crack length increment δa at γ = 41.09°

Based on the data presented in Figure 7, the determination of the strain energy release rate G(a) (7) was determined for the same three values of δa (Fig.8). It has been observed that a decrease in the assumed crack length increase δa causes an increase in the noise level of the release rate of the strain energy G(a).

In summary, it should be noted that in order to obtain the best possible fit of the force value P (Fig.7) and to achieve the lowest noise level in the relationship between the strain energy release rate G(a) (Fig.8) a mean value of δa = 0.525 mm was adopted in the experimental tests.

For the δa value selected in this way at γ = 41.09° from the P-δ graph, it was determined that the maximum force value is P_max = 72.235 N. This P_max value corresponds to the value of the released energy rate G_Pmax = 0.01898 N/mm (Fig.9). The average value of G_m of the function G(a) was found in the range from the value G_Pmax to the value G_Pl at which the linear decrease in force P (dashed line in Figure 9) is still maintained. The critical value G_crit was determined by finding the intersection of the average value G_m with the first value of the function G(a) (7) that differs the least from G_m (G(a) ≈ G_m). Then, the critical force value P_crit = 70.46 N corresponding to G_crit was read from the P-δ graph. It should be noted that in scientific literature, the critical value of the strain energy release rate G_crit is also known as the fracture energy at the first loading mode (MODE I) G_Ic [14].

Fig. 9. Graph of the G(a) (7) function for δa = 0.525 mm at γ = 41.09°

The critical length l_cz of the fracture process zone before the crack tip is determined using the following formula:

where a_crit is the crack length corresponding to G_crit, a₀ is the initial crack length (Fig.9).

The δ_n^c width of the crack tip opening displacement (CTOD) corresponding to the critical value of the released strain energy rate G_crit was determined using digital image correlation (DIC) of the ARAMIS 3D 4M measurement system. The measuring device was calibrated for a working area of 65mm x 48mm with an calibration deviation of 0.00081 mm. In order to synchronize the data, the sampling frequencies of the measuring system and the testing machine were identical. The CTOD determination scheme is shown in Figure 10. First, a pattern of spots was applied to the measurement area of the sample from which the CTOD was read (Fig. 10a). For the sample prepared in this way, a grid of points was created in the Gom Correlate 2018 software, against which the changes occurring on the surface of the sample are read at each successive time step. On the grid developed in this way, two measurement points 1mm apart, located symmetrically and perpendicular to the propagation plane, were selected at the apex of the incised crack (Fig. 10b). The origin of the Cartesian coordinate system was set at the tip of the crack (Fig. 10c). The width of the crack tip opening was determined using the following formula:

where δ_nig ,δ_nid are displacements perpendicular to the plane of propagation of the upper and lower halves of the sample, respectively (Fig. 10d).

Fig. 10. Schematic diagram of CTOD determination in the Gom Correlate 2018 program: a) creation of a point cloud on the sample; b) introduction of measurement points for measuring vertical displacements of the upper and lower layers of the sample; c) location of the coordinate system relative to the top of the crack; d) vertical displacements at the top of the crack

RESEARCH RESULTS

The results obtained were verified by comparing them with the corresponding data obtained for DCB samples in the case of RL and TL systems in [3,4]. The critical values of the strain energy release rate for the RL and TL systems obtained in [3] are almost identical and amount to G_Ic = 0.228 N/mm (s = 0.058 N/mm) and G_Ic = 0.229 N/mm (s = 0.048 N/mm), where s is the standard deviation (Fig. 11). In contrast, in this work, the corresponding values are G_Ic = 0.220 N/mm (s = 0.033 N/mm) and G_Ic = 0.477 N/mm (s = 0.036 N/mm). The nearly twofold difference in the results obtained for the TL system in this study and in [3] can be explained as follows. In [3], the point at which crack propagation initiates on the P-δ diagram occurs before the maximum stress value P_max is reached. In contrast, in this study, this point occurs after P_max is reached. It should be noted that in [3], 20 trials were conducted for each direction, whereas in this study, the number of trials was 6 for the RL direction and 7 for the TL direction. Hence, the standard deviations in [3] are significantly smaller than those presented in Figure 11.

Using the CBBM method in [4] for RL and TL systems, G_Ic = 0.218 N/mm (s = 0.034 N/mm) and G_Ic = 0.476 N/mm (s = 0.039 N/mm) were obtained, respectively. These values, with an accuracy of 0.1% are consistent with the results of this work presented above.

The main objective of this work was to develop a methodology for determining the fracture parameters of DCB samples not only for the RL and TL systems, but for any angle γ in the range from 0° (RL) to 90°(TL).The critical values of the strain energy release rate G_Ic depending on the parameter γ are shown in Figure 12a. The corresponding results for the width δ_n^c of the crack tip opening are presented in Figure 12b.

Increasing the angle γ does not change the critical length l_cz of the fracture process zone (Fig. 12c). The average values of l_cz obtained for angles γ = 0° and γ = 90° differ slightly and amounted to 14.263 mm and 14.477 mm, respectively.

Based on the experimental data shown in Figure 12, the following approximation relationships were developed:

where R² is the coefficient of determination. The G_Ic(γ), δ_n^c(γ) and l_cz(γ) functions are shown by the red line in Figures 12a, 12b, and 12c, respectively. It should be noted that the R² values in equations (16)–(18) differ. This can be explained by the fact that the critical value of the released deformation strain energy G_Ic is influenced solely by changes in the material properties of the wood, as well as the shape and dimensions of the specimen. In contrast, the values of δ_n^c are influenced by structural defects in the wood, such as radial cells and structural defects. Furthermore, in [4] it was observed that the difference in the value of the parameter l_cz for the RL and TL systems is negligible. Due to the low R² values in equations (17) and (18), these relationships can be treated as an approximate linear approximation.

SUMMARY

A methodology for determining the critical parameters for the fracture of pine wood with a crack inclined at any angle to the annual rings direction has been proposed. It consists of the following steps:

Tests conducted according to this methodology on DCB samples with a moisture content of 12% showed that:

The developed methodology and the research results obtained on its basis can be used to predict the cracking of finished GLULAM structural elements.