Study on the strength of glued laminated timber beams with round holes: proposal of the design formula for the splitting strength

When a glued laminated timber (GLT) beam with a round hole is subjected to a shear force and bending moment, the hole will crack and fail due to a large tensile stress perpendicular to the grain and shear stresses. If the stresses acting on the hole are known, it is possible to estimate the fracture load. However, it is necessary to obtain the stresses acting on the hole by finite element analysis, which is very time consuming. In this study, to easily estimate the fracture load, we proposed a formula to estimate the bearing capacity at the time of a hole fracture by obtaining the stress acting on the hole through finite element analysis and an approximate formula. The validity of the proposed formula was verified using the existing experimental data of a GLT beam. As a result, it was confirmed that the proposed equation can estimate the fracture load of GLT beams in Japan and that the proposed equation can estimate the fracture load of GLT beams in countries other than Japan with some accuracy.


Introduction
When a glued laminated timber (GLT) beam with a round hole is affected by a shear force or bending moment, a large tensile stress perpendicular to the grain acts on the cracking plane, and the tensile stress decreases exponentially with the distance from the hole (Fig. 1). Since the shear force and bending moment act at the same time, a crack occurs in the range of θ = 45° to 60° in Fig. 1, where the tensile stress perpendicular to the grain is the maximum tensile stress. When there is a stress gradient, as shown in Fig. 1, it is difficult to estimate the bearing capacity at the moment of split fracture based on material mechanics. Therefore, a mean stress method (a generalization of linear elastic fracture mechanics, as shown in Fig. 2) applicable to beams with notches or holes has been proposed [1]. However, the mean stress method has a problem associated with overestimating the bearing capacity at fracture when the hole is very small or large. Therefore, the authors proposed Eqs. (1) and (2), which take the size effect into account in the mean stress method. Equations (1) and (2) were validated against GLT beams of Scots pine and Japanese larch specified in the Japanese Agricultural Standard (JAS) [2-4]: (1) max (σ /F t90 ) 2 + (τ /F s ) 2 ≤ k vol , k vol = 30 max(0.2D, 30) 0.14 , (2) x = a ms 1 + a ms where D is the diameter of the hole, F t90 is the tensile strength perpendicular to the grain, F s is the shear strength, σ and τ are the mean values of the tensile stress σ and the shear stress τ perpendicular to the grain in the potential fracture area of length x (Fig. 2), k vol is the factor of the size effect, G Ic is the mode I fracture energy, E x is the modulus of elasticity parallel to the grain, E y is the Young's modulus perpendicular to the grain, G xy is the shear modulus, ν xy is the Poisson's ratio, α is a constant (if the hole is small, α = 1 [3] and for D/H ≤ 0.5, α = 1 because the shape of the stress distribution is the same), and H is the beam height. Since σ is the dominant stress at fracture, a ms is set to the value of pure mode I.
Equations (1) and (2) can be used to estimate the bearing capacity of a hole at the time of fracture, but this is very time consuming because the stresses acting around the hole must be determined by finite element analysis (FEA). Therefore, in this study, the maximum stress acting around the hole and the maximum value of the average stress over a certain length were formulated using FEA. Then, a simple calculation formula was proposed to estimate the bearing capacity at the time of a hole fracture. The validity of the formula was verified using existing experimental data. However, in the design code for laminated wood with round holes, the maximum value of the hole diameter (D)/beam height (H) is 0.5 at most [5], so D/H ≤ 0.5 is selected as the applicable range.

Finite element analysis
The stresses acting on the holes in the beam were determined by two types of two-dimensional finite element analysis (2D-FEA), "analysis of stress around the hole" and "analysis of the stress distribution in the direction parallel to the grain". The stresses acting on the holes in the beam were determined by 2D-FEA using ANSYS 18.2. The specification of the 2D-FEA was a case where only the shear force or bending moment acts on the center of the hole with D/H as a variable. Then, an approximate expression for the average stress for Eq. (1) was derived from the stresses obtained by the 2D-FEA. The details of the 2D-FEA are described below. The numerical finite element (FE) model was set up for a case where the stress at the center of the hole was shear force only (Q= τ 1 BH (N), M = 0, and τ 1 = 1 N/mm 2 ) and bending moment only (Q = 0, M = τ 1 BH 2 (N mm), and τ 1 = 1 N/mm 2 ), as shown in Fig. 3. "Edge and node parallel to the grain" is the edge and node provided to obtain the stress distribution in the direction parallel to the grain. "Edge and node parallel to the grain" is not placed in "analysis of stress around the hole", but "edge and node parallel to the grain" is placed in "analysis of stress distribution in the direction parallel to the grain". "Edges and nodes parallel to the grain" were located every 0.01D in   Fig. 1) perpendicular to the grain (y-direction) from the center of the hole. The D/H values were set to 0.05, 0.1, 0.2, 0.3, 0.4, and 0.5. However, analysis of D/H values of 0.6, 0.7, 0.8, and 0.9 was also performed for reference only in "analysis of stress around the hole". Since the stresses around the hole when the shear force of Q = τ 1 BH (N) or the bending moment of M = τ 1 BH 2 (N mm) acts on the beam are the same even if the beam height H is different, H = 300 mm was used. The material parameters of the elements were determined in the same way as in [3,4]. The results are the same even if the value of the Young's modulus E x in the fiber direction is different, but for the time being, E x = 10,000 N/mm 2 was used. The Young's modulus perpendicular to the grain was set to E y = E x /25, the shear modulus was set to G xy = E x /15, and the Poisson's ratio was set to ν xy = 0.4. The wood was treated as a four-node planar element, and the elements were less than 5 mm (H/60). The range of 10 mm (H/30) from the hole was divided into small sections so that the nodes were placed at less than one degree. The mesh size of the FE model for "analysis of stress around the hole" is shown in Additional file 1: Fig. S1. The mesh size of the FE model for "analysis of the stress distribution in the direction parallel to the grain" is shown in Additional file 1: Fig. S2. The elements were very finely divided around only the hole, so some 3-node planar elements were placed. For the numerical FE results, when all the elements were reduced in size and the number of elements was doubled, the error in the stress around the holes was less than 1%, so we concluded that the optimal solution could be obtained by the element division shown in Additional file 1: Figs. S1 and S2.

Experimental data for beams with round holes
The existing experimental data [3,4,[6][7][8][9][10][11][12][13] of a GLT beam with only one round hole with D/H ≤ 0.5 in the center were used to verify the validity of a simple calculation formula to estimate the bearing capacity at the time  of a hole fracture. Table 1 shows the experimental data for GLT beams in Japan [3,4,6,7], and Table 2 shows the experimental data for GLT beams in countries other than Japan [8][9][10][11][12][13]. Table 2 uses the values summarized in detail by Danielsson [14].

Material properties
E x , F s , F t90 , and G Ic used in Eqs. (1) and (2) were obtained from the density of the wood in the area where a hole was provided. E y , G xy , and ν xy were the same as those in the 2D-FEA, with E y = E x /25, G xy = E x /15, and ν xy = 0.4. The density of the wood in the hole area was set to 522 kg/ m 3 for OKAc [4], 458 kg/m 3 for PEN [9], 489 kg/m 3 for HOF [12], and 471 kg/m 3 for AIC [13], as shown in Tables 1 and 2. OKAa [3], OKAb [3], KAR [6], and HIJ [7] are GLT beams composed of heterogeneous grades, and the density of the wood in the area where a hole is provided is not known. Therefore, the density was estimated from the grade of the inner layer lamina (D/H = 0.3 at H = 150 mm for OKAa [3] and OKAb [3] is the middle layer lamina) in the area where the hole is provided. Figure 4 shows the relationship between the grade L and M is the bending moment at the center of the hole, Q is the shear force at the center of the hole, Q Ini. is the shear force at the center of the hole when the initial crack has propagated over the entire beam width, Q Max. is the shear force at the center of the hole at failure, E105-F300 and E120-F330 are GLT beams composed of heterogeneous glue, and ■represents the yield strength evaluated by the elastic-perfectly plastic model  M is the bending moment at the center of the hole, Q is the shear force at the center of the hole, Q Ini. is the shear force at the center of the hole when the initial crack has propagated over the entire beam width, Q Max. is the shear force at the center of the hole at failure, and ■ represents the average value of 5 test specimens density ρ of the specimens extracted from the existing experimental lamina data [15]. The density of the lamina was taken as obtained from the least-squares approximation line in Fig. 4, which is 446 kg/m 3 for the inner lamina of OKAa [3] and OKAb [3] with L80, 482 kg/m 3 for the middle layer lamina of OKAa [3] and OKAb [3] with L100, and 464 kg/m 3 for the inner lamina of KAR [6] and HIJ [7] with L90. BEN [8] and JOH [10] are assumed to be equivalent to PEN [9] of the same grade, so 458 kg/m 3 is used. HAL [11]  (3) E x N/mm 2 = 23.536 × ρ kg/m 3 .  3) with existing intensity data [17,18] and previous experimental data [3]. Since the lamina grade shows the lower limit, using the least-squares approximation line in Fig. 4 and Eq. 3 results in a value higher than the lamina grade.
F s was obtained from the density ρ of wood using Eq. (4) [19]: Figure 6 shows a comparison of Eq. (4) with existing intensity data [17,18] and previous experimental data [3,4]. F t90 was obtained by multiplying F s by a constant. According to the "Standard for Structural Design of Timber Structures" (Japan) [20], F s was multiplied by 1/3 to obtain F t90 , so the constant was set to 1/3. However, F t90 of Japanese larch is low [3], and a load of splitting fracture perpendicular to the grain against the density has also been confirmed to be lower in Japanese larch than in spruce, Japanese cedar, and Douglas fir [21]. Therefore, the constant was determined separately for Japanese larch. Figure 7 shows a comparison of the calculated values of F t90 with existing intensity data [17] and experimental data from the authors' previous studies [3,4,22]. Figure 8 shows a schematic diagram of the test method used to determine F t90 in the authors' previous studies [3,4,22]. In the studies of beams with round holes [3,4], a specimen height of h = 150 mm was used, but in the study of the tension perpendicular to the grain strength of the GLT beams of Scots pine [22], a specimen height of h = 100 mm was found to be optimal. Although the optimum value of specimen height h (4) F s N/mm 2 = 0.02256 × ρ kg/m 3 − 2.97. exp. Scots pine (Tensile test) [3] exp. Japanese larch (Tensile test) [3] ref. [ Fig. 7, the constant multiplied by F s for Japanese larch was set to 1/4. Equation (5) was used to obtain G Ic from the density ρ of wood [23]: Figure 9 shows a comparison of the results from Eq. (5) with the experimental data from the authors' previous studies [3,4,24]. Since Japanese larch has low G Ic and F t90 values, we decided to use the value obtained by multiplying Eq. (5) by 3/4 as the same reduction for F t90 .

FE results (1) Analysis around the hole
The tensile stresses perpendicular to the grain σ Q and σ M and the shear stresses τ Q and τ M acting around the hole obtained from the analysis in Fig. 3 and Additional file 1: Fig. S1 are shown in Fig. 10. The stresses when only shear force Q acts on the center of the hole are σ Q and τ Q , and the stress levels when only bending moment M acts on the hole are σ M and τ M . The vertical axis represents the stress/ maximum stress, with the maximum value normalized to 1.0. The horizontal axis represents y/D. The stresses acting on the hole are symmetric or inversely symmetric with respect to the center of the hole [4], so only a quarter of the hole is shown in Fig. 10. From Fig. 10, it was found that the (5) G Ic N · m/m 2 = 1.07 × ρ kg/m 3 − 162.
stress distribution is different from D/H of 0.6 or more but almost the same when D/H is 0.5 or less. This indicates that there is a possibility to approximate the stress distribution with a simple equation when D/H is 0.5 or less.
Therefore, the coefficients α Q and α M for finding Q = α Q τ 1 BH and M = α M τ 1 BH 2 acting when the maximum values of σ Q and σ M are almost 1.0 N/mm 2 were derived from σ Q, max and σ M, max in Fig. 10. α Q is shown in Eq. (6), and α Q is shown in Eq. (7). Equations (6) and (7) were obtained by trial and error to make them as simple as possible: σ Q and τ Q when Q = α Q τ 1 BH and σ M and τ M when M = α M τ 1 BH 2 are shown in Fig. 11. σ Q, max is at y ≈ 0.38D and σ M, max is at y ≈ 0.25D. Additional file 1: Table S1 shows the coefficients obtained by regressing the stress distribution of σ Q , σ M , τ Q and τ M in the range of 0.2D ≤ y ≤ 0.4D in Fig. 11 with the 5th-order polynomial in Eq. (8): The effect of σ on the splitting fracture of the hole is very large, while the effect of τ is very small (since k τ = 1.015-1.075 in Eq. (17), the effect of shear stress is approximately 1.5-7.5%). Therefore, it is assumed that the splitting fracture of the hole is in the range of 0.25 ≤ y/D ≤ 0.40 (30° ≤ θ ≤ 53° in Fig. 1) to derive the bearing capacity formula.

(1) Maximum value of σ around the hole
Since the maximum value of σ Q when Q = α Q τ 1 BH (τ 1 = 1 N/mm 2 ) is 1 N/mm 2 and the maximum value of σ M when M = α M τ 1 BH 2 (τ 1 = 1 N/mm 2 ) is 1 N/mm 2 , the maximum values of σ Q and σ M can be obtained using Eqs. (12) and (13): is not true. Therefore, the FE results in Eq. (9) and Additional file 1: Table S2 can be used to obtain the relationship between M/ (QH) and σ max by determining the relationship between M/ (QH) and σ for y = 0.25D to 0.40D, as shown in Additional file 1: Fig. S4 (a case of D/H = 0.3). Then, an approximation exp. Scots pine [3][4] exp. Japanese larch [3] exp. Japanese cedar [24] Eq. (5) [23] 3/4   14)) using Eqs. (12,13) was created by trial and error to make the formula as simple as possible: A comparison between the FE results for σ max and the calculation results of Eq. (14) is shown in Fig. 13. It was confirmed that the results can be estimated with good accuracy in the range of M/(QH) = 1-100.

(3) Maximum value of σ incorporating the effect of shear stress around the hole
The effect of the shear stress is determined by Eq. (16), which is the reformulation of Eq. (11) obtained using the FE results in Eq. (10) and Additional file 1: Tables S2 and  S3: The FE results in Eq. (10) and Additional file 1: Tables S2 and S3 were used to obtain the relationship between x/D and k τ · σ for y = 0.25D to 0.40D, as shown in Additional file 1: Fig. S7 (a case of D/H = 0.3, F s /F t90 = 3.0, and M/(QH) = 0.1, 5, 10, 100). We obtained the relationship between x/D and k τ · σ max for each M/(QH). Then, an approximation formula for k τ (Eq. (17)) was created by trial and error to make the formula as simple as possible. However, since we judged from Fig. 9 that F s /F t90 = 3 to 4, we only obtained k τ for the cases of