Ductility analysis of laminated timber beams of small section height

The purpose of this paper is to analyze and determine the applicability of the models that take into account the ductility of timber elements in bending. Beam elements whose span is many times higher than the height are experimentaly tested. In addition to the tested beam elements, paper presents the mechanical properties of the timber material and the results obtained by finite element method (numerical models). Numerical models are made in the software package Abaqus.


Introduction
Ductility is a desirable mechanical property because ductile structures have a certain "reserve" of bearing capacity, i.e. relatively high plastic deformations occur prior to failure.In addition, ductility is desirable from the standpoint of robustness and reliability.In that respect, it is generally considered that ductile structures may be more reliable and robust compared to brittle structures [1][2][3].While the analysis according to plastic theory is frequent in steel structures, the assumption of elastic behaviour until failure is always applied in the calculation of carrying capacity of timber structures.The real behaviour of wood subjected to bending is highly complex because wood has different values of compressive and tensile strength in the direction parallel to fibres (Figure 1).This complex behaviour is due to different values of elastic modulus, and to differences in constitutive laws for compression and tension in the direction parallel to fibres.If we consider an element subjected to bending moment only (without stability problems), which causes stress that is much lower that the tensile and compressive strength values (Figure 2.a), then the distribution of stress in transverse direction parallel to fibres is linear, and the neutral axis passes through the cross-section centre of gravity (traditional engineering approach).In case the real behaviour is considered (bilinear constitutive law) (Figure 2.b), the neutral axis is situated under the cross-section centre of gravity.If bending moment is further increased, the compressive strength parallel to fibres is achieved in edge fibres (Figure 2.c).Once the compressive strength parallel to fibres is achieved, the material starts to plasticise in compression zone, with a simultaneous increase in tensile stress.Finally, (Figure 2.d), once the tensile strength parallel to fibres is achieved, fibres yield in the tensile zone, and the element loses its cross-sectional resistance.

Current state-of-the-art
The bending moment M el at which the cross-sectional failure occurs according to the traditional elastic theory amounts to: where σ is the design bending strength, b is the cross-sectional width, and h is the cross-sectional height.As a rule, the design bending moment M el is smaller than the real carrying capacity because the design does not take into account the nonlinear stress distribution which is obvious because of the fact that elastic moduli and constitutive laws at compression and tension parallel to fibres are different.This behaviour is complex and so various authors [4][5][6][7][8][9] have assumed different stress distribution values along the element's cross-section.Ductility analysis of laminated timber beams of small section height The common feature of the examined models is the existence of parameter n, which is defined as follows: where f t || is the tensile strength parallel to fibres, while f c || is the compressive strength parallel to fibres.According to [4] the bending resistance moment M u is determined according to According to Moe [4], compressive stress values are constant up to a certain level, after which a stress jump occurs.After this jump, stress values are linearly variable, and the neutral axis coincides with the cross-sectional centre of gravity.The model depends on coefficients c and α whose values are not strictly specified in literature.
In the paper presented by Nwkoye [5], it is assumed that compressive stress values are constant up to a certain level, after which the stress curve is linearly variable.The neutral axis does not pass through the cross-sectional centre of gravity.
Zakić [6] assumes that the distribution of compressive stress values can be approximated with the second-order curve, while tensile stress values are linearly variable.In addition, one of the basic assumptions is that the neutral axis does not pass through the cross-sectional centre of gravity.Then the bending resistance is calculated according to: Bazan [7] assumes a different distribution of stress values.
According to this assumption, the compressive and tensile stresses are both linearly valuable, and the following expression is used for bending moment:  Dean Čizmar, Domagoj Damjanović , Krunoslav Pavković, Vlatka Rajčić Zaw et al [8] propose the distribution of compressive stress values that is constant up to a certain height of cross section, after which it corresponds to the second-order curve.Tensile stresses are linearly variable, and the neutral axis moves toward to the tensile edge.The following expression is used for bending moment: Buchanan [9] proposes the model similar to that given by Bazan [7], but introduces the coefficient c that is dependent not only on the relationship between the tensile and compressive strength values, but also on variable, i.e. decreasing elastic modulus in the compressive zone, which is described with the coefficient m.
Due to poor reliability and the lack of information about analysis of variable elastic modulus, this model is not analysed in this paper.Four studied models, differing by the value of n and the compressive to tensile strength relationship, are presented in Figure 5.The value M u is the plastic resistance model, while M el is the resistance moment according to elastic theory.Models presented by Bazan [7] and Nwkoye [4] are also valid for the case when n = 1, while other models are valid only in cases when n > 1. Models presented by Zakić [6] and Bazan [7] provide smaller resistance moments compared to other two models.All studied models are presented in Figure 5.

Introduction
Experimental tests were conducted at the Structure Testing Laboratory of the Faculty of Civil Engineering, University of Zagreb, in order to analyse the existing models and to define, based on this testing, ductility level for wood as construction material.The testing was organised in such a way that three large-scale samples of timber class GL24h (420 cm in total length) were exposed to bending.In addition to these tests, which form the central part of the testing campaign, some small-size samples made of the base material were also tested (Figure 6).Small-size samples were tested to bending (10 samples), compression perpendicular to fibres (9 samples), compression parallel to fibres (9 samples), tension parallel to fibres (9 samples), and shear (9 samples).

Base material testing
Results obtained by small-size sample testing were statistically analysed and modified using appropriate coefficients to take into account: timber errors, timber moisture, influence of sample size (size effect factor), and load increase rate.The following average timber strength values were obtained through the analysis of results: tensile strength parallel to fibres N/mm 2 , shear strength f sh = 2 83   , N/mm 2 , and bending strength f m = 37 68 , N/ mm 2 .In addition, the modulus of elasticity parallel to fibres of E L =12740 N/mm 2 was obtained by subjecting small-size samples to bending action.

Preparations for experimental testing
The experimental testing was conducted with three samples class GL24h, 600 cm in length (420 cm for bending test, and the rest for preparation of base-material samples), and crosssectional dimensions amounted to b/h = 10/12 cm.The objective of testing the beam model with two spans exposed to bending was to analyse timber behaviour in the support zone and in the span.The timber class was defined according to HRN EN 338:2006, and the testing was conducted according to EN 408:2006 [10].The diagram of the tested girder is presented in Figure 7.The moisture of each girder was measured at several positions prior to testing so as to check possible influences on mechanical properties.The moisture was tested at: left-side girder support, span between the left-side support and central support, central girder support, span between the central support and right-side support, and right-side girder support.Positions for measuring displacement and strains are shown in Figure 7.The measurement position marked with the ordinal number I. represents vertical displacement at the first girder support measured in order to check boundary conditions.Measurement positions II.and III.represent vertical girder displacement in halves of the span.Measurement sections marked with letters A-A and B-B represent sections in which girder strains were measured.Strain in compression and tensile zones was measured in each section.Thus, for section A-A, strains in the top zone were marked with numbers 1 and 2, and with numbers 3 and 4 in the bottom zone.In section B-B, points for measuring strain were marked in the top zone with numbers 5 and 6, and with numbers 7 and 8 in the bottom zone.It should be noted that the experimental testing was conducted in accordance with HRN EN 380:2006 [11].All samples were tested until failure.However, they were first subjected to 40 % of the failure force obtained in previous analyses.Once the 40 % of failure force was achieved, the load was released.In the next cycle, the sample was once again subjected to 40 % of the failure force, and this force was maintained for 60 seconds.After that, the load was gradually increased until failure.The load increase rate was related to the deflection level (displacement check) and was set to 0.2 mm/min.The girder position immediately before the testing is presented in Figure 8.The finite element method was used to analyse the model with geometrical dimensions corresponding to those of large-scale samples used in experimental testing.This analysis, based on the finite-element method, was conducted using the program package Abaqus/CAE, Ver. 10 with the UMAT subroutine.Mechanical strength properties of timber were taken from experimental testing fro smallsize samples.The tensile strength perpendicular to fibres was taken from the available literature [12], f t ⊥ = 0 38 , N/mm 2 , as this analysis was not made in the scope of experimental testing.The testing conducted according to paper [12] involved wood that grew under similar climatic conditions, and the quality and type of that wood corresponds to the wood used in this testing.The obtained elastic modulus of E L = 12740 N/mm 2 , and the bulk density of wood determined through experimental testing, were used as basis for adoption of elastic modulus, shear, and Poisson's ratio from the available literature [13][14][15][16][17][18].Mechanical properties for each direction were obtained as average values of E L = 12563,50 N/mm 2 , E R = 902,90 N/mm 2 , E T = 542,91 N/mm 2 , G LR = 742,68 N/mm 2 , G LT = 660,95 N/mm 2 , G RT = 67,32 N/mm 2 , ν LR = 0,41, ν LT = 0,505, ν RT = 0,495.Stress strain diagram for timber, described with the UMAT subroutine, is presented for normal stress values in Figure 9.

Boundary conditions, geometry and analysis of a beam element using the finite-element method
The timber beam element and steel plates were modelled using finite elements defined with twenty points.Previously defined mechanical properties of wood were implemented, together with the derivation of the Hill criterion, via UMAT subroutine, into the Abaqus software.The UMAT subroutine Ductility analysis of laminated timber beams of small section height (with the yield criterion derived from the Hill criterion), which provides a good-quality description of wood failure modes, is described in paper [19].The beam element model was defined in such a way that each slice 4 cm in thickness was modelled separately.Four modelled slices were interconnected by cohesive interaction between contact surfaces in such a way to form a beam element 20 cm x 12 cm in cross section.In addition, the bottom and top slices were divided into segments, and these segments were connected to one another via cohesive interaction.The described complex numerical model was created in order to obtain complex yield modes for the beam element, where the yield is due to: tensile stress, compressive stress, shear stress, and delamination.Beam element boundary conditions were modelled using steel plates with dimensions identical to those used in experimental testing (20 cm x 10 cm).Plates were modelled as being ideally elastic, and the elastic modulus amounted to E = 210000 N/mm 2 and the Poisson's ratio to ν = 0.3.The numerical model of beam element with boundary conditions is presented in Figure 10.An absolute stiffness for normal stress, and the possibility of tangential sliding of two surfaces with the friction coefficient, were used for modelling contact surfaces between steel plates and the wooden element.The coefficient of friction used for the wood-steel contact was taken from available literature [20] and amounts to µ = 0,25.The analysis of the model was conducted using the nonlinear analysis, including the geometrical and material nonlinearity.The Newton Method with automatic control of increase rate was adopted for displacement checking.The maximum displacement increase increment was limited to 0.1 mm, and the initial increment was set to 0.02 mm.

Cohesive interaction between two surfaces
The cohesive interaction between contact surfaces was used, together with the UMAT subroutine, to define the behaviour and yield of wood.The formulation of cohesive interaction is very similar to cohesive elements with the crack opening and spreading possibilities.The cohesive interaction is defined with the stiffness coefficient E i 0 (i = 1, 2, 3, which is related to normal directions, and to two shear directions), the cohesive interaction strength σ i 0 and the energy needed to achieve the yield of cohesive interaction G i 0 .
Figure 11 shows the opening and yield behaviour for cohesive interaction, which consists of elastic deformation, start of yield, and linear stiffness degradation of cohesive surface.The start of degradation of cohesive surface can be presented as δ σ In this paper, the cohesive stiffness coefficient was adopted as amounting to E i 0 = 10.000N/mm and it is not directly related to mechanical properties of wood.The cohesive interaction yield starts once the yield criterion has been met.For the problem under study, the criterion of the sum of stress values squared has been adopted, and it can be written as follows: where σ 1 means that the yield can occur in case of tensile stress only.The yield in cohesive interaction for one of main directions can be expressed by failure energy that is equal to the area under the curve presented in Figure 11.By using the failure energy, it can be written that the total deformation at the moment of full yield is equal to . The yield in cohesive interaction is the degradation of stiffness coefficient, and can be written using the damage coefficient: where δ m max is the total maximum deformation for all directions, δ m 0 is the deformation at the start of yield, and δ m k is the total deformation at the total yield of cohesive interaction.Constitutive equations for the complex crack opening state are defined as follows:

Results of experimental research
All girders were subjected to load amounting to 40 % of the estimated force of failure (i.e.20 kN).Three diagrams presented in Figures 12, 13 and 14 show the measured values of strain in the first span and near the left-side girder support (according to Figure 7 and Table 1).In addition to measured deformation values, mean values for individual zones were also calculated and presented.The diagrams show that, in the tensile zone, all three girders behave elastically until failure, which could have been expected.In the compression zone, the plasticisation in span and at one measurement point at support was registered for sample GR3.The force of failure for sample GR1 amounted to 56.94 kN.Once this value was achieved the force dropped suddenly and the sample broke at the sawtooth connection of the slice, and the crack spread to the connection zone of the bottom slice in the first span of the girder, in the tensile zone (Figures 18  and 19).The bending stress at failure amounted to 44.9 MPa, which greatly exceeds the medium standard bending strength value (calculated from typical value given in tabular form).
On the other hand, the shear stress at the force application point amounts to about 84 % of the mean standard shear strength (calculated from typical value given in tabular form), which points to the insufficient shear strength of this sample.Ductility analysis of laminated timber beams of small section height It can be concluded by analysis of results obtained by the finite-element method that the girder elements most often yield due to bending moment in the central support zone, and in the middle of two spans.In maximum moment zones, slices yield due to tensile stress parallel to fibres.The failure occurred due to tensile stress and this without a ductile character.The initial failure at the top slice on the support, i.e. at the bottom slice if the failure occurs in mid-span, results in delamination or separation of slices, as shown in Figure 27.Similar failure modes were also observed during experimental testing, which is why it can be concluded that the girder elements subjected to testing are described quite well with numerical models.Another beam model (marked as numerical model II) was made simultaneously with the beam model whose mechanical properties of wood were obtained by experimental testing on small samples.The tensile strength of this second mode was two times greater that the compressive strength parallel to fibres: 50 N/mm 2 .The purpose of this model was to enable comparison of ductility, and to verify theoretical expression for ductility of wooden beam elements.The results of this model with higher tensile strength are presented in Figure 24.The results obtained show that there is a small ductility D = 1.30obtained according to paper [19], which arises from deformation in the compression zone, but these values are much smaller than those obtained through theoretical expressions.Furthermore, it should be noted that the shear strength was also increased by two times in this numerical model ( f sh = 5 66 , N/mm 2 ), as it was determined that the model fails due to shear stress, with total delamination of the beam element (brittle fracture).

Conclusion
The purpose of the paper is to analyse and determine applicability of models that take into account ductility of wood subjected to bending action.The results obtained by experimental tests point to the fact that the bending strength exceeds the standard value, and that the failure of samples is caused by shear stress (2 samples).The failure of the third sample occurs in the force application zone and is due to bending.Failure modes obtained by experimental testing point to the fact that there is a very small difference in the failure force values and the corresponding modes of failure.Results obtained point to elastic behaviour until failure.Comparison of results obtained in this research with the data presented in [4][5][6][7][8][9] shows that the latter data do not realistically describe failure of wooden beam elements.Models analysed by the finite-element method also point to stiff behaviour.The numerical model I points to the lack of ductility.The minimum ductility (D = 1.3) is registered at model II, i.e. in case of girders made of a better-quality wood where the tensile strength parallel to fibres is two times greater than the compressive strength parallel to fibres.This model is also characterized by a very-high shear strength.It should be stressed that the brittle failure by delamination due to shear stress has also been registered during analysis of models with high tensile strength and with unchanged shear strength of wood.Consequently, it can be concluded that the models that take into account ductility at bending are not applicable, i.e. that the ductility at bending is negligible.At the same time, there is a great possibility that the girder will fail due to displacement, which undoubtedly results in brittle failure.

Figure 1 .
Figure 1.Idealised stress-strain diagram for wood, concrete and steel

Figure 2 .
Figure 2. Behaviour of timber element exposed to pure bending: a) traditional engineering approach; b) bilinear constitutive law; c) bilinear constitutive law -plasticisation of compression zone; d) bilinear constitutive law -cross-section failure in tensile zone

Figure 5 .
Figure 5. M u /M el relationship as a function of tensile to compressive strength ratio for wood parallel with fibres

Figure 9 .
Figure 9. Stress strain diagram for timber, defined with UMAT subroutine

Figure 8 .
Figure 8. Girder immediately prior to testing

Figure 10 .
Figure 10.Beam element model with boundary conditions

Figure 12 .
Figure 12.Force -strain diagram, sample GR1 (measurement positions: 1 and 2 top zone in span, 3 and 4 bottom zone in span, top zone at support 5 and 6 and bottom zone at support 7 and 8) Vertical displacements at measurement positions I, II, III and IV are presented in Figures 15, 16 and 17.It can be concluded from these diagrams that the relationship between force and displacement (measurement positions II and III) is almost linear until failure.

Figure 13 .and 8 )Figure 14 .
Figure 13.Force -strain diagram, sample GR2 (measurement positions: 1 and 2 top zone in span, 3 and 4 bottom zone in span, top zone at support 5 and 6 and bottom zone at support 7 and 8)