Mode I Critical Stress Intensity Factor of Beech Wood ( Fagus Sylvatica ) in a TL Configuration : A Comparison of Different Methods

The paper presents a comparison between various methods of mode I critical stress intensity factor KIC calculations of beech wood in the TL configuration. The first method is the stress intensity factor extrapolation to the distance of 0 mm from the crack tip; the second method is the use of the J integral; and the third method is based on the differences in deformation energies from which the strain energy release rate per unit of crack propagation length was obtained. The fourth method is the calculation of material deformation around the crack or the displacement of the triangle element node; and the fifth method uses a generally known equation for the CT specimen for plane-strain conditions in isotropic material. Using the finite element method, it was found that the J integral was least sensitive to the size and shape of the elements. It was used to calculate the critical stress intensity factor KIC for beech wood in a TL configuration. The average value is 0.56 MPa√m with a standard deviation of 0.047 MPa√m.


INTRODUCTION
1. UVOD Wood fracture has already been studied by several researchers.Porter (1964) measured critical strain energy release rate G IC by measuring the force and length of crack in white pine.The specimens were TL and RL-oriented, which means that he loaded them in tangential and radial directions, respectively, and the crack propagated in a longitudinal direction.He researched the influence of length, thickness, and height of the specimen, as well as the crack length.The method of determining the critical strain energy release rate G IC in pines was also studied by Stanzl-Tschegg et al. (1995) .They calculated the strain energy release rate G I by integrating the energy or the area under the curve, which describes the force depending on the specimen mouth opening.The energy obtained was divided by the size of the newly formed surface.Thuvander and Berglund (2000) researched pine fracture in the TR orientation.They stated that K IC of silver fir wood in the TR orientation was between 30 % and 50 % higher than that in the TL direction, while in the case of pine and spruce the difference was supposed to be even greater.As the newly formed surface proved to be equal in the TR and TL orientations, i.e., in the RL plane, they wanted to know the reason for such great differences in K IC .Fresh specimens were used because in dry specimens they encountered the problem of stable fracture due to microcracks resulting from the drying process.According to them, a specimen that has been dried and humidified again has a more brittle fracture due to microcracks formed during the process of drying.Similar tests in the TR orientation were also performed by Frühman et.al. (2003).
The influence of the moisture content of wood on K IC was also researched by many researchers (Ozyhar et. al. (2012), Reiterer and Tschegg (2002), Scheffler et.al. (2004), Vasić and Stanzl-Tschegg (2007), Yeh and Schniewind (1992)).They found that with increasing moisture content, the critical stress intensity factor in the RL and TR orientation decreases.Vasić and Stanzl-Tschegg (2007) report a value of 0.9 MPa•√m with 6 % moisture content of beech wood in the RL configuration, and the value of 0.62 MPa•√m with 12 % moisture content.
In their work, Stanzl-Tschegg and Navi (2009) sum up their research of wood fracture in the RL configuration under various conditions such as moisture and density of wood, combined fracture modes I and II, and loading rate.They mention the work of Beikircher who thermally modified wood and found that thermal modification of beech decreases K IC for the TL configuration from 0.8 MPa•√m to 0.6 MPa•√m.Likewise, Stanzl-Tschegg and Navi (2009) state that K IC in the RL orientation is higher compared to the TL orientation because of the 'bridging' effect that the parenchyma causes in a radial direction.Majano-Majano et.al.(2012) found that K IC of thermally modified beech wood in RL and TL configuration decreases and found that the K IC for unmodified beech in TL configuration form ranged from 0.44 to 0.63 MPa•√m.
To determine the fracture properties of wood, the majority of the aforementioned authors use the critical strain energy release rate G IC , determining the energy necessary for the formation of new surfaces on the basis of the force and mouth opening.A prerequisite for an experiment of this kind is a stable advance of the fracture, which means that the crack propagates in proportion to the crack mouth opening.However, the problem in the case of beech wood in the TL configuration is that after initiation the fracture process is distinctly unstable, during which the crack suddenly propagates to a certain unbalanced length.G IC does in fact express the energy to be put in per unit of the newly formed surface, but provides little information about the fracture initiation, which is of essential significance in cutting.Under certain conditions, the result can be a chip of type I, II, or III, as classified by Franz (Koch, 1985).Regarding the material which is turned into a chip, a type I chip is discontinued, formed by alternating fracture and bending failures (Merhar and Bučar, 2012).Whether the fracture under the chip will progress or the chip will break depends on the critical stress intensity factor for fracture mode I.It was, therefore, decided to determine the K IC of beech wood for the TL configuration in the manner enabling a direct determination of the value.First, the K IC was going to be determined on one specimen using the five most frequently used methods.The results obtained were going to be used to determine the most suitable method that yielded a satisfactory result in a simple manner.The method obtained in this way would be used to determine the fracture toughness of the remaining specimens.

MATERIJALI I METODE
Beech wood (Fagus sylvatica) specimens were taken from a peripheral part of one stem of 400 mm in diameter.The specimens moisture content was (9 ± 0.5) %, with a density of 630 kg/m 3 .A conventional compact tension CT specimen (Hertzberg, 1996) of 115 mm in length, 100 mm in height, and 10 mm thick was made.
The specimens were TL-oriented, which means that load was applied in a tangential direction, and the crack propagated longitudinally.Since the board used to make specimens was only 80 mm thick, an additional 10 mm thick strip of wood was glued to each side of the specimen.The glued-on strip was obtained from the immediate vicinity of the specimen so that it had similar mechanical properties.The obtained CT specimen was modified to enable the mounting of the crack mouth opening displacement meter as shown in Figure 1a.
A cut of approximately 56 mm in length was made in the specimen (Figure 1b).Then a razor blade was used to make a further 1 mm to 2 mm deep cut to obtain a sharp tip of the cut.After fracturing, the initial crack length was measured on each specimen.
The crack mouth opening displacement meter was made of two 1 mm thick and 90 mm long spring steel gauges.The length of the gauges was determined so as to make the meter measuring range of 10 mm, whereby the stress inside the gauge as a result of bending did not exceed 200-300 MPa.In this case each gauge was deformed by 5 mm.Strain gauges type 3/120LG11 produced by HBM, with the resistance of 120 W, were glued on the upper and lower side of both steel gauges, at the beginning or at the place of the maximum bending moment, and connected to the MES HPSC 3102 amplifier with full-bridge configuration.The displacement meter was calibrated by means of a reference dial gauge with the accuracy of 0.01 mm.
The specimen was placed on the tensile testing machine where the loading force was measured by a dynamometer, and the crack mouth opening displacement was measured by the previously described displacement meter.Data were captured by means of a personal computer, NI PCI-6014 measurement card and LabView software by National Instruments.The data capturing rate was 200 Hz, and the mouth opening velocity 15 mm/min.
The intersection of the measured data and a straight line with a 5 % smaller angle than the straight line representing a linear regression curve of the initial elastic part of the specimen loading, as laid down by the ASTM E 399 standard, was read from the measured data.At the same time the maximum force measured was read.When the maximum force was greater than the intersection of the measured data and straight line, or the difference between the value at the intersection and the maximum measured value was minimal, the maximum force measured was taken into account for the calculation.
The specimen was modelled by the finite element method using the Ansys program.The orthotropic properties of the wood were taken into consideration.
The measured modulus of elasticity in a longitudinal direction, which has the greater influence on the K IC calculation, was used in the model, while the data for moduli of elasticity in other directions, shear moduli and Poisson's ratios, which have minor influence on the K IC calculation, were taken from Kollmann and Cote (1984).Thus: Since the aim of the first part of the research was to investigate the accuracy of determining the mode I critical stress intensity factor K IC , the specimen was modelled with a linear elastic plane-strain state, where the modelled specimen thickness was 10.A PLANE183 higher-order 2D, 8-node brick element was used.The specimen had elements of 2mm and 1mm in size, with a combination of finer elements around the crack tip or special triangle elements for calculating the stress at the crack tip with an intermediate node at ¼ of the element's length.The tip of the crack was surrounded by two rows of 12 triangle elements each, whereby the elements length was 0.1 mm or 1/1000 of the crack length.The ratio between the size of the first row and the second one was set as 1.5.The mode I critical stress intensity factor K IC was calculated by five different and most frequently used methods of determining the critical stress intensity factor.The results of the comparison of these five methods were used to determine the most accurate and simple method, and this method was subsequently used to calculate the critical stress intensity factor for the remaining specimens.
The first method used to calculate the critical stress intensity factor K IC was the stress intensity factor extrapolation to the distance of zero using equation 2 (Broek, 1989) where r is the distance from the crack tip, and s y is stress in y direction as shown in Figure 2, and it is expressed as follows Since the specimen was TL-oriented, this means that it was loaded in a tangential direction, but due to tissue orientation the crack propagated in a longitudinal direction, the angle q in equation 3 equals 0. Only stresses in nodes lying in the crack propagation plane were thus taken into account in the calculation.The second method of the critical stress intensity factor calculation used the J integral according to equation 4 (Broek, 1989) and Figure 3 (4) G is integration path, and W d is deformation energy per unit of volume, T is stress vector acting perpendicularly on contour Γ, differential. (5) T is stress vector acting perpendicularly on contour G, u is deformation vector, and ds is the G path differential.The value of the J integral was calculated by means of the Ansys programme so that its contour or G contour nodes were defined.
Since the specimen was modelled as linearly elastic, the value of J integral can be equalled with the elastic energy release rate G (Smith, 2003) The following equation is also taken into account where E` is an equivalent modulus of elasticity (Sih et al., 1965) The value of the J integral was calculated by means of the Ansys Γ contour nodes were defined.
Since the specimen was modelled as linearly elastic, the va with the elastic energy release rate G (Smith, 2003) The following equation is also taken into account and b ij are compliance constants depending on the type of material.In the case of the plane-strain condition, the b ij constants must be calculated from coefficients a ij where The third method of the critical stress intensity factor calculation u strain energies, from which the strain energy release rate per unit of crack p (10) The compliance coefficient values for the value of a 11 were calculated from the measured modulus of elasticity, which amounted to 14 490 MPa, while the values Since the specimen was TL-oriented, this means that it was loaded in a tangential direction, but due to tissue orientation the crack propagated in a longitudinal direction, the angle θ in equation 3 equals 0. Only stresses in nodes lying in the crack propagation plane were thus taken into account in the calculation.The second method of the critical stress intensity factor calculation used the J integral according to equation 4 (Broek, 1989) and Figure 3  Γ is integration path, and W d is deformation energy per unit of volume, The third method of the critical stress intensity factor calculation used the differences of strain energies, from which the strain energy release rate per unit of crack propagation length was obtained.First a crack with a measured length was modelled, and then another with a longer crack.Differences in the lengths of modelled cracks ranged from 0.05mm to 0.35mm.The deformation energy of each modelled crack length was calculated by means of the programme.After that, the difference in the dW d energies between the specimen with the longer crack and the specimen with the measured crack was calculated.The difference was divided by the difference in lengths da and the thickness of the modelled specimen b.
The compliance coefficient values for the value of a 11 were calculated from the measured modulus of elasticity, which amounted to 14490 MPa, while the values for other coefficients were taken from Kollmann and Cote (1984) (equation 1).
The third method of the critical stress intensity factor calculation used the differences of strain energies, from which the strain energy release rate per unit of crack propagation length was obtained.First a crack with a measured length was modelled, and then another with a longer crack.
Differences in the lengths of modelled cracks ranged from 0.05mm to 0.35mm.The deformation energy of each modelled crack length was calculated by means of the programme.After that, the difference in the dW d energies between the specimen with the longer crack and the specimen with the measured crack was calculated.The difference was divided by the difference in lengths da and the thickness of the modelled specimen b.
Then the critical stress intensity factor was calculated considering equations 6 to 10.
The fourth method of critical stress intensity factor calculation was based on the deformation of material surrounding the crack, i.e., by the displacement of triangle element nodes as shown in Figure 4, and using the following equation (Sauoma and Sikiotis, 1986).
[ ] [ ] The matrix [ ] Then the critical stress intensity factor was calculated considering equations 6 to 10.
The fourth method of critical stress intensity factor calculation was based on the deformation of material surrounding the crack, i.e., by the displacement of triangle element nodes as shown in Figure 4, and using the following equation (Sauoma and Sikiotis, 1986).
[ ] [ ] The matrix [ ] The compliance coefficient values for the value of a 11 were calculated from the measured modulus of elasticity, which amounted to 14490 MPa, while the values for other coefficients were taken from Kollmann and Cote (1984) (equation 1).
The third method of the critical stress intensity factor calculation used the differences of strain energies, from which the strain energy release rate per unit of crack propagation length was obtained.First a crack with a measured length was modelled, and then another with a longer crack.
Differences in the lengths of modelled cracks ranged from 0.05mm to 0.35mm.The deformation energy of each modelled crack length was calculated by means of the programme.After that, the difference in the dW d energies between the specimen with the longer crack and the specimen with the measured crack was calculated.The difference was divided by the difference in lengths da and the thickness of the modelled specimen b.
Then the critical stress intensity factor was calculated considering equations 6 to 10.
The fourth method of critical stress intensity factor calculation was based on the deformation of material surrounding the crack, i.e., by the displacement of triangle element nodes as shown in Figure 4, and using the following equation (Sauoma and Sikiotis, 1986).
[ ] [ ] The matrix [ ] [ ] and matrix [ ] and matrix [ ] u B and u C are the displacements of nodes B and C in the x direction, and v B and v C are the displacements of nodes B and C in the y direction.L 1 is the length of triangle element, and s 1 and s 2 are complex zeros of the equation  In the case of plane-stress condition, the a ij constants are equation 10, and in the case of plane-strain condition by equati (17) In the case of plane-stress condition, the a ij constants are compliance coefficients yielded by equation 10, and in the case of plane-strain condition by equation 9. Since the calculation takes account of the material deformation on only one side of the crack, the method can be used only for symmetric specimens with a symmetric load application.In order for the calculation to be as accurate as possible, only half of the specimen was modelled and on the lower side the programme was set a boundary condition that the specimen was symmetric.The crack tip was surrounded by two rows of 6 triangle elements each, the elements being around 0. ond one was set as 1.5.The results were used to calculate the deformation of nodes in directions x and y, and to calculate the critical stress intensity factor K IC for fracture mode I.
The fifth method of the critical stress intensity factor calculation used a generally known equation (Broek, 1989) applying to the CT specimen for planestrain conditions in isotropic material on the lower side the programme was set a boundary condition that the specimen was symmetric.
The crack tip was surrounded by two rows of 6 triangle elements each, the elements being around 0.1 mm or 1/1000 of the crack length long.The ratio between the size of the first row and the second one was set as 1.5.The results were used to calculate the deformation of nodes in directions x and y, and to calculate the critical stress intensity factor K IC for fracture mode I.
The fifth method of the critical stress intensity factor calculation used a generally known equation (Broek, 1989) applying to the CT specimen for plane-strain conditions in isotropic material where a is the crack length, W is specimen length and B specimen thickness.
The specimens moduli of elasticity E L in a longitudinal direction were also measured.
Specimens -130 mm long, 10 mm wide, and 6 mm high -were subjected to a four-point bending load on the tensile testing machine.A linear variable differential transformer (LVDT) was used to measure the specimen displacement during loading.
From the graph of the measurements of force depending on deformation, the linear regression curve coefficient in the area of linear dependence between force and deformation was determined by means of the Excel programme.The equation describing the displacement in the middle of the specimen depending on the specimen geometric data and loading force was used to calculate the modulus of elasticity E L .

REZULTATI I RASPRAVA
on the lower side the programme was set a boundary condition that the specimen was symmetric.
The crack tip was surrounded by two rows of 6 triangle elements each, the elements being around 0.1 mm or 1/1000 of the crack length long.The ratio between the size of the first row and the second one was set as 1.5.The results were used to calculate the deformation of nodes in directions x and y, and to calculate the critical stress intensity factor K IC for fracture mode I.
The fifth method of the critical stress intensity factor calculation used a generally known equation (Broek, 1989) applying to the CT specimen for plane-strain conditions in isotropic material where a is the crack length, W is specimen length and B specimen thickness.
The specimens moduli of elasticity E L in a longitudinal direction were also measured.Specimens -130 mm long, 10 mm wide, and 6 mm high -were subjected to a four-point bending load on the tensile testing machine.A linear variable differential transformer (LVDT) was used to measure the specimen displacement during loading.
From the graph of the measurements of force depending on deformation, the linear regression curve coefficient in the area of linear dependence between force and deformation was determined by means of the Excel programme.The equation describing the displacement in the middle of the specimen depending on the specimen geometric data and loading force was used to calculate the modulus of elasticity E L .

REZULTATI I RASPRAVA (18)
Figure 5 shows forces dependent on the crack mouth opening displacement of the specimen used to determine the critical stress intensity factor.The figure shows instantaneous force drop as a consequence of sudden crack propagation or an instable fracture.The figure clearly shows the linear elastic part of force dependence on mouth opening, shifting to the nonlinear part just before the crack propagates.The continuous where a is the crack length, W is specimen length and B specimen thickness.
The specimens moduli of elasticity E L in a longitudinal direction were also measured.Specimens -130 mm long, 10 mm wide, and 6 mm high -were subjected to a four-point bending load on the tensile testing machine.A linear variable differential transformer (LVDT) was used to measure the specimen displacement during loading.
From the graph of the measurements of force depending on deformation, the linear regression curve coefficient in the area of linear dependence between force and deformation was determined by means of the Excel programme.The equation describing the displacement in the middle of the specimen depending on the specimen geometric data and loading force was used to calculate the modulus of elasticity E L .

REZULTATI I RASPRAVA
Table 1 indicates the measured values of the modulus of elasticity.The table shows that the specimens have an average modulus of elasticity E L of 14 487 MPa, with standard deviation of 1 246 MPa.line is a regression curve for the linear part of loading, while the inclination of the dashed line is by 5 % smaller than the continuous one.Since the intersection of a straight line with a 5 % smaller inclination and the measurements were practically equal to the maximum force, the maximum forces measured were used in the calculations.Table 2 indicates the results of the critical stress intensity factor calculation for five different methods.The results of calculations written in bold differ insignificantly from each other.Calculations using the J integral (Figure 6) have been shown as the least sensitive to the size of elements and to the range of integration.The calculation values were around 0.496 MPa•√m, regardless of the size of basic elements, the type of elements surrounding the crack, and the distance of contour around the crack tip, up to the distance of 0.4 mm.At this distance, the integration path included two types of triangle elements and at least one type of 8-node brick element.In the case of a shorter distance, however, it was demonstrated that a satisfactory result requires at least two types of 8-node brick elements, as is the case with 1 mm large elements around the crack tip.In the case of two rows of triangle elements the accuracy of result is not satisfactory.In this case, at least one more row of 8-node brick elements is required.Calculations based on difference of elastic deformation energy exhibit greater sensitivity to the size of elements as well as to the type of elements around the crack tip.The results are better in the case of triangle elements surrounding the tip.In the case of smaller elements, there is no deviation from the calculation using J integral and other two methods regardless of the size of the crack extension.When the crack is extended by 0.05 mm or 0.1% of the crack size, the results are no longer satisfactory whatever the size of the basic elements and type of elements surrounding the crack tip.In the case of 8-node brick elements where the refinement is made, the results are not satisfactory regardless of the size of elements.When elements of 2 mm in size were used and the crack was extended by 0.1mm, no calculation could be made because the elastic deformation energy decreased with an increased crack, which is contrary to other cases.
The calculation based on the stress intensity factor extrapolation at various distances of the crack tip to the distance of zero yielded equal results regardless of the size of the basic elements.The tip was surrounded only by elements with 8 nodes, since during the refinement of elements in the crack propagation line the triangle elements turned into 8-node brick elements.Stresses in the y direction and stress intensity factors as a consequence of the crack tip distance together with regression curve are shown in Figure 7.It clearly shows the linear dependence of the stress intensity factor on distance.
Likewise, the calculation based on the deformation of triangle element (Eq.12) nodes was equal for both sizes of basic elements, while the triangle element size was the same in both cases, equalling 0.118 mm.
The calculation of critical stress intensity factor by a generally known equation for an isotropic CT specimen (Eq.18) yielded about 5 % higher results.In our case, the reason for that was probably the orthotropic material, which has significantly lower modulus of elasticity in radial and in tangential direction compared to the longitudinal one, but in the equation 12 only the modulus of elasticity in longitudinal direction is used.Another reason could also be the specimen height-to-length ratio.This is to say that the equation assumes the mentioned ratio to be 1, while in our case it was less than 1.
Since the calculation of the critical stress intensity factor by means of the J integral, in which the contour is sufficiently far from the tip, is satisfactory, and the calculation method simple, the J integral was used to calculate K IC for the remaining specimens.
Table 3 shows the critical stress intensity factors calculated by means of the J integral.The integral contour was 6 mm away from the crack tip at the basic elements size of 2 mm and triangle elements around the crack tip.An average value of critical stress intensity factor is 0.56 MPa•√m with standard deviation of 0.047 MPa•√m, which represents less than 10 % of the mean value determined.A similar deviation can be found with the modulus of elasticity.Considering Vasić and Stanzl-Tschegg (2007), who stated the value of 0.62 MPa•√m for the RL orientation as the beech fracture toughness with a moisture content of 12 %, and the fact that fracture toughness in the TL configuration is lower than that of the RL configuration, as stated by Stanzl-Tschegg and Navi (2009), it can be ascertained that the values we obtained comply with their findings as well with Majano-Majano et.al. (2012).It should also be noted that the bigger specimens would probably give more representative K IC values, since they are dimension dependent as reported by Stanzl-Tschegg et.al. (1995).

ZAKLJUČAK
Comparing different methods for mode I critical stress intensity factor calculations, the J integral proves to be the most appropriate considering the simplicity and sensitivity to the size and shape of the elements.It was used to calculate the critical stress intensity factor for beech wood in a TL configuration, which means that the specimens were loaded in a tangential direction while the crack propagated in a longitudinal direction.The average value was 0.56 MPa•√m with a standard deviation of 0.047 MPa•√m.Comparing the results of Vasic and Stanzl-Tschegg (2007), who obtained K IC of 0.62 MPa•√m in the RL configuration at 12 % wood moisture content, and the fact that the value in TL configuration is lower than the value in RL configuration, we find that the values obtained agree with research results of other authors like Majano-Majano et al. (2012), who stated the value for K IC in TL configuration in the range from 0.44 to 0.63 MPa√m and Ozyhar et. al. (2012), who determined K IC in TL configuration to be around 0.406.Likewise, it was found that researchers predominantly investigate the so-called 'stable fracture', which means that the crack propagates in proportion to the crack mouth opening displacement.Beech wood in the TL configuration, however, has shown to be a distinctly brittle material, since after the maximum load is achieved, the crack propagates instantaneously and at a great speed to a certain equilibrium length (Merhar and Bučar, 2012), which various authors consider to be an unstable fracture.However, the RL configuration is exhibited as a more stable one, and therefore several authors prefer to use it as a model for determining the critical stress intensity factor.

Figure 3
Figure 3 J integral Slika 3. J integral is an equivalent modulus of elasticity (Sih et al., b ij are compliance constants depending on the type of materi condition, the b ij constants must be calculated from coefficients a ij

Figure 3 J
Figure 3 J integral Slika 3. J integral

.
The compliance coefficient values for the value of a 11 were calculated from of elasticity, which amounted to 14490 MPa, while the values for other c from Kollmann and Cote (1984) (equation 1).

u
B and u C are the displacements of nodes B and C in the x direction, displacements of nodes B and C in the y direction.L 1 is the length of triang are complex zeros of the equation 1,2).A zero with positive imaginary part for the solution,while p j and q j of plane-stress condition, the a ij constants are compliance equation 10, and in the case of plane-strain condition by equation 9.
the solution,while p j and q j are u B and u C are the displacements of nodes B and C in the displacements of nodes B and C in the y direction.L 1 is the len are complex zeros of the equation 1,2).A zero with positive im for the solution,while p j and q j

Figure 4
Figure 4 Triangle element with nodes Slika 4. Elementi trokuta s čvorovima 1 mm or 1/1000 of the crack length long.The ratio between the size of the first row and the sec-DRVNA INDUSTRIJA 64 (3) 221-229 (2013)

Figure 7 5 Distance
Figure 7Stresses in y direction s y and stress intensity factor K I with regression line as a consequence of distance from the crack tip Slika 7. Prikaz ovisnosti naprezanja u y smjeru s y i faktora intenziteta naprezanja K I o udaljenosti od vrha pukotine regresijskim krivuljama

Table 1
Measured modulus of elasticity in longitudinal direction E L Tablica 1. Izmjereni modul elastičnosti u longitudinalnom smjeru E L