VIBRATIONAL PROBLEMS OF TIMBER BEAMS WITH KNOTS CONSIDERING UNCERTAINTIES

Journal of the Brazilian Society of Mechanical Sciences and Engineering manuscript No. (will be inserted by the editor) VIBRATIONAL PROBLEMS OF TIMBER BEAMS WITH KNOTS CONSIDERING UNCERTAINTIES Diego A. Garcı́a · Rubens Sampaio · Marta B. Rosales Received: date / Accepted: date Abstract A stochastic model of the dynamic behavior of sawn timber beams of Argentinean Eucalyptus grandis is herein presented. The aim of this work is to study the influence of the timber knots in the dynamical response of timber beams. The presence of knots is known to be the main source of the lengthwise variability and reduction of bending strength and stiffness in timber beams. The following parameters of the timber knots are considered stochastic: position along the beam span and within the beam cross section, shape and dimensions. Experimental data, obtained from bending and density tests, are employed to find the timber modulus of elasticity (MOE) and density. On the other hand, the characteristic of the timber knots used in the stochastic model was obtained from a visual survey performed with 25 beams of the same species . The problem of the natural vibration frequencies of the timber beam is approximated with the finite element method. Numerical results are obtained using Monte Carlo Simulations (MCS). The uncertainties of the timber knots parameters and their influence are quantified with the probability density function (PDF) of the frequencies. Statistical results obtained by means of MCS are compared with experimental measurements in order to assess the accuracy of the stochastic model. The present approach that gives a more realistic description of timber structures, correlates better with experiments. Keywords Timber beams · Knots · Uncertainties · Natural vibration frequencies · Natural vibration modes · Composite materials D. A. Garcı́a · M. B. Rosales Department of Engineering, Universidad Nacional del Sur, Av. Alem 1253, 8000 Bahı́a Blanca, Argentina. CONICET, Argentina. Tel.: +54-0291-4595156 E-mail: garciadiego@fio.unam.edu.ar R. Sampaio Department of Mechanical Engineering, PUC-Rio, Rua Marquês de São Vicente, 225 22453-900 Rio de Janeiro RJ, Brazil. 2 Diego A. Garcı́a et al. 1 INTRODUCTION All trees have branches which start from the pith at the heart center of the trunk. Since the fibers on the upper side of branches are not interlocked with fibers of the trunk, the branches can break off, due to heavy fall or strong wind. Later, the wound caused by a broken off branch is, in some way overgrown by wood tissue and the fibers of the trunk become continuous again. Thus, the knots are unavoidable. When the tree trunk is converted into structural timber, cuts destroy the interlocking of the fibres and knots are created, leaving uncoupled fibres. Knots are related to grain distortion and disruption in the fiber continuity. Cracks are common in the surrounding material due to the stress concentrations caused by the difference in the physical properties of the knot and the normal wood tissue. Due to the specific structure of trunk and branches, structural timber can be characterized as a composition of clear wood and growth defects, a natural composite material. Clear wood is an anisotropic material, but its properties do not change considerably along grains. On the other hand, growth defects, such as knots, often related to localized grain deviations, are the main source of the lengthwise variability in the bending strength and stiffness in timber beams. The presence of grain deviations can decrease the MOE in the longitudinal direction. Since knots are unavoidable in structural timber, the effective MOE in the longitudinal direction varies along the main axis of a beam. The reduction of the timber strength and stiffness due to the presence of knots depends on the size of the knots, their type and their location. Eucalyptus grandis, which is mainly cultivated in the Mesopotamian provinces of Entre Rios and Corrientes, is one of the most important renewable species cultivated in Argentina. A simple method for visually strength grading sawn timber of these species has been developed by Piter [15]. According to Piter, the presence of pith, or medulla, often associated with other defects as fissures, significantly reduces the strength and the stiffness of this sawn timber. This feature is also considered the most important visual characteristic for strength grading this material by the Argentinean standard IRAM:9662-2 [11]. Other important features taken into account in grading are the knot ratio and the grain deviation [11, 15]. Due to the variability in the mechanical properties, a stochastic approach appears desirable to attain a more realistic structural model. The influence of the knots in the structural behavior of timber beams was first considered by Czmoch [6]. He studied the bending strength in sections with knots and determined the load carrying capacity of timber beams. The presence of knots was modeled through a Poisson process. Escalante et al. [7] studied the buckling of Eucalyptus grandis wood columns with the finite element methodology and the lengthwise variation of MOE was modeled as a Gaussian random process. They also applied a Karhunen-Loève (KL) expansion in order to discretize the random field. Köhler et al. [14] reported a probabilistic model of timber structures where the MOE and the mass density were represented by random variables with a lognormal and a normal PDF, respectively, assuming a homogeneous value within a structural element. Köhler [13] presented a discrete model of the lengthwise variability of bending strength taking into account the presence of timber knots following the model for bending moment capacity proposed by Isaksson [12]. In this work, the discrete sec- VIBRATIONAL PROBLEMS OF TIMBER BEAMS 3 tion transition was assumed to be Poisson distributed. Thus, the length between knots follows an exponential distribution. Probabilistic models of timber materials properties introduced by Köhler et al. [14] and Faber et al. [8] have been introduced in the context of the Probabilistic Model Code (PMC) of the Joint Committee on Structural Safety (JCSS) for the probabilistic design of timber structures. Baño [2] presented a study in which timber beams with defects are simulated and their maximum load in bending is predicted . The development of a bi-dimensional model of timber pieces free of defects in order to predict the performance of timber structural elements was reported by Baño [3]. Baño [4] analyzed the influence of the size and position of cylindrical knots on the load capacity of timber elements using a FE program. Guindos [9] studied a three-dimensional wood material model implemented in a finite element (FE) software which is capable of predicting the behavior of timber at the macro-scale taking into account the effect of any type of knot. They are modeled as oblique cones. Then, Guindos [10] analyzed the influence of different types of knots and fiber deviations on the bending of wood, using visual grading standards, by means of the FEM. The mechanical properties of the material in these works [2–4, 9, 10] were found in accordance to the guideline established by the standard UNE-EN 408 [18]. The aim of the present work is to quantify the influence of the timber knots and the variation that they produce in the mechanical properties and the effect in the dynamic behavior of sawn timber beams of Argentinean Eucalyptus grandis. This influence is quantified in the first three natural frequencies and in their associated modes of vibration. To accomplish this, the PDFs of the first three natural frequencies were found via Monte Carlo simulations (Rubinstein [16]). The lengthwise variability of the MOE and of the second moment of area of the beam cross section is introduced to account for the knots. Frequently, the presence of knots in structures made out of sawn timber is disregarded, maybe due to the lack of data or of an appropriate model. Timber knot parameters are modeled via the joint probability mass function (joint PMF) obtained with experimental data from visual survey of beams of Eucalyptus grandis with structural dimensions. In this study, timber knots are modeled as holes in the beam cross section, hence considered in the second moment of area of the beam cross section. This is due to the grain deviation produced by the knot presence and also because the timber fibers of the trunk have a different orientation than the fibers of the branches. The local reduction of the MOE due to the grain deviation is also considered. The lengthwise variability of the MOE, presented in this work, was developed starting from the weak-zone model [6,12] with modifications in the length of the weak-zone, which in this work, is considered proportional to the greater dimension of the knot. The model herein proposed also takes into account the modification in the beam cross section due to the knot presence. It should be noted that despite the knot is modeled as a hole when the bending stiffness is calculated, no reduction of mass is made in the inertial terms. The PDFs of the MOE and the mass density are obtained by means of the Principle of Maximum Entropy (Shannon [17]), and its parameters by means of the maximum likelihood method (MLM) applied to MOE values that were obtained experimentally. Additionally and in order to measure the fit between the experimental and 4 Diego A. Garcı́a et al. Fig. 1 A simply supported sawn beam with knots. theoretical PDFs of the MOE and the mass density, the Kolmogorov-Smirnov (K-S) and Anderson-Darling (A-D) tests of fit are used. Numerical results of simply supported sawn beams of Argentinean Eucalyptus grandis are presented and discussed. The PDFs of the natural frequencies and the modes of vibrations are reported. The influence of the knot modeling is evaluated. Finally, a comparison between results of numerical simulations and experiments is also reported. Thus, the accuracy of the stochastic model herein presented was assessed. 2 PROBLEM STATEMENT The natural vibration problem of a simply supported sawn beam of Argentinean Eucalyptus grandis with knots is herein described (Fig. 1). For the free vibration problem, the well-known differential equation for a EulerBernoulli beam is:   ∂ 2 v(x,t) ∂2 ∂ 2 v(x,t) =0 (1) + 2 e(x)i(x) ρ(x)a(x) ∂t 2 ∂x ∂ x2 where ρ(x) is the material density per unit of length, a(x) is the beam cross section, e(x) is the Modulus of Elasticity (MOE), i(x) is the second moment of area of the beam cross section, v(x,t) is the transverse displacement, x is the position within the beam span and t is the time variable. In the present work, the lengthwise variabilities of the MOE and of the second moment of area of the beam cross section are introduced to account for the presence of timber knots that produce a local reduction of both. Random variables are used and, in what follows, these stochastic quantities will be denoted with capital letters. The differential equation Eq. (1) becomes:   ∂ 2V (x,t) ∂2 ∂ 2V (x,t) ρa + 2 E(x)I(x) =0 (2) ∂t 2 ∂x ∂ x2 The results will be reported for pinned-pinned boundary conditions. At x=0 and x=L, the deflections and bending moment are zero. In this study, timber knots are modeled as holes in the beam that modify the second moment of area of the beam cross section and consequently the bending stiffness. However, no holes are considered in the mass since the knot mass participates in the inertial terms, assuming that the material inside the knots have similar density to the rest of the beam. VIBRATIONAL PROBLEMS OF TIMBER BEAMS 5 3 STOCHASTIC FINITE ELEMENT APPROACH If a set of admissible functions ψ is prescribed Eq. (1) can be written in a variational formulation context as:   Z L ∂ 2 v(x,t) ∂ 2 ∂ 2 v(x,t) ρ(x)a(x) φ (x)dx = 0 ∀φ (x) ∈ ψ (3) + 2 e(x)i(x) ∂t 2 ∂x ∂ x2 0 In particular, for the pinned-pinned beam, ψ = {φ : [0, L] → R, φ is piecewise continuous and bounded, φ (0) = 0, φ (L) = 0} (4) this formulation together with the boundary conditions lead to the following form of the variational problem: M(v, φ ) + K(v, φ ) = 0 ∀φ ∈ ψ (5) where M(v,φ ) and K(v,φ ) are the mass and stiffness operators respectively, defined as follows: M(v, φ ) = ∂ 2 v(x,t) φ (x)dx ∂t 2 (6) ∂ 2 v(x,t) ∂ 2 φ (x) dx ∂ x2 ∂ x2 (7) Z L ρ(x)a(x) Z L e(x)i(x) 0 and K(v, φ ) = 0 Now, Eq. (5) is discretized using the Galerkin Method. We define a N-dimensional subspace ψ N ⊂ ψ, where a function vN ∈ ψ N . The problem can be formulated as follows: Find vN ∈ ψ N such that: M(vN , φ ) + K(vN , φ ) = 0 ∀φ ∈ ψ N (8) Applying the standard finite element methodology (see for example Bathe [5]), the variational form Eq. (8) is discretized. Euler-Bernoulli beam elements with two nodes and two degrees of freedom per node (transverse displacement and rotation, respectively) are employed. These elements are based in the following shape functions:   3 2 1 0 L 2 L 3  e e      0 1 2 1  1    Le Le    x2  n(x) =  (9)  x   2 3  3 0 0 −  Le2 Le3   x   1 1  00 Le Le2 6 Diego A. Garcı́a et al. Fig. 2 Euler-Bernoulli beam finite element with four degrees of freedom. where Le is the element length. The spatial interpolation of the transverse deflection v(x) can be written in terms of the nodal variables as v(x) = nT (x)v. (10) where   vT = v1 θ1 v2 θ2 (11) is the nodal displacement vector of the beam element (Fig. 2). Through the application of the finite element method, the components of the beam element stiffness and mass matrices are obtained: d 2 ni (x) d 2 n j (x) dx dx2 dx2 Ke,i j = Z Le E(x)I(x) Me,i j = Z Le ρani (x)n j (x)dx 0 0 (12) (13) where the random (stochastic) quantities E(x) and I(x) represent the lengthwise variability within the beam. ρ presents variability among beams though its lengthwise variability within each beam is not taken into account. Next, the global stiffness and mass matrices can be obtained with the usual finite element assembling. The natural frequencies and modes are obtained solving the Eq. (14) below: [K −Vn2 M]Φn = 0 (14) where K and M are the n x n positive-definite global stiffness and mass matrices, respectively. 4 MECHANICAL PROPERTIES OF THE MODELS In this section, we present the assumptions and the way in which the mechanical properties that appear in Eq. (2) are represented. VIBRATIONAL PROBLEMS OF TIMBER BEAMS 7 Fig. 3 Geometry of the knots types considered in this work. 4.1 Timber-knot dimensional parameters In order to simulate the timber knots, we define the joint PMF of the timber knot shape parameters within the timber beam cross section and the probability mass functions of the distance between timber knots and of its length in the direction parallel to the longitudinal axis of the timber beam. To find the joint PMF of the knots parameters, experimental data obtained from visual survey of 25 sawn beams of Eucalyptus grandis of structural size with 180 timber knots were employed. The distance between knots, their dimensions perpendicular and parallel to the longitudinal beam axis, their depth and position within the beam cross section, are the knot features reported in the visual survey. Considering these visual parameters, the timber knots are classified in four types (Fig. 3): a. Timber knots with the depth equal to the beam width and with vertical position within the beam height. b. Timber knots with the depth less than the beam width and with vertical position within the beam height. c. Timber knots with the depth less than the beam width and with vertical position near to the edge of the beam cross section. d. Timber knots with the depth equal to the beam width and with vertical position near to the edge of the beam cross section. In order to simulate the dimensions of the timber knots and their position within the cross section, a joint PMF of the three random variables is defined taking into account the parameters which define the position and dimensional characteristic of the timber knots within the beam cross section: pX,Y,Z (x, y, z) = P[Z = z | X = x,Y = y]P[Y = y | X = x]P[X = x] (15) where the random variables are, (see Fig. 4): – x is the position of the knot centroid along the height of the beam cross section. – y is the knot size along the height of the beam cross section. – z is the knot depth along the width of the beam cross section. In the type b of timber knots, the timber beams are not cut through the cross section with the knot to determine the knot depth because they become useless to perform other experimental tests. Due to the lack of information about its depth in the beam cross section, the Principle of Maximum Entropy (PME) is employed to obtain the probability distribution. The PME states that, subjected to known constraints, the joint PMF which best represents the current state of knowledge is the one with largest 8 Diego A. Garcı́a et al. Fig. 4 Beam cross section with a knot. Random variables of the joint probability mass function. Table 1 Mean values and standard deviations of the timber knots parameters. Knot dimension µ (mm) σ (mm) δ= σ /µ X Y Z U R 61.62 23.95 20.22 288.62 40.19 38.77 11.46 11.49 175.52 21.36 0.63 0.48 0.57 0.61 0.53 entropy. The measure of uncertainties of a discrete random variable Z is defined by the following expression: n S (p) = − ∑ pi ln(pi ) (16) i=1 in which pi is the probability of the discrete random variable Z which assumes n different values. It is possible to demonstrate that the application of the PME when the random variable assumes a finite number of values within the interval [a, b] without further knowledge about the random variable, leads to a uniform PMF. The random variables that define the distance between timber knots and their dimension in the direction parallel to the beam axis are, respectively: – u is the distance between timber knots. – r is the length of the knot (dimension along to the longitudinal beam axis). and they are defined respectively by the following joint PMF: pU (u) = P[U = u] (17) pY,R (y, r) = P[R = r | Y = y]P[Y = y] (18) The last expression is assumed to simulate the knot shape resulting of the visual survey. This assumption implies that the dimension of the timber knots parallel to the longitudinal beam axis r are related with the dimension of the timber knots perpendicular to the longitudinal beam axis y. In Tab. 1 the mean values, standard deviations and coefficients of variation of the random variables that represent the characteristics of the timber knots are depicted. VIBRATIONAL PROBLEMS OF TIMBER BEAMS 9 Fig. 5 Two point load bending test according to UNE-EN 408 [18]. Table 2 Eucalyptus grandis strength classes, according to IRAM 9662-2 [11]. Strength class Presence of pith Knot ratio Grain deviation C1 C2 C3 No No Yes K ≤ 1/3 1/3 < K ≤ 2/3 2/3 < K gd < 1/12 gd < 1/9 1/9 < gd 4.2 Modulus of Elasticity (MOE) To find the PDF parameters of the MOE, experimental data presented by Piter [15], obtained by means of two point load bending tests performed with 349 sawn beams of Argentinean Eucalyptus grandis with structural dimensions, are employed. Bending tests were carried out according to UNE-EN 408 [18] (Fig. 5), and the worst defects were placed in the constant bending zone, between two concentrated loads and within the tensile region of the cross section. These values of the MOE were calculated taking into account the shear deformation (global MOE Eg ). Then, the MOE values obtained experimentally are classified according to the strength classes established for the visual grading of the Eucalyptus grandis cultivated in the Mesopotamian provinces of Argentina by the standard IRAM:9662-2 [11], see Tab. 2. MOE values obtained experimentally were corrected to a uniform moisture content of 12%, in order to make the 349 values comparable. The humidity content previously established, corresponds to measurements made at a temperature of 20◦ C and a relative humidity of 65%. The values of Eg obtained with experimental data have been corrected, increasing in 2% for each 1% in excess to the standardized condition of 12% of humidity content, and vice versa, in each timber beam. On the other hand, using the expression (19), 349 values of Eg have been calculated and reported by Piter [15]. Eg = L3 (F2 − F1 ) 4.7bh3 (w2 − w1 ) (19) where (F2 − F1 ) is the load increment and (w2 − w1 ) is the midspan deflection increment corresponding to the load increment. This load increment is within the linear elastic range of the material. To determine the PDF of the MOE, the Principle of Maximum Entropy was applied. The measure of uncertainties of a random variable X that represent the MOE is defined by the following expression S( fX ) = − Z D fX (X)log( fX (X))dX (20) 10 Diego A. Garcı́a et al. Table 3 Results of the K-S and the A-D tests of fit of the MOE PDF. PDF K-S √ D n Statistics Significance test A-D A∗ Statistics Significance test Gamma Lognormal Normal 0,64 0,74 1,00 0,64 < 1,36 0,74 < 1,36 1,00 < 1,36 0,23 0,81 1,32 0,23 < 0,75 0,81 > 0,75 1,32 > 0,75 Table 4 Parameters of the gamma PDF of the MOE. Parameters MOE sections with knots MOE sections free of knots a b µ σ δ = σ /µ 42,730 0,315 13,4879 GPa 2,063 GPa 0,153 31,301 0,507 15,889 GPa 2,84 GPa 0,178 in which fX stands for the PDF of X and D is its domain. It is possible to demonstrate that the application of the PME under the constraints of positiveness and bounded second moment, lead to a gamma PDF. This is due to the fact that the domain of MOE is the positive real numbers, Ei ∈]0, ∞[, and the interval is open, i.e. the boundaries do not belong to the interval. The parameters of the PDF of the MOE are estimated by using the maximum likelihood method (MLM). Finally, the Kolmogorov-Smirnov (K-S) and the AndersonDarling (A-D) tests of fit are used (e.g. Ang and Tang [1]). These tests have been employed due to the fact that the first is more sensitive to the values closer to the median of the distribution whereas the second gives more weight to the values in the tail of the distribution. The level of significance α of the parametric hypothesis is assumed to be 0,05. For α = 0, 05, the critical values for the K-S and the A-D test of fit have been obtained from [1]. Test statistics, critical values and the results of the test of fit are presented in Tab. 3. As can be observed, the best fit is attained with the gamma PDF in agreement with the Principle of Maximum Entropy. The gamma PDF of the MOE is (Fig. 6): f (x | a, b) = x 1 xa−1 e− b baΓ (a) (21) where a and b (shape and scale parameters, respectively) are depicted in Tab. 4. 4.3 Lengthwise Variability of the MOE The structural timber is composed of clear wood and wood with defects. The knots affect the mechanical properties considerably. However, the quantitative knowledge about this relation is very scarce. In the present study, the lengthwise variability of the MOE is represented following the model of the bending strength presented by Isaksson [12] and Czmoch [6]. In the weak-zone model, the timber beam is modeled VIBRATIONAL PROBLEMS OF TIMBER BEAMS 0.2 11 MOE sections with knots MOE sections free of knots Probability 0.15 0.1 0.05 0 7 9 11 13 15 17 MOE GPa. 19 21 23 25 27 Fig. 6 Gamma PDFs of the MOE in the free of knots sections (red) and in the sections with knots (blue). as a composite of short weak zones connected by longer sections of clear wood. Weak zones correspond to knots or group of knots and are randomly distributed. In the stochastic model herein presented, the length of the weak zones are proportional to the greater dimension of the knot. This is based on the fact that the local deviation of the fibres is greater than the knot dimension and with a length proportional to its size. This feature was observed in the visual survey of the timber knots and its surrounding fibres. The MOE in each of these zones is constant and is randomly assigned. The MOE of the clear wood is assumed constant along the beam span, analogously to the model of bending strength presented by Isaksson [12]. The dimensional parameters of the beam cross section affected by the knots presence, are modified only within the length of the knot. 4.4 Mass Density In this work, the mass density was considered constant along the beam span and the lengthwise variability due to the knot presence is not taken into account. In the Eucalyptus grandis beams, the knots are frequently composed of material with similar density than the clear wood. This does not mean that a mass lengthwise variability might be present but it is not considered in the present investigation. The variability in the density among beams is considered through the PDF obtained from experimental values. To find the parameters of the PDF of the mass density (ρ), experimental data presented by Piter [15] and obtained by means of density measurement performed with 50 sawn beams of Argentinean Eucalyptus grandis was employed. The density measurements were carried out according to UNE-EN 408 [18]. Values of mass density experimentally obtained were corrected to a uniform moisture content of 12 %, in order to make comparable the 50 values that have been calculated of beams with different humidity content. The values of ρ have been corrected, increasing in 0,5 % for each 1 % in excess to the standardized condition of 12 % of humidity content, and vice verse, in each timber beam. Köhler [14] represented the variability of the mass density between structural timber beams with a Normal PDF. In the present work, the parameters of the PDF of the mass density are estimated using the MLM. Finally, the K-S and the A-D tests of fit are used to choose the most adequate PDF. As before, a level of significance 12 Diego A. Garcı́a et al. Table 5 Results of the K-S and the A-D tests of fit of the mass density PDF. PDF K-S √ D n Statistics Significance test A-D A∗ Statistics Significance test Lognormal Gamma Normal 0,65 0,69 0,75 0,65 < 1,36 0,69 < 1,36 0,75 < 1,36 0,39 0,39 0,41 0,39 < 0,74 0,39 < 0,75 0,41 < 0,74 Table 6 Statistical parameters of the lognormal PDF of the mass density ρ . Parameters Values µ σ Mean Standard deviation Coefficient of variation 6,221 0,069 504,4 Kg/m3 34,84 Kg/m3 0,069 α=0.05 is assumed. Test statistics, critical values and the results of the tests of fit are presented in Tab. 5. As can be seen in Tab. 3, the PDF with the best fit is the Lognormal: f (x | µ, σ ) = 1 √ xσ 2π e −(ln(x)−µ )2 2σ 2 . (22) The parameters µ and σ are depicted in Tab. 6. 5 NUMERICAL RESULTS In what follows, some numerical results are presented. In all the simulations, 100 finite beam elements are used. The integrals of the components of the element stiffness matrix, Eq. (12), are computed by means of the Gauss quadrature using five points. The dimensional parameters of the timber knots are simulated with the inverse transform method (Rubinstein [16]). 5.1 First Model M1 In the first model, timber knot dimensional parameters are modeled as random variables through the joint PMF previously presented. The MOE in the free knots zone of the beam and in the zone with knots are considered deterministic. The length of the weak zone is assumed equal to seven times the greater dimension of the knot. This feature was observed in the visual survey of the timber knots and its surrounding fibres. The data of the simulation is depicted in Tab. 7. These dimensional parameters correspond to timber beams of structural size, that are often used in design practice and are within the dimensions of the timber beams used in the visual survey of timber knots parameters. MOEFKS stands for the modulus of elasticity of a free knots section and MOESW K for a section with knots. VIBRATIONAL PROBLEMS OF TIMBER BEAMS 13 Table 7 M1. Parameters used in numerical simulation. Parameters Values Length Nominal section MOEFKS MOESW K Mass Density 3m 50 x 150 mm 15,889 GPa. 13,489 GPa. 505 Kg/m3 0.1 F 1 F2 0.08 e% F 3 0.06 0.04 0.02 0 400 600 800 1000 1200 ns 1400 1600 1800 2000 Fig. 7 Convergence of the mean values E[Fn ] for the first three natural frequencies (Model M1). Error e% from Eq. (23) . A convergence study is shown in Fig. 7 where N is the number of independent Monte Carlo simulations that gives a prescribed accuracy. The adopted convergence criterion is the following: e% = E[Fnns ] − E[Fnns−200 ] E[Fnns−200 ] % < 0, 2% (23) where E[Fnns ] is the mean value for the number of simulations ns and E[Fnns−200 ] is the mean value for the number of simulations ns − 200. This criterion is adopted due to the simple shape of the natural frequencies PDFs. An acceptable convergence is achieved rapidly, when N=600, as can be observed in Fig. 7. Next, for this number of independent Monte Carlo simulations, the PDF of the first three natural frequencies (f[Fn ]) are obtained (Fig. 8). It can be seen that the shape of the three PDFs are approximately equal through the range of variation of the frequency which, in turn, increases from the first to the third natural frequency. In Tab. 8, the mean values, standard deviations and coefficients of variation of f[Fn ] are presented. As can be observed, the mean values and standard deviations increase from the first to the third natural frequency, while the coefficients of variations remain constant for the three PDF presented in Fig. 8. The aim of this first model M1 is to quantify the influence of the timber knots without taking into account the variability of the MOE between the beams, being the stochasticity of the model only originated by the knots parameters. Unlike, in the design practice, the influence of the knots is only taken into account in the value of the selected MOE (Tab. 2) through the strength class and not in the others features as the decrease of the second moment of area or the local reduction of the MOE derived 14 Diego A. Garcı́a et al. 1.8 0.45 0.2 1.6 0.4 0.18 1.4 0.35 1.2 0.3 1 0.25 0.16 0.14 0.8 f[F3] f[F2] f[F1] 0.12 0.2 0.1 0.08 0.6 0.15 0.4 0.1 0.2 0.05 0 35 36 37 F1 Hz 38 39 0 140 0.06 0.04 0.02 145 150 155 0 320 F2 Hz 330 340 350 F3 Hz Fig. 8 PDF of the first three natural frequencies found with M1. Table 8 Natural frequencies (M1). Mean values, standard deviations and coefficients of variations of f[Fn ]. Natural Frequency E[Fn ] σ [Fn ] δ [Fn ] = F1 F2 F3 37,435 Hz 149,677 Hz 336,707 Hz 0,170 Hz 0,716 Hz 1,576 Hz 0,004 0,004 0,004 σ [Fn ] E[Fn ] Table 9 Natural frequencies (M1). Differences between beams with knots and beams without knots. Natural Frequency E[di f f erence] σ [di f f erence] Max di f f erence F1 F2 F3 -0,792 % -1,055 % -1,249 % 0,516 % 0,638 % 0,708 % -3,794 % -3,847 % -4,369 % from the knot presence. The differences between the timber beam without knots and the same beam with knots is shown in Tab. 9. The mean values and the standard deviations exhibit a small increase from the first to the third natural frequency and the maximum difference obtained is depicted in the third column for the three first natural frequencies. The natural vibration frequencies of the timber beam without knots and with uniform mechanical properties were calculated from the following expression: s n2 π ei (24) fn = 2 2L ρa where ρ is the mass density; a is the beam cross section; e is the Modulus of Elasticity (MOE); i is the second moment of area of the beam cross section (Tab. 7) and n is the number of the natural frequency. Due to the stochastic variation of the mechanical properties of the beam in each section, the first three natural modes of vibration are represented by three stochastic processes parametrized by the position in the beam span Φn (w, x). After the finite −0.692 1 −0.968 240 −0.641 2 2 3 120 0 −0.701 −0.671 Φ 0,75 m. −0.994 Φ 1,5 m. −0.988 1 110 0 −0.03 0 Φ 1,5 m. 0.03 2 120 210 105 0 −0.712 220 −0.702 Φ 2,25 m. −0.692 0.984 Φ 2,25 m. 1 −0.676 Φ 2,25 m. −0.65 1 110 0 0.968 240 2 120 3 240 −0.984 Φ 0,75 m. 0 −1 220 110 0 −1 105 p(Φ 1,5 m.) 220 −0.702 Φ 0,75 m. 210 p(Φ 1,5 m.) 0 −0.712 15 p(Φ 2,25 m.) p(Φ2 2,25 m.) p(Φ1 2,25 m.) 1 105 p(Φ 1,5 m.) 210 3 p(Φ 0,75 m.) p(Φ2 0,75 m.) p(Φ1 0,75 m.) VIBRATIONAL PROBLEMS OF TIMBER BEAMS 0 0.953 3 0.9765 Φ 1,5 m. 3 1 0 −0.702 3 Fig. 9 Histograms of the natural shape vibration modes at three points of the beam span using M1. Table 10 Mean values and standard deviation of the natural shape vibration modes at three points of the beam span using M1. Statistic x=0,75 m x=1,5 m x=2,25 m E[Φ1 ] σ [ Φ1 ] E[Φ2 ] σ [ Φ2 ] E[Φ3 ] σ [ Φ3 ] -0,701 0,0022 -0,981 0,0028 -0,681 0,0061 -0,991 0,0011 8,44 e-04 0,0056 0,965 0,0036 -0,701 0,0021 0,981 0,0029 -0,683 0,0053 element calculation is performed, the shape vibration modes are obtained for all the MC realizations. Then, a statistical analysis permits to obtain mean of the modes and other quantities such as the histograms. In Fig. 9, the stochastic process Φn (w, x) is presented for the first three natural frequencies through the histograms of Φn at three different points along the beam span. In Tab. 10, the mean values and standard deviations of each histogram of Fig. 9 are presented. As can be observed, the mean values and standard deviations vary with the mode order and the point of the beam at which they have been obtained. This proposed model M1 considers a variability different to the usual approach employed in the design practice. Although small, some statistical variations of the frequencies are observed. 5.2 Second Model M2 In the second model, timber knots are modeled as in M1. Additionally, the MOE in the zones without knots of the beam and in the zones with knots are modeled with the PDF presented in Section 4, (Tab. 4). The condition imposed to the values of the random variables that represent the MOEFKS and MOESW K is MOEFKS > MOESW K . Furthermore, the values of the MOESW K are independent random variables. The length of the weak zone was assumed equal to seven times the greater 16 Diego A. Garcı́a et al. 0.6 F1 0.5 F 2 F e% 0.4 3 0.3 0.2 0.1 0 400 600 800 1000 1200 ns 1400 1600 1800 2000 Fig. 10 Convergence of the mean values E[Fn ] for the first three natural frequencies (Model M2). 0.18 0.045 0.02 0.16 0.04 0.018 0.14 0.035 0.12 0.03 0.1 0.025 0.016 0.014 0.08 f[F3] f[F2] f[F1] 0.012 0.02 0.01 0.008 0.06 0.015 0.04 0.01 0.02 0.005 0 30 40 50 F1 Hz 60 0 100 0.006 0.004 0.002 150 200 250 F2 Hz 0 300 350 400 F3 Hz 450 500 Fig. 11 PDF of the first three natural frequencies found with M2. dimension of the knot, similarly to M1. The variability of the mass density between timber beams is introduced by the PDF presented before (Eq. 22). The simulation data is depicted in Tab. 7. A convergence study is shown in Fig. 10 where N is the number of independent Monte Carlo simulations that gives a prescribed accuracy. The criterion of convergence adopted in this model is the same as in M1, Eq. (23). An acceptable convergence is observed after N=1800. Features analogous to M1 are viewed in Fig. 11 (C f . Fig. 8), though the skewness is slightly different. In Tab. 11, the mean values, standard deviations and coefficients of variation of f[Fn ] are presented. Again, the behavior of the mean and standard deviation is similar to M1. A distinct feature is that the standard deviation and the coefficient of variation increase with respect to M1. This could be justified by the inclusion of MOE and density uncertainties in M2. In Fig. 12, the stochastic process Φn (w, x) is presented for the first three natural frequencies through the histograms of Φn at three points of the beam span. In comparison with M1 (C f . Fig. 9) the shapes of the histograms are different, influenced by the uncertainties of the MOE and the mass density introduced in M2. A larger dispersion is apparent. VIBRATIONAL PROBLEMS OF TIMBER BEAMS 17 Table 11 Natural frequencies (M2). Mean values, standard deviations and coefficients of variations of f[Fn ] F1 F2 F3 41,658 Hz 166,678 Hz 375,065 Hz 2,425 Hz 9,772 Hz 21,782 Hz 0,058 0,058 0,058 −0.542 340 −0.764 0 −0.722 340 −0.632 Φ1 0,75 m. 3 2 170 0 −1 240 −0.882 Φ 0,75 m. 2 120 0 −0.721 −0.631 Φ 0,75 m. −0.541 110 0 −1 −0.891 Φ1 1,5 m. −0.782 170 2 110 p(Φ3 1,5 m.) 220 220 σ [Fn ] E[Fn ] p(Φ 2,25 m.) p(Φ 2,25 m.) p(Φ1 2,25 m.) 3 δ [Fn ] = p(Φ1 1,5 m.) σ [Fn ] p(Φ2 1,5 m.) E[Fn ] p(Φ 0,75 m.) p(Φ 0,75 m.) p(Φ1 0,75 m.) Natural Frequency 0 −0.08 240 0 Φ 1,5 m. 0.08 2 120 0 0.77 0.885 Φ 1,5 m. 3 3 1 220 110 0 −0.705 340 −0.629 Φ1 2,25 m. −0.553 0.8845 Φ 2,25 m. 1 −0.621 Φ 2,25 m. −0.521 170 0 0.769 240 2 120 0 −0.721 3 Fig. 12 Histograms of the natural shape vibration modes at three points of the beam span using M2. Table 12 Mean values and standard deviation of the natural shape vibration modes at three points of the beam span using M2. Statistic x=0,75 m x=1,5 m x=2,25 m E[Φ1 ] σ [ Φ1 ] E[Φ2 ] σ [ Φ2 ] E[Φ3 ] σ [ Φ3 ] -0,633 0,0221 -0,880 0,0310 -0,627 0,0251 -0,895 0,0309 -1,85 e-04 0,0127 0,886 0,0314 -0,663 0,0226 0,881 0,0311 -0,626 0,0252 The mean values and standard deviation of each histogram corresponding to Fig. 12 are presented in Tab. 12. Similarly to M1, the mean values and standard deviations vary with the mode order and the point of the beam. 5.3 Numerical Simulation of Experimental Test In this section, a comparison between numerical and experimental results of the first natural frequency is presented. Experimental data, obtained by means of the test represented in Fig. 13, is reported by by Piter [15]. 50 sawn beams of Argentinean Eucalyptus grandis, with a nominal section of 50x150mm and length of 3m were employed. The fundamental frequency of vibration was obtained mechanically exciting the beams through an impact at one end and placing the sensor at the center of the body, in the anti-nodal position. A piezo-electric accelerometer type Vibrator PZ-10, 18 Diego A. Garcı́a et al. Table 13 Classification of tested beams. Strength class Number of tested beams C1 C2 C3 12 3 35 Fig. 13 Experimental test carried out to determine the first natural frequency. a oscilloscope type Fluke 123 Scopemeter 20 MHz and a software that permits the identification of the fundamental frequency through a Fourier transform of the harmonic spectrum, were employed in the test. The tested beams belong to the strength classes presented in Tab. 13, according to IRAM 9662-2 [11]. A numerical study was carried out. The timber beams were discretized with 100 beams elements and 1800 independent Monte Carlo simulations. In Fig. 14, a comparison between numerical and experimental PDF and cumulative distribution function (CDF) of the first natural frequency is shown. The results correspond to different weak-zone lengths. The numerically found CDF differs with the experimental curve in the lower part (under 50 %) of the plot. Also, a good prediction of the upper percentile values of the experimental CDF of the F1 is obtained with the length of the weak zone equal to seven times the knot dimension. In Tab. 14, the numerical and experimental results for different lengths of the weak zone are presented and compared. These results show a good prediction of the mean value of the first natural frequency but the standard deviations present a lower value in the numerical models. In order to find out the source of this disagreement, the 95 % confidence intervals of both the experimental sample results and the numerical realization results are derived for the mean value and the standard deviation of F1 (Tab. 15). As can observed, only the two cases with the larger weak-zone length (9 and 11 knot dimension) fall within the interval of the mean value of the experimental sample. Furthermore, only the last case contains the experimental mean value. The length of the interval is wider in the experimental case, probably due to the smaller number of samples and a larger standard deviation value. If a close inspection of the experimental samples is carried out, one is able to observe that the cases with larger number of defects are found in the lower part of the CDF plot. Due to the simplicity of the numerical model herein presented, some of these defects (e.g. presence of pith in 26 of the 35 C3, Tab. 13) are not taken into account giving place to more discrepancies. The model could be improved to consider other type of defects such as pith and improvements in the assessment of the material properties according to the knot size. The authors are at present working in this direction. As an illustration, some realizations of the lengthwise variation of the MOE are shown in Fig. 15. VIBRATIONAL PROBLEMS OF TIMBER BEAMS 19 0.12 1 Experimental 5 knot dimension 7 knot dimension 9 knot dimension 11 knot dimension 0.1 0.9 0.8 0.7 0.08 F(F1) f(F1) 0.6 0.06 0.5 0.4 0.04 0.3 Experimental 5 knot dimension 7 knot dimension 9 knot dimension 11 knot dimension 0.2 0.02 0.1 0 60 70 80 F1 Hz 90 0 60 100 70 80 F1 Hz 90 100 Fig. 14 Comparison between numerical and experimental PDF and CDF of F1 for different values of the weak-zone length. Table 14 First natural frequency. Comparison between numerical and experimental results for different values of the weak-zone length. Statistic Exp. 5 knot dim. 7 knot dim. 9 knot dim. 11 knot dim. E[F1 ] σ [F1 ] Min F1 Max F1 77,334 Hz 5,037 Hz 66,780 Hz 87,600 Hz 80,034 Hz 3,795 Hz 67,131 Hz 95,915 Hz 78,843 Hz 3,792 Hz 66,846 Hz 92,510 Hz 77,886 Hz 3,77 Hz 66,105 Hz 91,133 Hz 77,315 Hz 3,805 Hz 66,052 Hz 89,128 Hz Table 15 Confidence Intervals (CI) of the mean value and standard deviation of F1 . Data CI of µF1 Lower Limit Upper Limit CI of σF1 Lower Limit Upper Limit Exp. 5 knot dim. 7 knot dim. 9 knot dim. 11 knot dim. 75,912 Hz 79,858 Hz 78,667 Hz 77,691 Hz 77,139 Hz 4,201 Hz 3,675 Hz 3,672 Hz 3,651 Hz 3,684 Hz 78,775 Hz 80,209 Hz 79,018 Hz 78,040 Hz 77,491 Hz 6,276 Hz 3,923 Hz 3,92 Hz 3,897 Hz 3,933 Hz 10 1.6 x 10 MOE Pa. 1.4 1.2 1 0.8 0 500 1000 1500 Length mm. 2000 Fig. 15 Realizations of the lengthwise variation of the MOE. 2500 3000 20 Diego A. Garcı́a et al. 6 CONCLUSIONS The Probability Density Functions (PDFs) of the first three natural frequencies of a timber beam with uncertain properties are obtained with numerical simulations. Also, histograms of the mode shapes at certain points of the beam along its length are reported. The stochastic analysis allows to obtain more information of the dynamic behavior of the structural component. The influence of the timber knots in the response is frequently disregarded. In the present study, its consideration derives in an improved representation of sawn timber structures. Two stochastic models to account for the presence of knots are proposed. The first model introduces a random variable for the second moment of area along the beam span (I(x)). This variable considers the geometric parameters of the knots. The Modulus of Elasticity (MOE) of the sections with and without knots, and the mass density are assumed deterministic. Meanwhile, in the second stochastic model, also the MOE and the mass density are considered with uncertainties. In the first model, the difference between the timber beam with and without knots was quantified. The response found with this numerical model shows the variability of the response due to the knot presence for a simple beam with mean values of the MOE and the mass density. In the design practice, the influence of the timber knots is only taken into account in the selection of the MOE value and not in the reduction of the beam cross section that affects the bending stiffness. The results of the second model were presented, discussed and the differences with respect to the first model were assessed. The PDFs of the first three natural frequencies found with the second numerical model vary the shape, the mean and the standard deviation compared with the first one. On the other hand, the histograms of the natural modes show that the variation in the mean value is small between the two models while that the standard deviation increases due to the greater uncertainty present in the second numerical model. The model of the lengthwise variability of the MOE was stated starting from the weak-zone approach proposed by others authors to study the bending strength. However, the model herein presented introduces the presence of knots in the sectional parameters and the length of the weak zone in a different way. Results of numerical simulations of experimental tests carried out to determine the first natural frequency were also presented. They show that the numerical simulations provide results relatively close to the ones obtained experimentally. The difference between numerical and experimental results might be influenced by timber defects which the numerical model does not include, such the presence of pith. The mean and upper percentile values of the numerically obtained CDF are close to the experimental CDF. Thus, the agreement is better when the timber beams are of superior quality. The stochastic models presented in the present study constitute a more realistic material approach, feasible to be applied to reliability studies of serviceability limit states of structural components made of Eucalyptus grandis timber. VIBRATIONAL PROBLEMS OF TIMBER BEAMS 21 7 ACKNOWLEDGEMENTS The authors acknowledge the financial support of CONICET, SGCyT-UNS, MINCyT from Argentina and CNPq, and FAPERJ from Brazil. The authors also acknowledge to Dr. J.C. Piter for the data employed in this work. References 1. Ang, A.H.S. and Tang, W.H., Probability Concepts in Engineering: Emphasis on Applications to Civil and Environmental Engineering. John Wiley and Sons, Inc., (2007). 2. Baño V., Arriaga F., Soilán A. and Guaita M., Prediction of bending load capacity of timber beams using a finite element method simulation of knots and grain deviation, Biosystems engineering, 109(4):241-249, (2011). 3. Baño V., Arguelles-Bustillo R., Regueira R. and Guaita M., Determinación de la curva tensióndeformación en madera de Pinus sylvestris L. para la simulación numérica de vigas de madera libre de defectos, Materiales de construcción, 62(306):269-284, (2012). 4. Baño V., Arriaga F. and Guaita M., Determination of the influence of size and position of knots on load capacity and stress distribution in timber beams of Pinus sylvestris using finite element model, Biosystems engineering, 114(3):214-222, (2013). 5. Bathe, K.J., Finite element procedures. Pentice hall Englewood Cliffs, volume 1, (1996). 6. Czmoch, I., Influence of structural timber variability on reliability and damage tolerance of timber beams. Ph. D. thesis, Luleå tekniska universitet, (1998). 7. Escalante M.R., Rougier V.C., Sampaio R., and Rosales M.B., Buckling of wood columns with uncertain properties, Mecánica Computacional, Volume XXXI, Number 14. Uncertainty and Stochastic Modeling, pages 2735-2744, (2012). 8. Faber, M.H., Khler, J. and Sorensen, J.D., Probabilistic modeling of graded timber material properties, Structural Safety, 26(3):295-309, (2004). 9. Guindos P. and Guatia M., A three-dimensional wood material model to simulate the behavior of wood with any type of knot at the macro-scale, Wood Science and Technology, 47(3):585-599, (2013). 10. Guindos P. and Guatia M., The analytical influence of all types of knots on bending, Wood Science and Technology, 48(3):533-552, (2014). 11. IRAM:9662-2, Madera laminada encolada estructural. Clasificación visual de las tablas por resistencia. Parte 1: Tablas de Eucalyptus grandis. Instituto Argentino de Racionalización de Materiales IRAM, Buenos Aires, (2006). 12. Isaksson, T., Modeling the Variability of Bending Strength in Structural Timber: Length and Load Configuration Effects. División of the Structural Engineering, Lund University, (1999). 13. Köhler, J., Reliability of timber structures. Vdf Hochschulverlag AG, 301, (2007). 14. Köhler, J., Sørensen, J. D. and Faber, M. H., Probabilistic modeling of timber structures, Structural safety, 29(4):255–267, (2007). 15. Piter, J., Clasificación por resistencia de la madera aserrada como material estructural. Desarrollo de un método para el Eucalyptus grandis de Argentina. Ph. D. thesis, Universidad Nacional de la Plata, (2003). 16. Rubinstein, R. Y., Simulation and the Monte Carlo method. John Wiley and Sons, Inc., (1981). 17. Shannon, C., A mathematical theory of communication, The Bell technical journal, 27:379–423, (1948). 18. UNE-EN 408, Estructuras de madera. Madera aserrada y madera laminada encolada para uso estructural. Determinación de algunas propiedades fı́sicas y mecánicas. Estructuras de Madera, AENOR Asociación Espaola de Normalización y Certificación, Madrid, (1995).