Modelling the hygro-mechanical creep behaviour of FRP reinforced timber elements

A fully coupled moisture-displacement finite element model has been developed to predict the viscoelastic, mechano-sorptive and swelling/shrinkage behaviour of FRP reinforced timber elements when stressed under long-term load and simultaneously subjected to changes in relative humidity. A DFLUX subroutine, to describe the changes in relative humidity with time, and a UMAT subroutine, implemented to describe the viscoelastic, mechano-sorptive and swelling/shrinkage behaviour, are presented. Additionally, an irrecoverable mechano-sorptive component is presented. This additional irrecoverable component occurs when the timber moisture content increases to a moisture content above levels previously attained. The model is found to be in agreement with experimentally determined deflection and strain results on both unreinforced and FRP reinforced beams subjected to constant and variable climates. 2020 Elsevier Ltd. All rights reserved.


Introduction
FRP (Fibre Reinforced Polymer) materials are increasingly being used to strengthen and stiffen structural timber products [1][2][3]. This occurs in both new and existing timber construction and has been used to great effect when retrofitting structures. Changes in use of the building or, indeed, changes in building regulations often require a higher load capacity than that of the existing members. FRP reinforcement has been used successfully to achieve such additional capacity requirements [1][2][3][4][5][6][7][8][9][10] but the influence of FRP reinforcement on the long-term or creep behaviour of timber elements has received less attention. Creep in wood is often separated into two main categories, namely viscoelastic creep and mechanosorptive creep. Viscoelastic creep is defined as the deformation with time at constant stress under constant environmental conditions and it has been the subject of many investigations, which have examined the influence of stress [11,12], temperature [13], moisture content [14] and flexural reinforcement on the longterm viscoelastic behaviour of timber elements [15][16][17][18]. Mechano-sorptive creep behaviour is defined as a deformation due to the interaction between stress and moisture content change in timber elements. Under variable environmental conditions, the mechano-sorptive effect can greatly accelerate the creep deflection of a timber element and ultimately lead to failure. The prediction of long-term effects and the influence of reinforcement on these timber beams are of crucial importance to structural engineers when designing timber structures.
Rheology is the study of the flow of matter and rheological models can be used to describe time-dependent viscoelastic behaviour in timber elements. The basic elements of these models are known as spring and dashpot elements. When combining springs and dashpots in series and/or parallel, rheological or viscoelastic models can be created [19]. In practical applications, generalised models are required to accurately model viscoelastic material behaviour. Generalised models combine a number of springs and dashpot elements, which increases the number of parameters and gives a better representation of the behaviour of the material being studied [11,20]. When modelling the long-term behaviour of timber elements, the surrounding environment must be considered given its significant effect on the material properties and creep behaviour. Many researchers have attempted to predict the effects of variable humidity and associated moisture fluctuations on timber elements and to model mechano-sorptive behaviour [21][22][23][24][25]. Leicester [21] and Ranta-Maunus [22,26] developed some of the first models in an attempt to quantify the mechano-sorptive effect in timber. Leicester [21] produced a uniaxial model consisting of an elastic component and an irrecoverable mechanosorptive component. The predicted results were compared with experimental data from a small sample of Eucalyptus during two weeks of drying. This rheological model was found to predict 85% of the total deformation recorded during testing. In a later study, Fridley et al. [27] developed a sophisticated model based on the Burger model, which included mechano-sorptive creep effects, but also included moisture and temperature dependent variations in stiffness. This adaption provided good agreement between modelled and experimental strain data. Mohager & Toratti [28] presented a mechano-sorptive model incorporating different tensile and compressive compliance functions and a model including partially irrecoverable compressive strains was proposed. Mårtensson [25] developed another model that decomposed the total strain into elastic, viscoelastic, shrinkage/swelling, and mechanosorptive creep strain components. Lu & Leicester [29] presented a mechano-sorptive model coupled with a simple approximate solution for moisture variation within a specimen based on Fick's law of diffusion and sinusoidal fluctuations in ambient moisture boundary conditions. Hanhijärvi & Mackenzie-Helnwein [30] presented an orthotropic material model based on a onedimensional model presented by Hanhijärvi [31] incorporating elastic, viscoelastic, mechano-sorptive and swelling/shrinkage strains. The mechano-sorptive strain component is altered to incorporate permanent deformations by introducing a hardening type plasticity element. This was found to accurately describe the deformations associated with high temperature drying in timber elements. Fortino et al. [32] presented a three-dimensional model for analysing timber structures under variable humidity and stress conditions. The viscoelastic/mechano-sorptive constitutive model is implemented in the user subroutine UMAT of the finite element modelling code Abaqus. They used Fick's law to describe moisture flow within the timber element. This Fickian behaviour was implemented into the Abaqus user-defined DFLUX subroutine. The model was successfully validated against experimental results by Leivo [33], Toratti & Svensson [34], Svensson & Toratti [35] and Jönsson [36]. Such advanced models have been used to examine a range of different challenges facing the timber industry and have contributed to increased knowledge and safer structures. Fragiacomo et al. [37] utilised an advanced numerical model to compute the moisture distribution and corresponding stresses or moistureinduced stresses in timber elements subjected to different climates which demonstrated that northern European climates can be more severe than southern climates. Massaro et al. [38,39] applied such models to examine timber elements subjected to compression perpendicular to the grain. This has been further advanced in recent years to included multi-Fickian transport theory which examines the influence of changes in both water vapour and bound water on timber structures [24,40]. Huč et al. [40] developed a model to successfully predict strains in timber elements subjected to sustained compressive loading. Fortino et al. [24] recently presented an advanced study incorporating a multi-Fickian hygro-thermal model capable of simulating severe changes in relative humidity and temperature. The study presented the associated moisture fluctuations and moisture-induced stresses on a bridge in the north of Sweden which is subjected to harsh climatic conditions over a 10-year analysis. While many models have been developed to examine the effect of FRP reinforcement on the mechanical properties of timber elements, the majority of these models focus on the short-term performance and are not concerned with viscoelastic creep, mechano-sorptive, temperature and moisture effects. A twodimensional model was developed by Tingley [41] to examine the stress-strain relationship of FRP reinforced glued laminated beams. This model did not include plasticity, creep, temperature or moisture effects and solely focused on the short-term elastic response. Serrano [42] modelled timber connections with bonded-in rods. The timber was modelled as a linear elastic orthotropic material and it was found that the stress concentrations were greatly influenced by the reinforcement stiffness and thickness. In another study, Alam [43] implemented anisotropic plasticity in a three-dimensional model to predict flexural properties of LVL beams reinforced with bonded-in FRP reinforcement. This model compared well with experimental findings. Raftery & Harte [44] also employed anisotropic plasticity theory in a model to predict the behaviour of unreinforced and glass fibre reinforced glued laminated beams. The model results agreed strongly with experimental results and a parametric study showed that as the reinforcement percentage is increased, the stiffness and ultimate moment capacity are enhanced.
A relatively low number of studies have attempted to predict the long-term performance of FRP reinforced beams under constant or variable climate conditions and these have been limited to the uniaxial case. In a study by Plevris & Triantafillou [15], an analytical model is employed to predict the creep behaviour of CFRP reinforced timber beams. The model, based on a Burger element, was developed to account for moisture effects. The experimental and analytical data showed good agreement in a constant climate; however, only limited experimental data was available. In Davids et al. [45], a uniaxial model was developed to model the creep behaviour of GFRP reinforced timber beams in a variable climate. The model, combining viscoelastic and mechano-sorptive creep strain parameters, was fitted using unreinforced specimen test results and then shown to accurately predict the relative creep deflection of glass FRP reinforced beams in a variable climate.
Numerical modelling can be an important tool when modelling complex materials such as timber. It has also been shown to adequately model various reinforcement materials in various configurations. These advances in numerical modelling have increased the reliability of results and as a result can aid future product development. The recent advances in modelling unreinforced timber by Hanhijärvi & Mackenzie-Helnwein [30] and Fortino et al. [24,32,46] are yet to be specialised for reinforced timber elements under changing relative humidity conditions.

Objectives of the current study
The objective of the current numerical study is to develop a model to predict the influence of flexural reinforcement on the long-term creep behaviour of timber beams in constant and vari-able climates. The numerical model presented is developed using Abaqus finite element analysis software [47]. As mechanosorptive creep is not explicitly available in Abaqus, a userdefined UMAT subroutine is required to define mechano-sorptive behaviour during loading and simultaneous moisture diffusion within timber. The moisture content distribution within the timber and associated strains are determined in a fully coupled timedependent hygro-mechanical analysis. A DFLUX subroutine is utilised to describe the interaction between the surrounding environment and the surface of the timber in the model. Such models have been shown to successfully replicate deflection and stress development in fluctuating environments by Mirianon et al. [48] and Fortino et al. [32]. In this study, the formulation of the UMAT subroutine follows the approach implemented by Santaoja et al. [49], Mårtensson [50], Mirianon et al. [48] and Fortino et al. [32]; however, a significant alteration was made to the irrecoverable mechano-sorptive creep matrix to include irrecoverable deformations in the longitudinal direction when stressed under loaded conditions and subject to moisture contents not previously attained. The coupled hygro-mechanical model is compared to the experimental creep deflection and strain tests results on unreinforced and FRP reinforced elements subjected to four-point bending tests presented by O'Ceallaigh et al. [17,[51][52][53][54].

Materials
In this study, timber is modelled as an orthotropic elastic material in Abaqus. The Sitka spruce material tested experimentally by O'Ceallaigh et al. [17,[51][52][53] was graded to strength class C16. The standard EN 338 [55] provides a mean value of 8000 MPa for the elastic modulus of C16 timber in the longitudinal direction; however, the experimentally determined elastic modulus in the longitudinal direction of 9222 N/mm 2 , determined from short-term static bending tests, is used. The elastic moduli in the radial and tangential directions were not measured experimentally; however, Bodig & Jayne [56] expressed a general relationship between the elastic modulus, E and shear modulus, G in the three material directions. These ratios are reproduced in Eqs.
These ratios are widely accepted for modelling the orthotropic properties of timber. However, it is noted that these ratios are not without their inaccuracies. The size of a timber element and the location it was cut from a log are factors in the accuracy of the assumption. More information on the timber material properties is given in Table 1. Additionally, these material properties are also dependent on the density, temperature and moisture content of the timber. Eqs. (4) and (5) are implemented to adjust the elastic and shear modulus properties when the environmental conditions are different from the reference conditions [49].
and u refer to the current density, temperature and moisture con-tent, respectively, at any given time step and a 1 = 0.0003, a 2 = À0.007 and a 3 = À2.6 are material constants provided by Santaoja et al. [49] and Mirianon et al. [48].
The FRP rod reinforcement and epoxy adhesive in reinforced members are modelled as linear elastic materials. Experimental tests have shown that the reinforcement is negligibly affected by swelling and shrinkage or creep and have not been considered in the study [51]. The material properties implemented in the model are presented in Table 2.

Finite element formulation
In this section, the constitutive models used to simulate the moisture diffusion and deformation of timber beams in both constant and variable climates are described. The formulation presented for deformation follows the approach implemented by Santaoja et al. [49], Mårtensson [50], Mirianon et al. [48] and Fortino et al. [32] but is adapted to represent Irish-grown Sitka spruce using experimental data. In the fully coupled 3-dimensional finite element model, a UMAT subroutine to determine the total strain experienced in timber elements when subjected to mechanical stress and simultaneous moisture content changes with time is implemented. This simulated total strain is subdivided into five separate strain components. These strain components are solved for in an incremental analysis with time. The total strain e Total , is

Elastic strain component
The elastic component (e e ) follows the generalised Hooke's law as seen in Eq. (7) where r is the stress and C e is the elastic compliance matrix of an orthotropic material. In 3-dimensional matrix form, this equation defines the orthogonal elastic stress-strain behaviour of timber (Eq. (8)).
The elastic compliance matrix is symmetric, and the Poisson's ratios are related according to Eq. (9).
The subscripts L, R and T refer, respectively, to the directions of assumed elastic symmetry, longitudinal, radial and tangential directions. E and G represent the elastic modulus and shear modulus, respectively. Bodig & Jayne [56] expressed a general relationship between the moduli E (elastic) and G (shear) in the three orientations as presented in Eqs. (1)-(3). These material properties of timber are dependent on the moisture content, temperature and density of the timber. Eqs. (4) and (5) are implemented to adjust the elastic and shear modulus properties when the environmental conditions are different from the reference conditions [49].

Viscoelastic strain component
The additional deformation with time in timber elements loaded under constant stress is known as viscoelastic creep. In this study, the viscoelastic strain component (e ve ) is an extension of the one-dimensional model presented by Toratti [12] consisting of a sum of Kelvin type elements as illustrated in Fig. 1.
The stress-strain relationship describing the time-dependent behaviour of a Kelvin element is described in Eq. (10) [19,30,32].
where e ve,i, = viscoelastic strain of the i th Kelvin element and Dt = the time increment (t n+1 À t n ). By integrating this expression for viscoelastic strain (Eq. (11)) over time for the i th Kelvin element, the following increment of elemental viscoelastic strain is obtained [48,57]. where: where J ve,i represents the dimensionless viscoelastic compliance ratio of the i th Kelvin element.

Total mechano-sorptive strain component
The total mechano-sorptive strain component (e ms,Total ) consists of two individual components as seen in Eq. (15). Both of these components occur under coupled stress and moisture content change associated with relative humidity variations in the surrounding environment. The increment of the irrecoverable component of mechanosorptive strain is calculated using Eq. (16) [32,35]. It is independent of stress rate and is a permanent deformation proportional to the current stress state and change in moisture content [58]. De ms;irr;nþ1 ¼ C ms;irr r n DU j j ð16Þ where De ms,irr = irrecoverable mechano-sorptive strain increment, C ms,irr = irrecoverable mechano-sorptive compliance matrix (Eq. (17)), r n = stress state at time t = t n and DU = increase in moisture content above levels attained during previous moisture cycles. The irrecoverable mechano-sorptive matrix previously presented by Fortino et al. [32] accounts for irrecoverable deformations in the radial and tangential directions only. As a novelty of this research project, it was shown experimentally that the development of irrecoverable mechano-sorptive creep strains are also observed in the longitudinal direction. As a result, Fortino's formulation has been modified to account for irrecoverable strains (m L ) in the longitudinal direction as seen in Eq. (17). The value of the irrecoverable coefficient in the radial/tangential direction, m v presented by Fortino et al. [32] and the longitudinal direction, m L have been determined through careful calibration against experimental results on unreinforced beams.
where m L = Longitudinal irrecoverable strain coefficient and m v = Radial/tangential Irrecoverable strain coefficient. When describing the recoverable mechano-sorptive strain associated with repeating cycles of relative humidity, the schemes previously implemented by Santaoja et al. [49] Mårtensson [50] and Ormarsson [59] have been implemented. The solution is found using the absolute value of moisture content, which is assumed to be constant over the time increment. De ms;nþ1 ¼ C ms r n Du j j ð18Þ where: De ms = mechano-sorptive strain increment, C ms = mechanosorptive compliance matrix (Eq. (19)) r n = stress state at time t = t n and |Du| = the moisture increment |u n+1 -u n |. The elemental mechano-sorptive compliance matrix assumes the form given in Eq. (19) [49,50,59].
The coefficients of the matrix (C ms ) are based on the mechanosorptive compliance coefficients (m ms,i ) in the respective material directions and, using the ratios v ms,i , to describe the coupling of the mechano-sorptive strain between the different directions. The values used have been based on previous research by Santaoja et al. [49] and Ormarsson et al. [60] and finally matched to experimental test results.

Swelling/shrinkage strain component
Swelling/shrinkage strains (e s ) occur within the timber during moisture content change. This change in strain depends on the material direction and magnitude of the change in moisture content. The increment in e s due to a moisture content change, Du, is described by Eq. (20).
where the moisture swelling/shrinkage vector a u is described by Eq.
The components of the swelling/shrinkage vector (a u ) used in the study were determined experimentally [61]. It has been found that different loading conditions (tension, compression or bending) result in different longitudinal swelling/shrinkage behaviour of timber. Wood tends to swell/shrink less when in tension and swell/shrink more in compression [62]. Toratti [23] implements an adapted formula (Eq. (22)) to account for variations in longitudinal strain when in a state of tension and compression as is the case in bending.

Diffusion modelling
Fick's law, which describes the moisture diffusion process, is analogous to Fourier's law of heat transfer. Using this analogy, the governing equations allow diffusion coefficients and moisture content to replace thermal conductivity and temperature. This allows coupled temperature-displacement elements in Abaqus to be exploited in a hygro-mechanical model to simulate coupled moisture diffusion and mechanical behaviour. The environmental load or relative humidity load is implemented in Abaqus through a DFLUX user subroutine in order to apply a defined timedependent relative humidity boundary condition. The transport of moisture across the surface of the timber is proportional to the difference between the equilibrium moisture content of the surface of the timber and the equilibrium moisture content of the timber corresponding to the relative humidity of the surrounding environment. In the DFLUX subroutine, the rate of moisture transport between the surface of the timber and the surrounding environment is governed by Eq. (23).
where q n = rate of moisture flow, q = density at 0% moisture content, u eq = equilibrium moisture content of timber corresponding to the relative humidity of the surrounding environment and u surf =-moisture content of the timber surface. The surface emissivity S u defines the rate of moisture content exchange across the boundary and is a function of the environment, the characteristics of the timber surface and is also dependent on the velocity of the circulating air [58]. The value used in this model of 3.2 Â 10 À8 e 4.0u m/s was previously reported by Hanhijärvi [58] where u is moisture content. More information related to the diffusion coefficients and the influence of relative humidity on the equilibrium moisture content (u eq ) of the timber used in this study is provided by O'Ceallaigh et al. [61], which focuses on moisture-induced strains in fast-grown timber. In the current study, the environmental load is applied to the exposed surfaces of the beam. In line with the experimental tests [17,[51][52][53][54]], the exposed surfaces include both sides of the glued laminated beam and the top surface only. The BFRP reinforcement in the tension zone of the reinforced beams impedes the movement of moisture through the bottom and end grain surfaces. As a result, the bottom and end grain surfaces were sealed for both unreinforced and reinforced beams to ensure all beams are subject to the same exposure conditions.

Model geometry
The geometry of unreinforced and FRP reinforced beams, which are subjected to a creep test in four-point bending, [17,51,53] is modelled using Abaqus finite element software. As seen in Fig. 2, the 98 Â 125 mm 2 unreinforced and reinforced beams comprise four laminations. Each of the four laminations is assigned material properties in a local cylindrically orientated coordinate system. The glue-line between laminations was not modelled in this study due to the good performance of the glue-line under variable climates and its limited influence on the moisture diffusion behaviour through the cross-section of the timber as shown by Raftery et al. [64] and O'Ceallaigh et al. [61], respectively. The FRP reinforced beam model is similar to that of the unreinforced beam; however, in the bottom tensile lamination, two routed grooves are created to accommodate the two 12 mm diameter basalt fibre reinforced polymer (BFRP) rods and the 2 mm structural epoxy adhesive. Both models use 8-noded coupled thermaldisplacement C3D8T elements. As seen in Fig. 3, symmetry is utilised allowing only half of the 2300 mm long beam to be modelled to reduce the computational time. The respective mesh sizes for the unreinforced and reinforced beam models were determined from mesh sensitivity studies to provide accurate results in a reasonable time frame. The beam is simply supported on 80 mm Â 100 mm plates which are free to rotate about their central axis. These plates are modelled using 2-dimensional shell S3 elements. Hard contact is defined between the surface of the beam and the steel plate with a tangential friction coefficient of 0.4. Similar plates are used to apply the dead load in tests. These are also modelled using S3 shell elements and the load is applied as a uniform pressure over the plate area.
Regarding the numerical models, a load of 5749 N was applied to the unreinforced beams and a load of 6241 N was applied to the reinforced beams to induce a common maximum compressive bending stress of 8 MPa in each beam as was done experimentally by O'Ceallaigh et al. [17,51].

Material data
This section presents the material data for the timber, epoxy adhesive and BFRP implemented in the numerical models. The timber material data is summarised in Table 1. The BFRP rod reinforce-ment and epoxy adhesive material data implemented into the mechanical model are presented in Table 2.

Climate conditions
The numerical model is compared to experimental creep results on beams in both a constant and variable climate [17,[51][52][53]. The constant climate remains at a relative humidity of 65% and a temperature of 20°C for the duration of the test. The variable climate condition induces mechano-sorptive creep and swelling/shrinkage strains due to hygro-expansion. This occurs due to the change in the relative humidity and subsequent change in the moisture content of the beams. Testing commenced at a relative humidity of 65% and at a temperature of 20°C for a period of three weeks, after which, the relative humidity in the variable climate chamber was changed to 90%. After three weeks, the variable relative humidity cycle commenced with an increase in the relative humidity to 90% for four weeks followed by four weeks at a relative humidity of 65%.

Unreinforced beam model results
In this section, the results of the long-term numerical simulation of the unreinforced beam model in constant and variable cli-mates are compared to the mean experimental result presented by O'Ceallaigh [51]. The experimental groups presented comprised of nine beams each [51]. The predicted total vertical deflection is compared to the experimental results in Fig. 5. The predicted initial elastic deflection of 6.09 mm compares well with the mean measured values of 6.27 mm (Group UC) and 5.95 mm (Group UV). Over the 75-week test period, it can be seen that the total mean long-term vertical deflection of Group UC beams is simulated accurately with a final predicted deflection of 8.08 mm compared to a measured mean value of 8.03 mm. This is to be expected as the viscoelastic compliance coefficients defined in Table 1 and used in this unreinforced model have been calibrated against the mean experimental deflection data from the unreinforced beams. It was important to accurately define the viscoelastic compliance coefficients of the unreinforced beams prior to simulating creep in reinforced beams.
In the variable climate condition, during the first three weeks, the relative humidity remains constant and, as a result, this period is solely associated with elastic and viscoelastic creep. Consequently, the numerical result for Group UV produced the same elastic deformation and viscoelastic creep as the constant climate model during this three-week period. After the first three weeks, the relative humidity increased resulting in an increase in moisture content which in turn leads to a significant increase in deflection due to irrecoverable mechano-sorptive creep (Fig. 4.) Significant increases in deflection and strain in the first relative humidity cycle have been previously reported due to changes in moisture content higher than that previously attained under a loaded condition [68]. After week 7 the relative humidity reverts to 65% from 90% and the relative humidity continues this 8 week cycle for 75 weeks.   result of À864 le compared well to the numerical result of À929 le. The total mean strains measured on the tension and compression faces of Group UV in a variable climate have also been measured and are compared to the numerical model results in Fig. 5. The mean initial elastic strain has been slightly over estimated on both the tension and compression faces of the Group UV beams. After the initial elastic component, the model appears to adequately predict the total strain on the tension and compression face of the Group UV beams in a variable climate. The large increase in longitudinal strain during the first moisture content change attributed to the irrecoverable mechano-sorptive component has been simulated.
The fluctuation in strain with repeated cycles has been successfully simulated with the numerical model. On the tension face, the change in strain with each relative humidity cycle is more gradual than that experienced on the compression face. This is a result of the different exposure conditions of the tension and compression faces. The bottom (tension) face is sealed from the surrounding environment resulting in a reduced rate of moisture content change and more steadily increasing and decreasing rate of strain development in the tension face. After 75 weeks, the experimental tensile strain of 1294 le compared reasonably well to the numerical result of 1535 le and the experimental compressive strain of À1069 le compared well to the numerical result of À1174 le. It  can be seen that the results of the numerical model slightly over predict the mean experimental strain results measured on the tension and compression faces of the unreinforced beams. This over prediction is partially due to initial differences in the elastic strain.

Reinforced beam model results
The simulated numerical deflection for the reinforced beam model in constant and variable climates is compared to the mean experimental result of Group RC and Group RV, respectively in Fig. 6. Group RC and Group RV comprised nine beams each [51]. The addition of the BFRP reinforcement has resulted in a reduced initial elastic deflection and reduced long-term deflection after 75 weeks when compared to the unreinforced Groups UC and UV. The reduction in initial elastic deflection has been accurately predicted using the model. In the constant climate condition, the mean experimental elastic deflection of Group RC was 5.72 mm and the numerical model predicted a value of 5.76 mm. Thereafter, the numerical model closely matches the experimental creep deflection up to 75 weeks.
For Group RV, the measured elastic deflection of 5.44 mm compared well with the numerical prediction of 5.76 mm. Once the relative humidity was increased, the deflection increased rapidly. Experimentally, the initial increase in creep deflection was less in the reinforced than the unreinforced beams and this is also predicted numerically. With additional relative humidity cycles, the fluctuations in deflection have been successfully predicted; however, greater fluctuations in the predicted numerical deflections have been observed. After 75 weeks, the mean experimental Group RV deflection result of 10.10 mm has been accurately predicted by the numerical model with a simulated deflection result of 10.09 mm. It is noted that the peak deflection of each cycle matches the numerical model well, but the troughs are under predicted but the general trend of increasing deflection with repeated cycles has been accurately modelled.
The total measured and predicted mean longitudinal strains on the tension and compression faces of the reinforced beams are shown in Fig. 7. In the constant climate, the maximum tensile and compressive strains in the reinforced beams are significantly lower than in the unreinforced beams, demonstrating the beneficial effect of the reinforcement. For the reinforced Group RV, the elastic tensile strain of 597 le predicted by the model is slightly higher than the mean experimental tensile strain of 577 le. On the compression face, the predicted elastic strain of À718 le is larger than the experimental result of À643 le. Although the elastic strain on the compression face is larger than that measured experimentally, the difference is not deemed significant.
A large increase in longitudinal strain during the first moisture content change can be seen. This irrecoverable mechano-sorptive component has been simulated accurately in the reinforced beam. The reduction in irrecoverable mechano-sorptive creep due to the reinforcement, observed experimentally, is also predicted by the numerical model. The fluctuations in strain with repeated cycles simulated using the numerical model match well on the tension face but are overpredicted on the compression face. After 75 weeks, the measured maximum compressive strain was À866 le while the model predicts a conservative value of À1068 le. The reinforced beam model has adequately simulated the mean experimentally measured deflection and strain behaviour in a variable climate condition. The addition of the BFRP rod reinforcement into the model has been shown to simulate the beneficial reduction of vertical deflection in reinforced beams and has been shown to adequately predict the total longitudinal creep strain which comprises the viscoelastic, mechano-sorptive and the swelling/ shrinkage strain, on both the tension and compression faces of the reinforced beams.

Summary and conclusion
Unreinforced and reinforced glued laminated beams have been subjected to long-term creep tests in a variable climate condition [51]. The experimental creep tests confirm that reinforcing timber with a material of superior properties has a positive effect on the mean deflection and longitudinal strain results in a variable climate. This reduction is currently not accounted for in the design of timber structures. A coupled hygro-mechanical finite element model has been developed to predict the behaviour of FRP reinforced timber elements when stressed under long-term load and simultaneously subjected to changes in relative humidity. The model incorporates elastic, viscoelastic, mechano-sorptive and swelling/shrinkage behaviour described using a UMAT user- subroutine in Abaqus finite elements software. Included in the formulation is an irrecoverable mechano-sorptive component. This irrecoverable component characterises the creep behaviour of stressed timber elements when the moisture content increases to a moisture content not previously attained. The climate conditions are implemented using a DFLUX subroutine and can be easily adapted to mimic other environmental conditions.
The following conclusions can be formulated based on the investigations presented on the long-term behaviour of FRP reinforced timber elements loaded in a constant or variable climate.
The hygro-mechanical creep model has been shown to accurately predict the creep deflection of unreinforced and reinforced beams in both constant and variable climates. The model has also accurately predicted the reduced creep deflection behaviour of reinforced members in both constant and variable climates due to the addition of FRP reinforcement in the tension zone of the beams. The hygro-mechanical model has also accurately predicted the longitudinal creep strain measured on the tension and compression face of unreinforced and reinforced members. A beneficial reduction in longitudinal creep strain on the tension face of reinforced members, which was observed experimentally in both a constant and variable climate, has been simulated by the model. The irrecoverable mechano-sorptive creep matrix proposed by Fortino et al. [32] has been modified here to include the effects of irrecoverable creep in the longitudinal direction. This greatly improves the accuracy of the predictions of the long-term deflection and strain in both the unreinforced and reinforced beams in variable climates. It is recommended to include this modification in future studies, but it is noted that the parameters provided in this study have been calibrated against experimental test results on fast-grown Sitka spruce and should be checked when implemented on different species of timber.
Long-term testing can be expensive and time-consuming. This validated hygro-mechanical model can now be used to examine other engineered wood products in different climates more quickly and at a significantly lower cost. In relation to reinforced elements, this model will be used in future parametric studies to investigate the influence of FRP type, percentage reinforcement area, and stress level on the creep behaviour of FRP reinforced beams with a view to determining creep modification factors for design.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
neering, NUI Galway, in particular, Peter Fahy, Colm Walsh and Gerard Hynes, is acknowledged.