Experimental study on glued laminated timber beams with well-known knot morphology

Abstract

Nowadays, the impact of knots on the failure behaviour of glued laminated timber (GLT) beams is considered by subjecting the single lamellas to a strength grading process, where, i.a., tracheid effect-based laser scanning is used to obtain information about knot properties. This approach single-handedly defines the beam’s final strength properties according to current standards. At the same time, advanced production processes of such beams would allow an easy tracking of a scanned board’s location, but, at this point, previously obtained detailed information is already disregarded. Therefore, the scanning data is used to virtually reconstruct knot geometries and group them into sections within GLT beams. For this study, a sample of 50 GLT beams of five different configuration types was produced and tested under static four-point-bending until failure. As for each assembled lamella the orientation and position within the corresponding GLT beam is known, several parameters derived from the reconstructed knots can be correlated to effective GLT properties. Furthermore, the crack patterns of the tested beams are manually recorded and used to obtain measures of cracks. A detailed analysis of the generated data and their statistical evaluation show that, in the future, dedicated mechanical models for such timber elements must be developed to realistically predict their strength properties. A potential approach, using fluctuating section-wise effective material properties, is proposed.

1 Introduction

Wood is a naturally grown material and therefore is gaining popularity for its sustainability as well as its ecological and often also aesthetic advantages over other construction materials. However, the natural growing process also leads to fluctuations in both morphological features and mechanical properties. The morphology is dominated by the fibre course, which in turn is strongly influenced by knots (Phillips et al. 1981). On the microscale, wood shows an orthotropic material behaviour: the stiffness and strength values are orders of magnitude higher in longitudinal direction of the fibre than in the radial or tangential direction. Since the material orientation follows the fibre course, complex stress and strain fields are observed within naturally grown wood even under simplest loading conditions (Hankinson 1921; Goodman and Bodig 1978). Thus, the morphological fluctuations have direct impact on the mechanical behaviour of wood.

The impact of knots and other growth irregularities within mechanical systems obviously depends on the dynamic boundary conditions and the length scale of the corresponding system. The concept of GLT aims at a homogenisation of the material and, in that way, a reduction of the negative aspects of the natural growing process. Nonetheless, for predicting the structural strength of wooden structures, the knot morphology, particularly in the regions exposed to tension stresses, plays an important role.

One of the first more advanced models for the strength of GLT was presented by Foschi and Barrett (1980). Therein, each lamella was divided into equally sized cells, which were assigned a random mass density and knot diameter, from which the stiffness and strength were estimated through regression.

In Ehlbeck et al. (1985), a numerical model for describing the strength of GLT beams was presented. In Colling (1990), important parameters for the GLT bending strength are discussed and a statistical model for predicting the strength of GLT beams was developed. Furthermore, Colling (1995) investigated the relation between characteristic values of GLT bending strength, board tensile strength and finger joint strength. In accordance thereto, in Falk and Colling (1995) the lamination effect – describing the homogenisation of material fluctuations from individual timber boards to GLT structures—was investigated in detail and its ranges were calculated.

Serrano and Gustafsson (2007) discussed different fracture mechanics-based approaches for assessing the strength of timber structures. The concepts are implemented into a probabilistic framework in the work of Danielsson and Gustafsson (2011), where a combination of Weibull weakest-link theory and mean stress method is proposed.

Brandner and Schickhofer (2008) investigated the relation between tensile stiffness and strength characteristics of single boards and bending strength of GLT. As result, requirements in regard to both strength and stiffness for the board material are proposed for specific GLT strength classes.

Also Frese and Blaß (2009) proposed a strength model for spruce GLT and simulated GLT beam tests to investigate the influence of single board and finger joint strength characteristics on GLT bending strength. Thereby, it was concluded that the requirements for boards and finger joints are insufficient to ensure the quality of strength classes according to EN 1194 (1999).

In Fink et al. (2013, 2014), bending tests on glued laminated timber beams with well-known material properties are reported. Subsequently, Fink and Frangi (2013) presented a model for the bending strength of GLT considering the growth characteristics observed within timber boards. Also, Fink and Kohler (2014) presented a model for predicting tensile strength and stiffness of knot clusters in structural timber.

Recently, Blank et al. (2017) discussed two GLT models: One considering perfectly brittle material behaviour and one considering quasi-brittle behaviour. Thereby, it is concluded that besides the Weibull-type size effect, which accounts for the variability of the material, a second size effect is influencing the load bearing capacity of timber. This second size effect is of deterministic nature and is characterised by the quasi-brittle failure behaviour of glued laminated timber.

The work of Kandler et al. (2015) was based on an experimental study performed at Linnæus University, Sweden. While Kandler et al. (2015) covered the linear elastic behaviour, this work focuses on structural failure mechanisms observed during the experimental study. It is the aim of this work to relate (a) the stiffness profile information presented in Kandler et al. (2015) and (b) the knot geometries reconstructed using the approach presented in Kandler et al. (2016) to the failure mechanisms observed in the experimental study.

In Sect. 2, material parameters and methodologies are presented. Section 3 discusses the results and in Sect. 4 conclusions and final remarks are given.

2 Material and methods

2.1 Material

Five different types of GLT beam setups, with ten specimens each, were produced. The five GLT beam types A, B, C, D, and E differ in number of laminations, grading class, and length. In Table 1, the parameters of the different types are summarised.

Table 1 Grading class, number of laminations, distance between the supports L and height h of the different GLT types as well as the minimum, mean and maximum values of the observed bending strength f_b and the coefficient of determination R² between maximum load bearing capacity and system stiffness

Since the study is devoted to the impact of knots on the mechanical behaviour, all beams were produced of full length lamellas. In this way, finger joints, which would have represented additional weak-sections in the GLT assembly, could be avoided.

Norway spruce boards of grades LS15 and LS22—delivered by a sawmill located in Töreboda, Sweden—were used. The original dimensions were 35 mm height, 95 mm width and 5400 mm length. Each board was marked with a unique ID, so the measured data could always be linked to the corresponding board.

For each board, the dynamic modulus of elasticity was obtained from Olsson et al. (2012)

$$\begin{aligned} E_{dyn} = 4 \rho (f_1 L)^2, \end{aligned}$$

(1)

where $\rho$ is the density, $f_1$ is the lowest axial eigenfrequency and L is the length of the board. The eigenfrequency $f_1$ was obtained using a hammer, a microphone (GRAS, type 26 CA) and a DataPhysics DP700-60 FFT-analyser.

The used boards were inspected using a tracheid effect-based laser scanning device, yielding the fibre angle on all surfaces (Nyström 2003; Simonaho et al. 2004; Petersson 2010). For a more detailed description of the procedure it is referred to Olsson et al. (2013) and Kandler et al. (2015), where the fibre angle measurements have been used to obtain a stiffness profile for each board.

After inspection, all boards were sent to a GLT manufacturer for gluing using a MUF-adhesive. Prior to gluing, all boards were planed to a height of 33 mm. Using the unique board number, the position of each board within the GLT assembly could be traced. The orientation of the boards was chosen to meet the Swedish standards for GLT production. That is, the boards were oriented such that the local pith of the top-most and bottom-most boards would be located above and below of the beam, respectively. After gluing, the beams were planed to the final width of 90 mm.

The lamella strength classes T14 (LS15) and T22 (LS22) correspond to GLT strength classes GL 24h and GL 30h, respectively (EN 14080 2013).

2.2 Reconstructing the knot geometries within wooden boards

In a first step, a geometrical description of the knot morphology is determined by employing the automated approach described in Kandler et al. (2016). Therein, fibre angle measurements were experimentally obtained from laser scans based on the so-called tracheid effect, which describes the light propagation patterns in wood. From the fibre angle measurements, the knot areas are estimated. In combination with the pith location, all knots were reconstructed using rotationally symmetric cones. It has been observed that this algorithm also works well for small knots. In Fig. 1a–c, the results for an exemplary board are shown. In Fig. 1b, all recognised knots are displayed, whereas in Fig. 1c only the significant knot groups according to a modified version of the criterion presented in Kandler et al. (2016) are shown.

The modified criterion is based on omitting knots with an opening angle $\alpha <1.5^\mathrm {o}$. Then, for all knots, the corresponding knot areas visible on the board surface are evaluated. Subsequently, all knot area values are sorted according to their value, resulting in a histogram. Then the value situated at the 70% quantile is chosen to be the critical value. All knot areas larger than this value are from now on defined as large knots, forming the 30% largest-knots quantile. These knots are assumed to be mechanically significant and are retained. In case the mutual distance in longitudinal direction x of two such knots is below 200 mm, these knots are grouped together. It is assumed that smaller knots are not important for the mechanical model and are therefore omitted. However, in case a smaller knot is in close vicinity to one of the larger knots, it is retained in the model, since it is assumed to lead to interactions in terms of fibre angle course and failure modes. As measure for the vicinity, a distance value in x-direction of 100 mm has been observed to yield good results.

The automatic nature of the reconstruction algorithm lends itself particularly useful for big sample sizes. In this study, a total number of 350 boards was investigated using this procedure.

Fig. 1

From board to stiffness profiles. a For a regular board, b the knot fitting algorithm (Kandler et al. 2016) automatically reconstructs the 3D knot geometry. c After the regrouping procedure has been applied, only the most important knot groups are retained in the model. The resulting stiffness profiles using the approach outlined in (Kandler et al. 2015) as well as the 3D FE approach (Lukacevic and Füssl 2014) for each knot group are displayed in d. In e, different strength profiles according to IPs 1 to 4 are displayed

2.3 Stiffness profiles

In this study, two different types of stiffness profiles were assessed. A stiffness profile describes the fluctuation of the local modulus of elasticity E in the longitudinal direction x, i.e. $E \equiv E(x)$. The procedure for the first type of stiffness profile is presented in Kandler et al. (2015), where fibre angle measurements obtained through laser scanning are directly used to transform the clearwood stiffness tensor. This approach is similar to the approach presented in Olsson et al. (2013), where the clearwood stiffness tensor was obtained by combining dynamic eigenfrequency measurements with a dynamic finite element analysis. Therein, the components of the stiffness tensor were calibrated so that the lowest eigenfrequency of the finite element model would match the experimentally obtained eigenfrequency. Subsequently, in each point the longitudinal stiffness component is extracted from the transformed stiffness tensor. Employing a nearest-neighbour interpolation, the stiffness values are homogenised over each cross section, resulting in a stiffness profile for each board. In Olsson and Oscarsson (2017), it has been shown that the weakest value in such a stiffness profile correlates well to the strength of a single board.

In Kandler et al. (2015), the approach is based on the same principle of transforming the clearwood stiffness tensor according to the fibre angle. Therein, in addition to the in-plane fibre angle, an empirical model for the dive-angle was employed and the clearwood stiffness tensor was obtained from a micromechanical model based on a multiscale homogenisation procedure (Hofstetter et al. 2007). Based on knowledge of the stiffness profiles in combination with the position of each board, an accurate mechanical model for assessing the stiffness of GLT beams could be developed. In Kandler et al. (2015), a linear 2D finite element code was employed to predict the effective stiffness of GLT beams. In this study, on the other hand, it is investigated how well these data can be correlated to failure mechanisms in GLT.

Complementing the approach described in the previous paragraph, the reconstructed 3D knot geometries as described in the previous section are also considered. The 3D knot geometries of each knot group were used to compute the 3D fibre angle course within the volume of the board using the procedure outlined in Hackspiel et al. (2014) and Lukacevic and Füssl (2014). This procedure follows the so-called grain-flow analogy (Foley 2003). Therein, the fibre angles in the longitudinal-tangential plane are assumed to follow the trajectories of a fluid in laminar flow where knots constitute elliptical obstacles. For the dive-angle, polynomials fitted to photographs of knot sections are used. Subsequently, the knot geometries and numerically obtained 3D fibre angles are used within a linear elastic 3D finite element analysis. Within this analysis, each knot group (in total more than 3000 groups) was loaded in tension to estimate its effective stiffness. At one end, the knot group is fixed in longitudinal direction and at the other end all nodes are exposed to a predefined displacement in longitudinal direction. The remaining boundary conditions were applied such that the specimen could freely expand in the other directions. From the resultant forces, the effective elastic modulus in longitudinal direction of each knot group could be computed. Based thereon, an alternative stiffness profile was developed, as can be seen in Fig. 1d. From Fig. 1d, it can be seen that for the knot sections a notable stiffness difference exists. The main reason is that the approach by Kandler et al. (2015) considers the local fibre deviation in each measured point but it does not consider interactions between knots. The 3D approach, on the other hand, inherently takes account of interactions between knots. Therefore, the 3D approach generally leads to lower stiffness values.

2.4 Strength profiles

For the purpose of assessing the strength properties of individual boards, a set of indicating properties is computed from the reconstructed knot geometry. Based on the study presented in Lukacevic et al. (2015), different combinations of the following parameters were considered:

• Knot-area-ratio (KAR), which is the ratio between the knot area of all knots of a knot group projected on the cross section and the cross section area.

• Knot area, which denotes the knot area visible on the surface of the board. In addition thereto, the weighted knot area distinguishes between top-bottom and left-right sides of the board,

• Knot volume.

• Interfaces between knots and the surrounding wood matrix.

• Foley-area-ratio (FAR), which—in analogy to the KAR value— describes the ratio of the area disturbed by fibre deviations to the cross section area.

In total, the four parameter sets

• 1. IP 1 $\sim$ KAR + weighted knot area + knot volume + interface + FAR + interaction terms,

• 2. IP 2 $\sim$ KAR + knot area + knot volume + interface,

• 3. IP 3 $\sim$ KAR, and

• 4. IP 4 $\sim$ KAR + knot area + knot volume.

were employed to predict the effective tensile strength $f_t$ of each knot group.

For example, in Fig. 2, the total knot volume in the maximum bending moment region (between the two points of load application) of all tension lamellas is displayed.

Fig. 2

Histogram of knot volumes in all tension lamellas in the region between the two load application points

Similar to the stiffness profiles, for each knot group the strength was estimated employing indicating properties 1–4. For the clearwood sections, the tensile strength is computed following the approach proposed in Hackspiel et al. (2014), which is based on scaling empirically based material properties according to results obtained by the micromechanical model for elastic behaviour.

2.5 Experiments

For each GLT beam, the dynamic stiffness was obtained using the same procedure as for single boards described in Sect. 2.1. Thereafter, each beam specimen was tested under quasi-static four-point-bending until failure. The four-point-bending setup was chosen to closely resemble the procedure described in EN 408 (2010). The bending span was adjusted to the length of each type, see Table 1 and Fig. 3. At load application locations, a load cell recorded the corresponding forces $F_\mathrm {left}$ and $F_\mathrm {right}$, respectively. The displacement w was recorded at the center of the beam at the compression (bottom) side. The rate of loading was approximately ${5}{\hbox {kN} \hbox {min}^{-1}}$.

Fig. 3

Four-point-bending setup and geometry. For type C, the length x was adjusted such that a knot group was located as close as possible to the centre point within the tension-side lamella. All dimensions are in (mm)

For all GLT specimens, continuous deflection measurements at the midspan point of the tension side of the beam were conducted during the loading process.

Fig. 4

Exemplary crack documentation for a GLT beam, divided into knot sections (shaded areas, indicating the sections’ knot volumes), and the manually recorded cracks (dotted and solid thick lines): a 3D view of GLT beam representation with single lamellas and their reconstructed knots, b top (=tensile loaded side), c front and d photograph of experimentally tested beam

After testing, the geometry of all cracks occurred was manually recorded and documented for all sides of the 50 GLT beams. For cases where the point of crack initiation was not straighforward to identify, all possible crack initiation locations weres documented. In Fig. 4, the documented crack pattern is displayed exemplarily for one beam. Figures of all documented results for types A to E are available as electronic supplementary material, where within each type, beams are numbered from 1 to 10.

The recorded crack measurements were used to estimate the crack area $A_\mathrm {crack}$. Assuming that all cracks range from one side of the beam to the other side of the beam, the crack area was computed according to $A_\mathrm {crack} = L_\mathrm {crack} b/2$, where $L_\mathrm {crack}$ is the accumulated crack length on all sides of the beam. For all beams, the width was b = 90 mm.

2.6 Simplified models for prediction of GLT stiffness and strength

All collected data was assembled in a database and statistically evaluated in the commercial software tool optiSLang$^\mathrm {TM}$.

In addition, an approach of directly employing the stiffness profiles E(x) was pursued. This approach is inspired by the promising procedure presented in Olsson et al. (2013), Oscarsson et al. (2014), and Olsson and Oscarsson (2017). Therein, the minimum value of the stiffness profile E(x) of a board is assumed to present the weakest defect. Consequently, this value is used as indicating property for the effective strength of the whole timber board, resulting in remarkable prediction results with values of $R^2 \approx 0.7$.

Since for most GLT beams failure is initiated in a weak section on the tension side of the beam, here the minimum value of E(x) of the outer-most lamination on the side exposed to tension is considered. It should be noted that here the stiffness profiles computed in Kandler et al. (2015) are used, which, although based on a different procedure, should yield similar results. Furthermore, the minimum stiffness value not of the whole board length but in the maximum moment range is considered. For the present load case of four-point bending this means that the minimum stiffness value between the two load application locations is considered. To account for the size of the beam, also the beam height h and the beam bending span L are considered, resulting in the linear model

$$\begin{aligned} f_b \sim 1 + \min (E(x)) + h + L. \end{aligned}$$

(2)

The model was trained using stepwise regression, where at each iteration step, predictor terms are added or removed according to a given criterion. Herein, for each term an F-test of the change in the sum of squared errors (SSE) through addition or removal of the corresponding term is used. In case the p-value of the F-statistic is smaller than 5%, the term is added to the model. Conversely, if the p-value exceeds 10%, the term is removed from the model. The bending span L and quadratic terms have been considered but did not result in significant improvements and, therefore, were omitted, resulting in the estimated bending strength

$$\begin{aligned} \hat{f_b} = c_0 + c_1 \min (E(x)) - c_2 h, \end{aligned}$$

(3)

where the regression coefficients, standard errors (SE) and p-values are

$$\begin{aligned} c_0&= {29.2}{\hbox {MPa}}\ (\text{ SE: } \text{3.6, } p\text{: } 2.0e-10) \nonumber \\ c_1&= 0.002\ (\text{ SE: } \text{0.0003, } p\text{: } 0.0004)\ \mathrm {and}\nonumber \\ c_2&= {-0.034}{\hbox {MPa}\,\,\hbox {mm}}\ (\text{ SE: } \text{-3.8, } p\text{: } 1.3e-9). \end{aligned}$$

The stepwise regression approach helps in systematically reducing the number of parameters to the statistically significant ones. In this example, the stepwise regression procedure revealed that for the investigated sample, the bending span L does not have a statistically significant impact on the resulting bending strength.

Fig. 5

Total load $F = F_\mathrm {left} + F_\mathrm {right}$ vs. displacement at the center point $w_c$ as well as maximum load levels vs. system stiffness $k = \varDelta F / \varDelta w$. a, b The results for types A and B, c, d the results for type C and e, f the results for types D and E, respectively, are shown

3 Discussion of results

3.1 System properties/effective properties of GLT specimens

From the plots shown in Fig. 5a, e, it can be seen that the two grading classes reach different maximum load peaks.

In accordance with EN 408 (2010), the system stiffness $k = \varDelta F/ \varDelta w$ is computed from a linear regression of the load displacement curve in the range of $0.1 F_{\max }$ and $0.4 F_{\max }$, where a coefficient of determination $R^2> 0.99$ is guaranteed. In Fig. 5b, d, f, the relation between system stiffness and maximum load are displayed. While a trend can be clearly observed, the values within each type are not very strongly correlated as the $R^2$ values given in Table 1 suggest.

For four-point-bending tests, the bending strength is computed according to (EN 408 2010)

$$\begin{aligned} f_b = f_m = \frac{3 a F_{\max } }{b h^2}, \end{aligned}$$

(4)

where $F_{\max }$ is the maximum total load, a is the horizontal distance between support and load, b is the width of the beam, and h is its height. In Fig. 6, it can be seen that the bending strength $f_b$ decreases with increasing lamination number. Not surprisingly, for beams of same lamination number, the bending strength of grading class LS22 is higher than the bending strength of grading class LS15. As can be seen from the confidence intervals indicated by the dotted lines, the difference between the median values of the two grading classes is significant at the 5%-level for the beams of 10 laminations. For the beams of 4 laminations, the 95% confidence intervals show a slight overlap. Therefore, strictly speaking, the null-hypothesis of different medians cannot be rejected at the 5%-level. However, looking at the individual sample values, still a notable difference is observed for the 4-lamella beams. The boxplots indicate three outliers: C4, C10 and E9.

Fig. 6

Boxplots of bending strength $f_b$. The lower edges of the boxplots indicate the 25% quantile, the upper edges of the boxplots indicate the 75% quantile and the middle line marks the median. The horizontal dotted lines indicate the 95% confidence interval for the median

The system-related stiffness k can be translated to the material-related effective modulus of elasticity (EN 408 2010)

$$\begin{aligned} E_\mathrm {GLT} = \frac{3 a L^2 - 4 a^3}{2 b h^3 \left( \frac{2}{k} - \frac{6 a}{5 G b h}\right) }, \end{aligned}$$

(5)

where L denotes the distance between the supports and G is the elastic shear modulus. In Kandler et al. (2015), the value for G was obtained from a micromechanical model. However, the study of Kandler et al. (2015) and more recently also Balduzzi et al. (2018) revealed that the shear modulus only has minor effect on the outcome of Eq. (5) for the investigated beams. For this reason and also to avoid introducing unnecessary error, here a constant value of $G={650}{\hbox {MPa}}$ in accordance with EN 408 (2010) is used.

For both stiffness and strength, the transition of the system-related quantities $F_\mathrm {max}$ and k to the material-related quantities $f_m$ and $E_\mathrm {GLT}$, respectively, ‘compresses’ the data. This ‘compression’ leads to a reduced coefficient of determination $R^2$ for the whole sample, as can be observed in Fig. 7. However, within each individual type, $R^2$ remains the same as a linear relationship is maintained.

Fig. 7

Transition of system-related quantities ($F_\mathrm {max}$, k) to material-related quantities ($f_m$, $E_\mathrm {GLT}$), leading to a ‘compression’ of the data

Fig. 8

Boxplots of the steepness measure $L_{\mathrm {crack},z}/L_{\mathrm {crack},x}$

3.2 Encountered failure mechanisms

In Fig. 5, the load deflection curves $F = F_\mathrm {left} + F_\mathrm {right}$ of all types are displayed. It can be seen that after an initially linear curve, 12 beams show nonlinear behaviour before the system load bearing capacity $F_{\max }$ is reached. These nonlinearities are on the one hand cracks on the tension side, resulting in a spike in the load-displacement curve and, on the other hand, plastifications on the compression side of the specimen, leading to a reduction of the load-displacement gradient. Thereafter, brittle system failure due to initiation of cracks is observed. The transition from a linear to a nonlinear curve can be explained by local plastification-like effects in the compression zone of the beams, as can be seen in Fig. 4a. Computing $f_b$ according to Eq. (4) does not reflect these local plastifications which lead to a non-linear normal stress distribution over the cross-section height, therefore rendering Eq. (4) inaccurate. Rather, $f_b$ has system character and represents a stress quantity corresponding to conventional brittle strength theory (Blank et al. 2017).

After a first crack has formed, some beams reach a higher load bearing capacity. This behaviour is observed for 2 beams of type A, 4 beams of type B and 3 beams of type C. For type E, one specimen showed abnormal behaviour. While the sample size may be too small to give a reliable assessment, the bigger beam types seem to show a more brittle failure behaviour than the smaller beam types, as was also found in Blank et al. (2017) (Table 2).

Table 2 Name, average mass density $\rho _\mathrm {avg}$, dynamically obtained stiffness $E_{\mathrm {GLT,exp}}^\mathrm {dyn}$, effective stiffness $E_{\mathrm {GLT,exp}}$ under quasi-static four-point-bending, ultimate bending strength $f_b$ and number of failed lamellas $n_\mathrm {lam,failed}$

To identify patterns in the crack directions, for each segment of the recorded crack geometries the height difference $\varDelta z$ between end and start point were computed. Subsequently, for each beam the sum of those values was computed to obtain a measure for the steepness of the cracks: $L_{\mathrm {crack},z} = \sum \varDelta z$. Similarly, the component related to the x-direction, $L_{\mathrm {crack},x}$, was computed from the sum of differences $\varDelta x$. In Fig. 8, the ratio of those two results $L_{\mathrm {crack},z}/L_{\mathrm {crack},x}$ is displayed for each beam. Therein, it can be seen that this ratio lies in the same range for beams of grading class LS22 and seems to be independent of the number of laminations. For example, a crack that spans 1000 mm in x-direction is, on average, accompanied by a z-increment of 40 mm. Thus, such a crack crosses approximately one lamination (recall that all laminations have a thickness of 33 mm). Conversely, for grading class LS15, the ratio $L_{\mathrm {crack},z}/L_{\mathrm {crack},x}$ is significantly larger and a crack with $\varDelta X = {1000}{\hbox {mm}}$, on average, crosses at least 2 lamellas.

This behaviour can also be observed by comparing the visualizations of the crack patterns of the two grading classes for the same numbers of laminations, i.e. Figures E.1 with E.2 and Figures E.4 with Figures E.5, respectively. For the lower grading class of LS15 crack patterns seem to remain more localized with respect to their extension in longitudinal direction, which can be explained by the higher propability of adjacent weak sections compared to the higher grading class of LS22, which is emphasized by the higher density of colored patches within the plots of the former, showing the locations of knots, and also by the higher amount of blueish/darker colors, denoting higher single knot volumes and, thus, bigger knots.

The comparison of GLT beams for the lower grading class of LS15 (cf. Figures E.1 and E.4) shows that the difference in dimensions and bending spans, about twice the length and bending span for the bigger GLT beams, leads to almost twice the steepness of the cracks. This behaviour might be explained by the fact that for the smaller dimensions, the extension of the crack in vertical direction is limited by their height, as after failure of just two laminations already half of the beam’s cross section is cracked. For the bigger GLT beams of LS15, the crack, which as explained for this grading class tends to be more localized and, thus, to propagate more likely in vertical direction, leads to comparatively more failed laminations.

Interestingly, as mentioned above, the steepness measure for the beams of LS22 (cf. Figures E.2, E.3 and E.5) seems to be the same for all dimensions and numbers of laminations. This means in the present case that as the bending span gets bigger and also the extension of the crack in longitudinal direction $L_{\mathrm {crack},x}$ is getting bigger, more laminations fail. But still even for the biggest beam type only the first few tensile-loaded laminations fail, such that the steepness measure stays the same and the size of the crack only depends on the dimensions and/or bending span.

3.3 Models for the bending strength

For the linear regression model based on minimum values of the stiffness profiles shown in Eqs. (2) and (3), an adjusted coefficient of determination of $R_{\mathrm {adj}}^2=0.6$ and root mean squared error (RMSE) of $\sqrt{\mathrm {MSE}}={5.62}{\hbox {MPa}}$ were obtained. In Fig. 10a, the values predicted from the regression model are plotted against the actual values. It can be seen that lower strength values tend to be overestimated while higher strength values tend to be underestimated by the criterion.

In addition, a ‘stiffness-profile-curvature’ was introduced to also model the spatial vicinity of neighbouring weak spots (see Fig. 9). Starting at the topmost lamella 0 (at the tension side), the lowest stiffness value in the region of the maximum bending moment is determined. For the next lamella 1, all local minima are determined, and the one located closest to the original weak spot is selected. Starting from lamella 1, the next weak spot in lamella 2 is searched and so on. Finally, a gradient is estimated by linear regression through the determined points. The idea of this approach is that the gradient represents the crack pattern responsible for failure.

Fig. 9

Examplary result of the computed stiffness profile curvature. In this case, the curvature (marked over top four layers in stiffness profile plot) is computed from the 4 topmost lamellas

Fig. 10

Predicted vs. experimentally obtained bending strength $f_b$, a using the linear regression model given in Eq. (2), b using stiffness and strength profiles (for three different IPs) in combination with a Tsai-Wu failure criterion and c using stiffness and strength profiles (for IP 3) in combination with a Tsai-Wu failure criterion and a mean stress approach

Alternatively, a more elaborate model employing a 2D finite element analysis is employed. For this, an approach similar to the mechanical model in Kandler et al. (2015) was chosen. Instead of the continuous laser-scan based stiffness profiles, the 3D FE based stiffness profiles according to Fig. 1d are used for describing the longitudinal stiffness of each lamella. In addition, the strength profiles are used to describe the tensile strength of each lamella.

The material properties are extracted from the set of stiffness and strength profiles that are provided for the FE procedure. The stiffness profiles provide the value for the longitudinal stiffness $E_L$. The remaining parameters of the transversal isotropic material are defined in accordance with the values found in EN 14081 (2009) and Neuhaus (1981):

$$\begin{aligned} G= {650}{\hbox {N}\hbox {mm}^{-2}},\nonumber \\ E_R= \frac{E_L(x)}{15},\nonumber \\ \nu _{RL}= 0.41, \nonumber \\ \nu _{LR}= 0.027. \end{aligned}$$

(6)

The values for $E_L(x)$ are obtained from the stiffness profile of the corresponding lamella. Thus, for the plane stress problem, in each integration point the elasticity matrix $\mathbb {C}$ is computed from

$$\begin{aligned} \mathbb {C} = \left[ \begin{array}{lll} 1.011 E_L(x) \quad & 0.027 E_L(x) \quad & 0 \\ 0.027 E_L(x) \quad & 0.067 E_L(x) \quad & 0 \\ 0 \quad & 0 \quad & 650.0 \end{array} \right] , \end{aligned}$$

(7)

where $E_L(x)$ is the stiffness profile value of the corresponding lamella. Similarly, each integration point is associated with strength parameters which are obtained from the corresponding strength profile. The results returned by the FE solver comprise the displacement values of all nodes as well as the stresses in all integration points. In addition, a Tsai-Wu criterion (Dorn 2012) is used to assess the utilization and load bearing capacity of the beam. Retaining the terms for the plane stress state results in the Tsai-Wu criterion as follows

$$\begin{aligned} \varPhi (\varvec{\sigma }) = a_{LL}\sigma _{LL} + a_{RR}\sigma _{RR} + b_{LLLL} \sigma _{LL}^{2}\nonumber \\ + b_{RRRR} \sigma _{RR}^2 + 2 b_{LLRR} \sigma _{LL} \sigma _{RR} \nonumber \\ +4 b_{LRLR} \tau _{LR}^2 - 1 \le 0. \end{aligned}$$

(8)

Thereby, L corresponds to x and R corresponds to the z-direction. Since the tensile strength values, represented by strength profiles, vary spatially, the corresponding parameters are dependent on the location of the integration point. The components in L-direction are computed in each integration point according to

$$\begin{aligned} a_{LL}&= \frac{1}{f_{t,L}(x)}+\frac{1}{f_{c,L}}, \end{aligned}$$

(9)

$$\begin{aligned} b_{LLLL}&= -\frac{1}{f_{t,L}(x)\ f_{c,L}}, \end{aligned}$$

(10)

where $f_{c,L}= -52.0{\hbox {N}}/{\hbox {mm}^{2}}$ and $f_{t,L}(x)$ is obtained from the strength profile of the corresponding lamella. The coefficients which are constant within the beam are (Dorn 2012)

$$\begin{aligned} a_{RR}&= 0.08564\ \frac{1}{{\hbox {N}/{\hbox {mm}^2}}}, \nonumber \\ b_{RRRR}&= 0.06588\ \frac{1}{({\hbox {N}/{\hbox {mm}^2}})^2},\nonumber \\ b_{LLRR}&= 0.0\ \frac{1}{({\hbox {N}/{\hbox {mm}^2}})^2},\nonumber \\ b_{LRLR}&= 0.01181\ \frac{1}{({\hbox {N}/{\hbox {mm}^2}})^2}. \end{aligned}$$

(11)

In accordance with the findings presented in Serrano and Gustafsson (2007), a mean stress approach is pursued. Consequently, the stress components $\sigma _{11}$, $\sigma _{22}$, and $\sigma _{12}$ are not compared directly within each integration point. Rather, the mean values of these components are computed for each cell of a rectangular grid (cell height of 43 mm and cell length of 79 mm). The cell dimensions were found by employing an optimization scheme, where the chosen dimensions were found to yield the highest coefficient of determination $R^2$. The resulting mean values are subsequently used within the Tsai-Wu failure criterion. In comparison to a strictly point-related evaluation, the mean stress approach leads to higher estimations of the total system load bearing capacity.

A comparison of the corresponding numerical and the experimental results for the bending strength $f_b$ is given in Fig. 10b. Therein, the results using the procedure with the four different IPs for the tensile strength property are shown. Results for IP 1 were omitted as the results did not show usable agreement. However, while the correct trend can be observed for IPs 3 and 4, the prediction quality is unsatisfying, the highest coefficient of determination being $R^2=0.40$ (IP 3). Using the mean stress approach (see Fig. 10c), higher numerical bending strength values are achieved, leading to a coefficient of determination of $R^2=0.54$, which is still not reliable enough. Therefore, it can be concluded that although the system failure behaviour can be interpreted as brittle failure, such brittle mechanical models do not concur well with the experimental observations. This observation also agrees with the findings presented in the work of Blank et al. (2017).

3.4 Statistical evaluation of the data

Fig. 11

Linear correlation coefficients $\delta$ and anthill plots for input and output parameters and combinations. The general input parameters are height h, bending span L, moisture content (MC), mass density $\rho$. Stiffness-based parameters are the minimum stiffness profile value following the approach by Kandler et al. (2015), the combined minimum stiffness profile and height model proposed in Eq. (3) as well as the minimum stiffness profile value based on the 3D FE calculation. $E_\mathrm {GLT,exp}^\mathrm {dyn}$ is the dynamically obtained GLT beam stiffness measure. The response parameters are the experimentally obtained effective elastic modulus $E_\mathrm {GLT,exp}$, the bending strength $f_b$ and the number of failed lamellas $n_\mathrm {lam,failed}$

Subsequently, a statistical evaluation of the data is performed using the commercial software package optiSLang$^\mathrm {TM}$. A general overview of linear correlation coefficients $\delta$ and distributions of both parameters and results $E_\mathrm {GLT,exp}$, $f_b$ and $n_\mathrm {lam,failed}$ is given in Fig. 11. The linear correlation coefficient between two quantities X and Y is computed according to

$$\begin{aligned} \delta _{XY}&= \frac{\mathrm {COV}(X,Y)}{\sigma _X \sigma _Y} \nonumber \\&\approx \frac{1}{N-1} \frac{\sum _{i=1}^{N} (x_i-\hat{\mu }_X) (y_i-\hat{\mu }_Y)}{\hat{\sigma }_X \hat{\sigma }_Y} \end{aligned}$$

(12)

where $\mathrm {COV}(X,Y)$ is the covariance between two variables, $\sigma _X$ the standard deviation, $x_i$ the i-th measurement of variable X, N is the sample size and $\hat{\mu }_X$ the estimated mean value and $\hat{\sigma }_X$ the estimated standard deviation of the corresponding variable. As for the medium-sized GLT beams (Type C) no experiments for the lower grading class were conducted and, in addition, not all parameters for the higher grading class were available, in Fig. 11, which correspond to the 3D FE and the strength profile parameters, the results only for Types A, B, D and E are shown. The data is clustered into general parameters and specific parameter groups as follows:

• General parameters General parameters cover bending span L and height h of the beam and also the average moisture content (MC). Also, the average mass density $\rho$ of the topmost (tension) lamella is included. Regarding the correlation within this group of parameters, mass density $\rho$ and moisture content show a linear correlation coefficient of 0.78. This can be explained by increasing weight (and, therefore, increasing values for mass density measurements) of wood with increasing MC. The relationship between these parameters is visualised in Fig. 11b.

• Knot morphology parameters The investigated knot morphology parameters comprise the knot volume, the knot area visible on the board surface and the knot interface area to the surrounding wood matrix. Herein, for each parameter the total sum of all knots of the topmost (tension) lamella occurring between the two load application points is used. The linear correlations between knot volume, visible knot area and interface area are between 0.87 and 0.99. Consequently, the correlation with the quantities of interest $E_\mathrm {GLT,exp}$ and $f_b$ is approximately the same for these parameters, as can be seen in the three rightmost columns in Fig. 11c. It can be noticed that all three parameters show correlation to the beam length L and beam height h. The reason for this behaviour is the increasing distance between load application points with bigger beam dimensions, see also Fig. 3, which, in turn, leads to an increasing total sum of knot morphology parameters. The distance to the pith did not yield any notable linear correlation to the remaining parameters.

• Stiffness related parameters The stiffness related parameters represent the stiffness profiles computed according to the model presented in Kandler et al. (2015) as well as the 3D FE approach. For both stiffness profile types, the minimum value occurring in the tension lamella between the load application points is used as parameter. Also the stiffness profile curvature corresponding to Sect. 3.3 as well as the regression model in Eqs. (2) and (3) belong to this parameter group. Unsurprisingly, the regression model parameter is strongly correlated to the stiffness profile parameter. A notable correlation can be observed between the two stiffness profile parameters and mass density and moisture content measurements. The reason for this observation lies in the micromechanical model (Hofstetter et al. 2005) which was employed for computing the clearwood stiffness tensor within Kandler et al. (2015). For the micromechanical model, mass density and moisture content are the two main input parameters. Also, the two stiffness profile parameters show notable correlation to the knot morphology parameters. The knot morphology can be interpreted as a latent factor influencing both the knot morphology parameters as well as the stiffness profile computation. While the knot morphology is not directly used within the computation of the stiffness profiles, it influences the fibre deviations (Foley 2003) and thus an important aspect of the stiffness profile computation presented in Kandler et al. (2015).

• Strength related parameters Strength related parameters represent the strength profiles computed in accordance with Sect. 2.4. In analogy to the stiffness related parameters, for each beam the corresponding strength parameter was defined as the minimum value occurring in the tension lamella between both load application locations. Since all four indicating properties are based on the same knot morphology parameters, they are—with exception of IP 1—very closely correlated. For the same reason they show notable correlation to the 3D FE stiffness parameter. However, the linear correlation values $\delta$ are rather low, ranging from 0.36 (IP 1) to 0.51 (IP 4). The correlation between the 3D FE stiffness profile-based criterion and the strength $f_b$ is 0.47.

• Dynamic stiffness measurement Between the dynamically obtained stiffness value $E_\mathrm {GLT,exp}^\mathrm {dyn}$ and the remaining input parameters, the highest linear correlation value is observed to the stiffness profile parameter.

• Output/result parameters For the results, the input parameters are sorted according to their linear correlation coefficients with the quantities of interest in Fig. 12. For the effective bending stiffness $E_\mathrm {GLT,exp}$, the parameter with the highest correlation of 0.95 is the dynamically obtained stiffness $E_\mathrm {GLT,exp}^\mathrm {dyn}$, followed by the stiffness profile parameter according to Kandler et al. (2015). For the remaining parameters, e.g. moisture content, mass density etc., the correlation is relatively low and it is not expected to gain reliable predictions from those parameters alone. Regarding the bending strength $f_b$, the regression model (2) achieves a correlation coefficient of 0.79, corresponding to a coefficient of determination (CoD) of $R^2 = 0.62$. The stiffness profile criterion yields a correlation coefficient close to 0.71, corresponding to $R^2=0.50$. While exposing a clear trend, these results indicate that for the reliable prediction of bending strength, more sophisticated models have to be employed. Interestingly, the knot morphology apparently shows better correlation to the bending strength than the applied indicating properties. For the number of failed lamellas, $n_\mathrm {lam,failed}$, no meaningful correlation could be identified.

Fig. 12

10 highest linear correlation coefficients between parameters and results for a effective bending stiffness $E_\mathrm {GLT,exp}$, b bending strength $f_b$ and c number of failed lamellas $n_\mathrm {lam,failed}$. h_Kandler2015 denotes the approach given in (3)

4 Conclusion

To study the impact of knots on the mechanical behaviour of GLT beams, several parameters were investigated with regard to their ability to predict effective GLT beam properties.

The use of the laser scanning information of all boards, obtained before their assembly to GLT beams, allowed a virtual reconstruction of the knot morphology and led to knowledge of the knots’ locations, orientations and sizes within such beams. Furthermore, real failure mechanisms of experimentally tested GLT beams were investigated, by recording crack patterns, which were used to obtain measures of cracks.

The main findings can be summarized as follows:

• 1. With an increasing number of laminations, the bending strength $f_b$ decreases. In addition, the grading class also had the expected influence of lowering the bending strength. It was also observed that more laminations lead to more brittle failure modes. With respect to this effect it would be interesting to investigate beams with constant height and varying number of laminations in the future.

• 2. The load-displacement curves confirmed that bigger beam types show a more brittle failure behaviour than smaller ones. In addition, the steepness of the observed cracks ($L_{\mathrm {crack},z}/L_{\mathrm {crack},x}$) seemed to be independent of the number of laminations for the higher grading class of LS22. For LS15, the crack patterns were more localized, with smaller extensions in longitudinal direction, such that the steepness of the crack was controlled by the number of laminations. Thus, a higher number of such laminations led to steeper cracks. This behaviour can be explained by the higher probability of adjacent weak zones in adjacent laminations for the lower grading class.

• 3. The best estimate for the effective GLT beam bending stiffness, besides a correlation with experimentally obtained dynamic stiffness values for assembled beams, was achieved by correlating them to minimum values of stiffness profiles of the outer-most lamination on the side exposed to tension.

• 4. This stiffness profile parameter also led to the best correlations to the effective beam bending strength ($R^2=0.62$). Here, an overestimation of lower strength values and an underestimation of higher strength values was observed. In the application of strength profiles to linear elastic FE models the correct trend for the estimation of strength values could be found, but the resulting coefficient of determination was rather low.

It can be concluded that any simple correlation model is unable to reliably predict the failure behaviour of GLT beams and that more sophisticated models are needed.

Especially the last-mentioned approach of using stiffness and strength profiles, i.e. of fluctuating section-wise effective material properties for all boards, based on virtual reconstructions of boards and the evaluation of indicating properties, within simple linear-elastic FE simulations, showed the problems of such a method: complex morphologies and failure mechanisms of knot sections cannot be easily replaced by sections with effective strength values if failure mechanisms are not taken into account. Therefore, in the future, distinct mechanical models for timber elements with section-wise effective properties must be developed. This could for example be achieved by additionally assigning fracture energies to the individual sections, as also suggested by Blank et al. (2017). They already showed, by using linear traction-separation laws in combination with constant fracture energies, good fits to experiments could be obtained. In future, more realistic, variable fracture energies for knot sections might be obtained by simulating those sections with a multisurface failure criterion (Lukacevic and Füssl 2016; Lukacevic et al. 2017), whose combination with the mentioned fiber deviation model shows promising results. Furthermore, in such an approach, additional strength reducing or failure inducing factors, like finger joints, could be easily implemented by introducing additional weak sections.

Particularly, such an approach is needed if the structure’s failure mechanism is rather brittle. In contrast, for ductile failure mechanisms, like observed within cross-laminated timber plates under point loads, it could be shown that the application of strength profiles in the framework of limit analysis (Füssl et al. 2017; Li 2018), i.e. in combination with perfect plasticity models, leads to reasonably effective strength estimates.