CREEP IN THE FUNDAMENTAL FREQUENCY AND STABILITY OF A SLENDER WOODEN COLUMN OF COMPOSITE SECTION

Creep is a phenomenon that can occur in wooden structures since wood is a viscoelastic material. Creep may change the purely elastic parameters determined in wood characterization initial tests, as its behavior depends on the rheology of the material, even under a constant stress level. Mathematically, creep can be characterized by models in which the immediate elastic deformation is increased by a viscous deformation, resulting in a temporal function. For this reason, the calculation of the natural frequency of vibration and the stability verification of a slender column should include the reducing effects of stiffness both of axial force and creep. The first one can be considered through the geometrical portion and the second one by the introduction, in the conventional portion, of a variable elasticity modulus over time, obtained in relation to the adopted rheological model. A numerical simulation was performed to evaluate the aspects above, considering a bar compressed by a force at the free end equivalent to 10% of the Euler critical force, plus its own weight, adopting a rheological model with three parameters for the variation of the elasticity modulus. The results show differences of 60% and 50% for the frequency and elasticity modulus, besides defining the exact instant of column collapse in the case of its non-observance.


INTRODUCTION
The slow increase of deformation under constant stress over time is called creep.Mathematically, creep can be represented by a time-dependent function associated with viscoelastic rheological models capable of describing the phenomenon, according to Findley; Lai; Onaran (1989).The slow deformation, or creep, in wooden structures, is a time-dependent phenomenon, which is also related to loads and deformations, defined as the increase of deformation under the action of loads, or permanent stresses, over time.Creep produces residual deformation and can change the characteristics of the materials and their mechanical properties, and it is even able to produce the failure of a structure, as reported by Gottron;Harries;Xu (2014).When wooden structures are used as horizontal components and loaded, they undergo long-term deformation and, in different humidity conditions, the wood has a tendency to increase its deformation due to creep, which can lead to operational problems due to excessive deformation, or may cause safety problems due to loss of resistance, which can decrease the load capacity and, in the case of columns, reach the buckling stress, emphasize Epmeier; Johansson; Kliger; Westin (2007).
In general, two types of criteria are used to represent the creep within the structural design.The first one is related to a deformation limit value, restricting it to validity intervals, and the second one, less common, takes into account a monotonic growth of deformation.From a practical point of view, the technical standards take into account the creep phenomenon in the design of structures by two considerations: by proposing a coefficient of increase or decrease of stiffness, or proposing a coefficient of increase or decrease in resistance as the loading acting time and humidity class, according to Celia-Silva;Calil Junior (1992).The decrease in stiffness is related to the fact that the axial force is responsible for this effect and can lead to loss of stability of a column.In this sense, Gambhir (2004) refers to the influence of the compressive force highlighting that once there is a reduction in the structural stiffness, the effect of the loads on the structure increases, which also increases the forces on the elements, so the structure resistance capacity decreases.For this reason, in the case of compressed columns, a premature analysis can produce undesirable consequences, and any failure causes catastrophic effects since it involves the balance of structures, Timoshenko (1992).
Regarding the study of wood mechanical properties, it is a typically elastic or elastoplastic material, depending on the level to which it is subjected.Its behavior, concerning the elasticity, is considered as anisotropic, i.e., a material whose elastic properties vary with the considered direction, differently from an isotropic material, which keeps its properties constant regardless of direction.During tree growth, the internal structure of the solid material which consists of wood becomes highly oriented; for this reason, it is an anisotropic material.On a macroscopic level, the arrangement of the wood structure results from both length and diametral growth processes.Therefore, wood structure is usually referred to two privileged axes: the longitudinal direction (L) along the fiber and on the radial direction (R) in relation to the annual growth rings.The direction (T) is tangential to the annual growth rings.The three orthogonal planes (RT), (RL) and (LT) are then three planes of symmetry to the internal structure of wood, as explained by Oudjehane;Raclin (1995).Wood anisotropy also results from its cell structure, which is most composed of fibrous material, which will influence its apparent mechanical properties.The anisotropy changes the vibrational properties of viscoelastic materials, i.e., it modifies the structural dynamics since it depends on the elastic or Young's modulus (E 0 ) and shear (G 0 ), changing the standard of modes and vibration frequencies, said by Brémaud; Gril; Thibaut (2011).
Wood anisotropy is well represented by the mentioned orthotropic model, which takes the longitudinal, radial and tangential directions, with elastic parameters for each of them.That means there is a modulus of elasticity and poison coefficient for each of these directions, making impossible to operate the known linear elastic relationship, specific to isotropic material, between these parameters to determine the transversal elasticity modulus.Accordingly, creep, given the considerations presented in this study, is characterized by specific deformation analyzed in the longitudinal direction of the column, which varies with time, and whose elastic reference is the elasticity modulus parallel to fibers.In line with this approach, Schniewind; Barrett (1972) state that wood can be considered as an orthotropic viscoelastic material.They observed that for certain purposes and under certain conditions, the analysis of stress on wood can be considered as linear elastic and viscoelastic and that the region of viscoelastic behavior is large and has considerable practical importance.A mathematical approach on creep which includes the three directions of wood orthotropy was developed in the sixties of the twentieth century (BHATNAGAI; GUPTA, 1967).
Because of its plastic behavior, even with a constant level of tension, the deformations in structural elements of wood tend to increase with time, i.e. after the initial elastic deformations, additional deformations occur, which can be reversible or not, but never complete.The slow creep is partially reversible.For an unloading, after elastic recovery, there is a later recovery, which is called slow recoverable deformation, slow reversible deformation, or delayed elastic deformation, and only a remaining portion of the deformation is residual or irreversible, being this portion of deformation called creep, Leohard;Mong (1977).
In Brazil, there are some records of studies regarding creep on wooden structures carried out by Almeida (1990) and Celia-Silva;Calil Junior (1992).The studies conducted by the first researcher, on tensioned parts, revealed high levels of creep after the first twentyfour hours of the assay, corresponding to a rapid creep phenomenon, and a convergence after 90 days.The measurements were made daily and involved the recording of deformations and environmental conditions (temperature and relative air humidity).On the other hand, the second researchers conducted experiments on slow deformation in wooden columns.Their studies evidenced a tendency for creep stability after 60 and 90 days of loading.Studies on creep in wood, considering it under several aspects, were more recently developed by George et al. (2003), Wang, J.B.;Foschi, R.O.;Lam, F (2012;2012), Reynolds, T.; Harris, R.; Chang, Wen-Shao (2013), Honfi,D. et al. (2014), Sharapov, E.;Mahnert, K.C;Militz, H. (2016) andSaifouni et al. (2016).
Commonly, the representation of creep is based on rheological models, which associate deformations deferred at the time.These models, according to Gomes et al (2007), allow to simulate the response of the material to forces and tensions applied to it.The inclusion of these rheological models can be made both in static analysis and dynamic of the structures.In the case of dynamic analysis, the structure stiffness must be composed of two terms; the first one corresponding to the portion of conventional stiffness and the second to the portion of geometric stiffness, according to Clough; Penzien (1993).Thus, it is possible to adapt the first portion of stiffness by introducing a variable elasticity modulus over time, making possible to monitor the increase of deformations, according to the adopted rheological model, keeping a constant stress level.Therefore, the total stiffness assumes the form introduced by the first portion, via elasticity modulus, in the rheological model intended to represent the creep; and the second portion is geometrical, which is a function of normal operating force, which includes the own weight of the structural element.
A numerical simulation was performed to evaluate the creep in natural frequency and in the stability of a slender wood column, taking into account a threeparameter model to represent its viscoelastic behavior.The formulation developed to calculate the column natural frequency of vibration was based on the method of virtual works and allowed the evaluation of the buckling critical load and loss of system instability, besides calculating the frequency of the first mode.The results indicated significant differences in the fundamental frequency and elastic modulus, considering values obtained by a purely elastic analysis and those nonlinear with the creep, in relation to the material, and geometric value considering the geometric stiffness.It was also possible to define the exact moment of the column collapse in the case of non-observance of creep, considering a force of 10% of the critical load of Euler concentrated at the free end of the structural element.Comparison with design criteria and verification for parts in flexo-compression and stability of the NBR 7190 -Wooden Structures Design (1996) were adopted for this level of loading, as well as the definition of the limit percentage of the Euler critical load, which leads the system to the collapse, using these same criteria and mathematical tools developed in the present study.

Mathematical simulation support
The rheological models used to represent the viscoelastic behavior of fluids and solids are usually based on the association of springs and dampers that predict total deformations and try to describe more appropriately the behavior of each material or group of materials, as presented in Figure 1.The model of Group I (Figure 1a), for example, represents bodies that have an initial elastic deformation and viscous deformation that occurs over time, so this model is more suitable for the description of certain solids.The model of Group II (Figure 1b) describes the behavior of liquids simulating an initial viscous outflow and a slow elastic deformation over time.The model of Group III (Figure 1c) shows an instantaneous elastic response followed by a viscous response with a slow deformation as well.The model of Group IV (Figure 1d) simultaneously presents two slow viscous deformations over time.These models can be used alone or in associations, forming chains, in an attempt to obtain the best representation for each situation.
One of the models more used to represent the creep in several materials is the three-parameter model, in which an elastic parameter E 0 is associated with a viscoelastic model of parameters E 1 and  1 called Kelvin-Voigt model, which is a simplification of the Burgers model of Group type I (Figure 1a).According to Keramat;Shirazi (2014), a single model of the Kelvin-Voigt is enough to describe the viscoelastic nature of many solid.Based on this consideration an adaptation of the Burgers model was used by Kränkel; Lowke; Gehlen (2015) to analyze the non-linear viscoelastic deformation in bonded anchors.In relation to the wood, Mukudai (1983) mentions that the functional form represented by the Voigt and Maxwell models is very efficient since they can be conveniently applied to assess the viscoelastic behavior of that material through mathematical models.In this sense Hering;Niemz (2012) have used the Kelvin-Voigt model, together with experimental activity to descript the creep in bended pieces.A Burger chain of four elements was used by Gao;Wang;Shao (2016) to study the behavior of Xuan paper, making possible to description convenientily the elastic, viscoelastic and plastic deformations presented by the material.
The adoption of Kelvin-Voigt model has repeatedly been made to study phenomena in various science fields.Hackney; Aifantis; Tangtrakarn; Shrivastava (2012) used this model to examine the creep behavior of nanostructures of composites.In the context of the finite element method, Chung; Tamma; Namburu (2000) used the Kelvin-Voigt model to represent the primary creep to investigated the thermo-viscoelastic of composite structures.Puertas; Gallego (2014) used the same model, together with the boundary element method, to consider the viscoelasticity in the solution of the three-dimensional dynamic problems.
Mathematically, the total deformation of the Kelvin-Voigt model is given by equation ( 1), where e e is the deformation in the elastic model and e v is the deformation in the plastic model.Take into account its derivative relative to time, the total deformation of the model is found by equation ( 2), Ec 0 ef Figura 1 -Viscoelastic rheological models.
Figura 1 -Modelo reológicos viscoelásticos. (1) Creep in the fundamental frequency and stability... which is the constitutive equation of the Kelvin-Voigt model.
The equation ( 3), for this specific case, is obtained by assuming From the previous equations, it is obtained the differential equation ( 4) for the Kelvin-Voigt model, with t representing the moment of load application.
As the stress is constant, the derivative of stress relative to time is null.By applying the above stress conditions, the equation ( 4) is reduced to equation ( 5), which is an ordinary differential equation, whose general solution, for t > 0, with the initial condition (0) = is given by equation ( 6).
Obviously, if the stress level remains constant, the elasticity modulus must simultaneously reduce to as the deformation increases, then: After finding the variable elasticity modulus over time, it is necessary to establish the conditions for an analysis of vibrations that, from the point of view of structural dynamics, make possible the obtainment of the temporal variation of the system natural frequency.One of the classic methods is the Modal Analysis.This method allows the obtainment of information about the behavior of a structural system and therefore reveals important aspects related to its dynamics, such as frequencies and natural modes of vibration, according to Carrion; Mesquita; Ansoni (2014).In addition, it is possible to verify the stability of a structural system and the moment of its collapse, based on the nullability conditions of the fundamental frequency and, consequently, its stiffness.
The formulation developed in this simulation to consider the creep in the vibration of a column is based on the principle of virtual works associated with the Rayleigh technique, which suggests that for a system containing infinite degrees of freedom can be associated with another system with a single degree of freedom to bring its frequency closer, as described by Leissa (2004).The process takes into account the suitable choice of the generalized coordinate that describes the structure deformation, considering the first vibration mode.In the end, the equation of motion appears regarding the generalized system properties such as stiffness and mass, which are necessary to calculate the frequency.It is important to note that the technique developed by Rayleigh aimed to calculate the fundamental frequency of vibration of elastic systems, and the precision obtained by this method depends directly on the function chosen to represent this vibration mode, as also said by Leissa (2004).
The basic concept of the method is the principle of conservation of energy, and can, therefore, apply to linear structures (or not), according to Clough; Penzien (1993).Temple;Bickley (1933) consider that the Rayleigh technique is applied both to systems with infinite degrees of freedom and continuous systems, and serves both to determine the fundamental period of vibration and to verify the stability of mechanical systems, but it is not limited to these purposes.As seen in studies carried out by Villasenor; Farfán; Guzmánb; Romero; Castellanos; Sesma (2014), the Rayleigh technique was used in the detection of subsurface cracks in solids.Pena (2015) studied the asymmetric vibration of rigid bodies, considering the damping of Rayleigh in the context of small rotations, allowing the mathematical simplification of the problem.Cámara and Astiz (2014) used the form function technique as recommended by Rayleigh in the dynamic modal analysis of cable-stayed bridges, within the region of elastic deformation of the material.
It is considered the model of the bar in Figure 2, of length L, in a free motion non-damping free motion, with the generalized coordinate q(t) defined on the upper end of the element.The function v(x,t) provides the horizontal displacements for each x position along the height, in relation to a function (x) that describes the form of the vibration and meets the boundary conditions of the problem.It is assumed this system is composed of a prismatic bar, made of a viscoelastic material, clamped on its the base, supporting, in addition to its own weight m, a mass m 0 at the free end, representative of bodies fixed to its top.N(x) is the normal generalized force and e(t) is the vertical displacement of the bar.
It is also considered that the system motion does not change the direction of the normal force.Thus, the represented structure is a bar in bending, so that the virtual work of internal forces W I is performed by the bending moment M(x,t), acting on the deflection of the virtual bar.It is assumed that the section remains flat after deformation.The Principle of Virtual Works requires that the virtual work of external forces is equal to the virtual work of the internal forces.The virtual work of the external forces is obtained by equation ( 8), in which the is f I (x) = m 1 (x) v(x,t) the inertial force.In turn, the virtual work of the internal forces is given by equation ( 9).
An infinitesimal element ds of the curved bar is required to find the axial displacement e(t).The shortening of the axis due to the axial displacement is given by equation ( 10).
From the binomial development, considering that the higher order terms are smaller compared to the first order term, the initial series can be reduced, as observed in equation ( 11).
The reduction of the series by the binomial development allows rewriting equation ( 10) into a more compact form, as shown in equation ( 12).
The total displacement e(t) throughout the column comes from the integration of equation ( 12), as indicated in equation ( 13) in which the first derivative of that function is represented by an upper row to the right.
The real and virtual displacements and their derivatives, with the same notation established in equation ( 13), which are expressed in function of the generalized coordinate and the chosen shape function to represent the considered vibration mode, are calculated by equation ( 14).
By equating the expressions ( 15) and ( 16) and adjusting the terms, the equation of the free undamped L v x t x q t v x t x q t v x t x q t v x t x q t v x t x q t v x t x q t v x t x q t v x t x q t e v x t v x dx

_
Creep in the fundamental frequency and stability... motion appears regarding the generalized coordinate q(t), as equation ( 17), in which M, K 0 and K g are the generalized mass and the generalized conventional and geometric stiffness, described from the chosen form function . By assuming the trigonometric function to describe the vibratory motion, the generalized conventional stiffness is given by equation ( 18), where E(t) is the variable modulus of elasticity with time, as shown in equation ( 7), and I is the inertia of the section in relation to the considered motion.The generalized geometric stiffness is calculated by equation ( 19), in which N(x) = [m 0 +m (L -x)]g is the normal force function, which includes the own weight of the column and the force concentrated at the free end.The generalized mass is obtained by equation ( 20), where m 0 is the mass concentrated at the top of the bar and m is the mass per unit length obtained by multiplying the section area by the density of the material.As the natural frequency depends directly on the total stiffness and inversely on the mass, it should be calculated by equation ( 21).
Therefore, the equation of the frequency of the first mode of vibration, in Hertz, which includes the geometric effect and creep, is finally provided by equation ( 22), in which L is the length of the column, g is the gravity acceleration and E(t) the temporal elastic modulus as found in equation ( 7).In equation ( 22), it can be observed the conventional portion of the stiffness of the column to the left of the negative sign in the numerator, in which is included the variable elasticity modulus due to creep, and the stiffness geometric portion to the right of the same sign, in which it is considered the force concentrated at the free end of the column by mass m 0 in addition to its own weight by the parameter m, which is the mass distributed per unit length, both multiplied by the gravity acceleration.In the denominator is the generalized mass of the system, where there is a factor of 0.227 multiplying the product mL.This fact induces which buckling critical load of the column is seen decreased from this value, when the own weight of the structural element is introduced, which is in line with Timoshenko (1992) who predicted, in static analysis, that this factor would be of the order 1/3.By applying the equation ( 22), it is performed, using a single mathematical operation and at once, a non-linear material analysis, with the creep, and a geometric non-linear analysis, with the effect of the normal force, dispensing the use of computational tools or iterative calculations.It is worth mentioning that the equation ( 22) was comparatively evaluated by experimental activity, in physical laboratory, and by modeling using the finite element method -FEM, by Wahrhaftig; Brasil; Balthazar (2013), and again using the FEM by Wahrhaftig and Brasil (2016), studying the frequency of a real structure, considering for both studies, linear and isotropic elastic material, presenting excellent convergence of results in both cases.

Evaluated model problem
For this simulation, a composite section wooden column was idealized, with the parameters of interest shown in Table 1.The adopted gravity acceleration g was 9.81 m/s 2 .The characteristic parallel resistance to fibers (f c0k ) and the elastic modulus (E co ) were adopted considering class wood C-60.The effective elastic modulus and calculating the resistance of the material were defined according to the recommendations of NBR 7190/1996-Wood Structural Design (1996), according to equation ( 23).The used weighting and modification coefficients were: K mod1 = 0.7, K mod2 = 1.0,K mod3 = 0.8 and  wc = 1.4,then: The numerical simulation was performed based on a column of height L and cross section as represented in Figure 3.The arrangement of the composite section was defined within a design reality typical of wooden structures and in function of the commercial limitations of laminated and glued wood, which had already been adopted in a design carried out by Wahrhaftig and Carvalho (2016).It is important to mention that the viscous parameter was adjusted so that the deformations stabilized at 90 days, as shown in Figure 4(a), as indicated by Almeida (1990), Celia-Silva;Calil Júnior (1992) and Kataoka;Bittencourt (2014).Thus, the variation of the elastic modulus E(t) was obtained, represented by Figure 4 (b).
It is emphasized the need to reduce the section inertia moment by multiplying it by a reduction factor of 0.7, since it is a composite section, constructed with screws, which is an important procedure to provide protection against plastic accommodation, specific of wooden connections, which was studied by Kharouf;Mcclure;Smith (2003) and Wilkinson;Rowland;Cooks (1981).

RESULTS
The column section was verified in relation to the flexocompression and stability according to the criteria of NBR 7190/1996-Wooden Structures Design (1996), proving to be able, as the results presented in equation ( 24), in which K M = 0.5, and equation (25), respectively.
In equations ( 24) and ( 25),  Nd ,  Md and f c0d are the calculation stresses due to the normal stress and bending moment, and calculation resistance to parallel compression to the fibers, obtained according to equations ( 26 In equation ( 26) subindices c for N, E and  mean compression; E in F is the buckling critical force of Euler; for e, the subindice i means initial, a accidental, c creep; for w in g means weighting.Regarding the stress conditions in the design g and q in N, they indicate variable and permanent actions; k means characteristic, d for the design and ef effective value, calculated with the modification coefficients, shown in section 2.2.

DISCUSSION
The frequency of the column was then numerically calculated at the instants zero and 90 by equation ( 22), as shown in Figure 4(c) and Figure 4(d).In the first case, the frequency variation was obtained by considering the original height of the column, adopted to obtain the limit slenderness established by the NBR 7190/1996-Wooden Structures Design (1996).In the second case, the height of the column was designed to initiate its collapse at 90 days, defining thus the limit height it can have when the creep is introduced in the frequency calculation.In the first condition, the column is stable, whereas in the second case, the creep effect leads to loss of stability at the considered instant (90 days).Any height of the structure which is defined between the stability limit at 90 days (second condition) and the limit established without creep, i.e., stability at instant zero (first condition), puts at risk the system equilibrium, which can be observed in Figure 4(e).For a height of 8 m, for example, the collapse would occur on the 16 th day.

CONCLUSION
The elasticity modulus of 6865.265MPa calculated by equation ( 7), at 90 days, represents a 50% decrease compared to the initial value of 13720 MPa.The structure frequency calculated at the initial time of loading is 0.262 Hz, and at 90 days is 0.106 Hz, representing a decrease of 60%.The column reaches its stability limit at the height of 7.35 m (slenderness of 156.19), collapsing at 90 days.Without considering the creep effect, the height limit is 10.37 m, at instant zero (slenderness of 220.41); a height of 29.14% higher than the previous one.This aspect is of great importance because, since if the height of the structure had been defined between the limit established without the creep and the limit determined considering the creep, the structure would collapse even before completing 90 days in service.For a height of 8 m (slenderness of 169.99), e.g., the collapse would occur close to the 16 th day.
It can be concluded that the slenderness limit of 140, established by the Brazilian Standard, keeps the structure in safety condition regarding the loss of stability, considering a force 10% of the Euler force concentrated at the free end of the column.Additionally, it can be added that the safety limit condition is established by equation ( 24) for a critical load of 44.41% of the Euler Force at the end of the structure when it is simultaneously evaluated by equation ( 22).
(a) Kelvin-Voigt viscoelastic model of three parameters -Group I (b) Kelvin-Voigt viscoelastic model of three parameters -Group II (c) Burges viscoelastic model of four parameters -Group III (d) Viscoelastic model of four parameters -Group IV e  Elastic

Figure 2 -
Figure 2 -Mathematical model.Figura 2 -Modelo matemático. 0 ), with the results shown as MPa (10 6 Pa; Pa = N/ m 2 ).The force concentrated at the free end of the bar, which induces the stresses  Ncd and  Md was established to represent one-tenth of the Euler force (F E ) expected in equation (26)(c) where = 0.8,  1 = 0.6 and  2 = 0.4 are the creep coefficients appropriate to the case under study.
Creep in the fundamental frequency and stability...