A one-dimensional elasto-viscoplastic model coupled to damage for the description of creep in wooden materials

This work focuses on the development of a model for the description of the tertiary creep phenomenon in wooden materials. We stared from an extended standard solid body model capable of best describing primary and secondary creeps. We then modify this model by introducing a damage variable to explain and model the rapid growth of viscoplastic strain during tertiary creep. We obtain a model comprising a reduced number of parameters (05) all physically interpretable; which can be easily determined from the results of creep tests. The proposed model has been tested using the experimental results of creep-rupture tests and it has been shown to be very suitable for describing the three phases of creep, with a relative error of less than 1%. The breaking time proposed by the model is lower, but very close to the experimental breaking time (Err = 0.01). The time to failure is easily accessible, thanks to the simplicity of our model, without necessarily going through heavy algorithms. This represents a significant advantage of our model, which in sum offers both a more realistic way of describing the three phases of creep by fully accounting for the phenomenon of damage during the tertiary phase, and a simple and fast way to analyze the rupture time, compared to other models in the literature. Our model is therefore presented as a good alternative for modeling the behavior of wood material under creep stress.


Introduction
The context of sustainable development and respect for the environment have in recent years aroused an evergrowing interest in the use of biological materials, especially those based on wood. These biological materials, due to their cellular composition (lignin, cellulose, hemicelluloses), can experience a state of deformation evolving over time even under stresses kept constant [1,2]. This behavior is known as creep and is defined as the growth of strain over time under constant loadings.
It appears to be of strategic importance for engineers, with the aim of offering better safety conditions as well as better predictions of the lifespan of the materials structures, to be able to describe and characterize creep deformations. Much work has been undertaken in this direction [1][2][3][4][5][6][7][8][9] to name only those. Tests show that wood up to a certain load limit exhibits linear viscoelastic behavior. It is possible to determine a stress level below which the strain under creep evolves towards a finite value [1,3,5]. In this case, the creep curve has only two parts namely the primary creep and secondary creep. The work of Foudjet and Guitard, [1,3] show that the mechanical behavior of the wood material is appreciably linear for loading rates lower than a third of the breaking load. For high stress levels, the material experiences irreversible deformations which remain even after stress cancellation [3,5]. Wood can therefore under certain stresses show plastic behavior. Also, under high stress and for fairly long periods of time, wood experiences the onset of the damage phenomenon which is responsible for the rapid growth of deformation during tertiary creep. The various experimental observations made on the mechanical behavior of wood material therefore show that it has, depending on the type and level of stress, an elastic, viscous and plastic behavior.
Wood is therefore a complex material, a complexity which is increased by the sensitivity of its mechanical properties to ambient climatic conditions [1]. The complexity of the wood material makes it difficult to model its behavior; so that it is almost futile to imagine a model describing the behavior of wood without resorting to some simplifying assumptions.
However, for the most faithful modeling, the simplifying hypotheses must not only be reduced in number, but also be based on experimental observations. Many models have been proposed for modeling the creep behavior of wood. They can easily be grouped into two classes, namely: mathematical models and rheological models [8]. Mathematical models use mathematical functions whose parameters are determined from laboratory test results for the description of the creep curve. They have the advantage of the reduced number of parameters required for modeling and suffer from the disadvantage of the lack of physical significance of these parameters. The mathematical functions generally used are the powers functions, the logarithmic and exponential functions [6,7,10]. The rheological models use combinations of springs, viscous pots and pad for the modeling of elasticity, viscosity and plasticity, respectively. They require a relatively large number of parameters, which is their major disadvantage, but on the other hand, their advantage is that they have parameters that are physically interpretable. The classic models commonly encountered are those of Maxwell (association of series of spring and viscous pot) Kelvin (parallel association), Burgers [11][12][13]. Some of these rheological models have shown their efficiency in capturing primary creep (Kelvin) and the onset of secondary creep (Burgers) [13,14], but have difficulty in describing the end of secondary creep and tertiary creep. Modifications are usually made to these models to describe the end of secondary creep and material failure. In the wake of these approaches, we find the work of authors such as [13,15] to quote a few.
One of the crucial points in modeling wood creep is the capture of tertiary creep and failure. Tertiary creep remains insufficiently described. Indeed, the models proposed for the modeling of tertiary creep and the fracture of wood are generally focused on the determination of the final rupture time [16][17][18][19][20]; and ignore the process of deformation evolution during tertiary creep phase. However, understanding and describing this deformation process is just as important for a full understanding of the phenomenon of creep. This problem has recently been the subject of the work of authors such as Wang et al. [16], Kaliske [21]. Wang et al. modified the Burgers model by introducing a damage variable to describe the creep failure of wood materials, their model was shown to be effective in predicting the failure time, but was not able to capture the rapid growth of the strain which precedes the failure. Moreover, the model proposed by Wang et al. does not take into account the existence of a threshold plasticity stress such as the work of Foudjet, Guitard, and others [1][2][3]5] has shown. The model proposed by Kaliske et al. [21] takes into account the existence of a linear viscoelastic limit, and describes tertiary creep and in particular creep failure. The approach used by Kaliske et al. to define the failure criterion is an approach based on strain energy density, introduced and developed by authors such as [20,[22][23][24]. This approach, although it has shown its effectiveness in defining the failure criteria for certain types of stress such as compression and tension [24], nevertheless remains pregnant for a large number of stresses. To our knowledge, there is no work allowing an unambiguous description of an energy failure criterion for all types of stresses. In addition, the description of the evolution of the deformation of the tertiary creep as proposed by Kaliske et al. does not take into account the phenomenon of damage, damage which nevertheless has a fully established responsibility in the occurrence and development of tertiary creep [2,25,26].
The objective of the present work is to propose a model which describes the elastic, viscous and plastic behavior of wood materials during creep, while taking into account the phenomenon of damage for a better description of tertiary creep and fracture by creep. The proposed model is a rheological and one-dimensional model. The model will be validated by the results of tests carried out by Wang et al. [16]. Our work includes the following sections: "Development of the model", where we present the hypothesis and equations which have led to the development of the model as well as the method of determining the parameters, Validation of the model, where the results are presented and commented on and "Conclusion" section.

Development of the model
The model is developed with the aim of describing as faithfully as possible the experimental observations made on the behavior of the wood material under creep stress. These observations are [1][2][3][24][25][26]: -Instant elasticity, -Viscoelasticity, -Damage and rupture.
As the wood material is a biological material, its mechanical behavior is strongly influenced by its humidity level. This humidity level is also responsible for another phenomenon observed experimentally, namely mechano-sorption. This phenomenon will not be taken into account in the present model. In addition, the humidity level will be assumed to be constant for our material.

Instantaneous elasticity
This is the instantaneous deformation of the material upon application of the stress. This deformation disappears entirely as soon as the stress is canceled. It is modeled by a purely elastic spring of stiffness E. The stress-strain relation is given by Hooke's law:

Viscoelasticity
It is the growth of strain under constant stress. This deformation does not cancel instantly upon removal of the stress; it gradually disappears during the so-called recovery phase. The evolution of viscoelastic deformation is a function of the level of stress. [1,3] In the literature, it is frequently modeled by simple or generalized fields of Kelvin-Voigt or Maxwell [14,24]. One of our objectives being to take into account the existence of a finite limit in creep for low stress levels, and this with a reduced number of parameters, we carry our choice on the model of Kelvin-Voigt; because, Maxwell's body reflects the behavior of a fluid without the possibility of equilibrium for constant and non-zero stress [27]. The stress and strain relation for this model is given by: where the elastic modulus of the spring E 1 is, η 1 is the viscosity coefficient of the viscous pot, and ε ve (t) is the viscoelastic strain.

Viscoplasticity
The works of Van der put, Gressel and Hasanni [5,13,14] show that wood, for certain stress thresholds has a viscoplastic behavior where the stress is a function of the strain rate, which strain does not cancel out when the stress disappears even after a long enough time, thus leaving a residual strain or even plastic strain. Plastic behavior is frequently modeled using a pad, and viscoplasticity by a series and/or parallel association of a pad and a damper.
One of our objectives being to take into account the existence of a threshold conditioning the birth of nonlinearity phenomena, we adopt a simple Bingham element, which is a parallel association of a damper and a pad. The stress strain relation for this model is: where η vp is the viscosity coefficient of the damper, σ y is the threshold plastic stress, and ε vp is the viscoplastic strain.
We have so far individually modeled each of the phenomena (elasticity, viscoelasticity, viscoplasticity) observed experimentally by single and/or combined elements. It is now a question of combining all these individual elements for the modeling of the general behavior of the material. The combination must, however, have a reduced number of parameters and be thermodynamically valid (possibility of partitioning the total deformation among others) in addition to having physically interpretable parameters and to describe fairly faithfully the behavior of the material. A series association of the spring, the Kelvin-Voigt body and the Bingham body leads to an extended standard solid body model ( Fig. 1), [21].
This model admits the possibility of the partition of the deformation and presents a reduced number of parameters which are all physically interpretable. Note that if we take σ y = 0 , we obtain the Burgers model, which is known [8,13] for its ability to best describe primary creep and the onset of secondary creep. We will therefore adopt the model for the rest of our modeling.
The partition of the total deformation makes it possible to write: And the serial association makes it possible to write:

Damage and failure
For fairly long durations, and stresses greater than a certain threshold, the creep deformations of the material experience an exponential growth which leads to the rupture of the material. This is the last phase of creep or tertiary creep. The literature offers several

Fig. 1
Extended standard solid body model [21] approaches for the description of this phenomenon. Authors such as Wang, Cai, Pierce et al. [13,15,16] have modified the components of the model by introducing a damage variable to explain tertiary creep.
This approach, although yielding significant results, encounters a certain number of difficulties, namely firstly the too large number of parameters and the difficulty in finding a physical interpretation of each parameter. The other difficulty remains the basic rheological model chosen for this modeling, namely the Burgers model. If the Burgers model is known for its ability to properly describe primary creep and the onset of secondary creep, it remains insufficient to describe the existence of a finite limit of creep strains under low stress levels. However, any model which must serve as a basis for modeling damage and therefore tertiary creep must first be able to model as faithfully as possible the entire behavior of the material in the absence of damage. The model proposed by Reichel and Kaliske and adopted as a basic model within the framework of the present work captures very well the primary and secondary creep, and materializes well the finite limit for low levels of stresses.
However, for the description of tertiary creep, they used an approach based on strain energy density. This approach measures the strain energy e(t) defined by: This strain energy measured at each increment of the strain is compared with a threshold value characteristic of each material, threshold value beyond which the material breaks. This approach has been shown to be effective for the description of creep failure in that it allows a failure criterion to be defined. However, according to Becker [18], it remains insufficient to describe the rapid growth of strain during tertiary creep and is therefore unable to capture tertiary creep. To solve this problem, Reichel and Kaliske modified the viscoplastic viscosity coefficient by adding to it a factor governed by an arc-tangent function which decreases this coefficient when the strain increases. This approach, although having led to the capture of tertiary creep, does not take into account the phenomenon of damage, damage which however has a fully established responsibility in the occurrence of tertiary creep [2,25,26]. Indeed, the rapid growth of deformation during tertiary creep is due to the progressive decrease in rigidity, itself caused by the multiplication of gaps within the material, and therefore by damage. It therefore appears more reasonable and more realistic to take into account a damage variable in the modeling of tertiary creep. Moreover, this approach used by the authors introduces two parameters in the tangent arc function, parameters which are not easily accessible.
The originality of our approach in our modeling consists in modifying the component of the rheological model responsible for Viscoplasticity, by introducing a damage variable to better describe the rapid growth of the strain during tertiary creep, and thus model the tertiary creep. Our approach will make it possible to take into account the damage phenomenon so observed in the tertiary phase of creep, and thus lead to a model that is more faithful to reality and at the same time simpler. This modification will be made based on a set of assumptions all based on experimental observations.

Hypothesis 1
The damage begins after a certain threshold plastic deformation. It has been observed that the initiation of the damage is consecutive to the large plastic deformations [25,26]. There is therefore a threshold of plastic deformation conditioning the onset of the damage phenomenon. The damage variable therefore acts on the component of the solid responsible for viscoplasticity. We therefore assume that there is a defined threshold value ε D such as:

Hypothesis 2
The damage affects the coefficient of viscosity. The rapid growth of strain during tertiary creep is due to the gradual decrease in the coefficient of viscosity; the decrease is caused by the loss of rigidity, itself caused by the multiplication of gaps in the material, and therefore by damage. We therefore consider that the coefficient of viscosity of our material experiencing damage D is given by:

Hypothesis 3
The damage, once initiated knows an exponential growth, according to the stress level. We therefore admit that the damage is given as a function of time and of the stress by the relation: where t 0 is the instant at which the damage begins, and β a parameter conditioning the effect of the stress on the damage. The form of Eq. (9) is chosen to ensure that the damage evolves in the interval [0, 1], [27]. The instant t 0 is defined as: By putting Eqs. (8) and (9) in Eq. (3), we obtain the following modified rheological model (Fig. 2), with D given by Eq. (9).

Expression of the total strain
The total strain ε(t) is given by relation (4) respecting the hypothesis of the partition of the total strain, in elastic, viscoelastic, and viscoplastic strain.
The elastic deformation is obtained from Eq. (1) taking into account relation (5) and is given by: The viscoelastic strain is obtained by solving Eq. (2) and is given by: The development of the viscoplastic strain being dependent on the stress level, to express them, we consider the situation where the stress is greater than the plasticity threshold. So σ ≥ σ y ; thus, we can distinguish two situations. The first situation is where the plastic deformations evolve, but have not yet reached the threshold causing the damage to start. In this case, and in accordance with relation (7) obtained from hypothesis 1, Eq. (3) allows us to write: An integration of (14) makes it possible to write: In Eq. (14), the integration constant c is determined by considering the condition: This gives: The viscoplastic deformation before the onset of damage therefore has the following expression: The second situation is when the plastic deformations have evolved to the point of causing the onset of damage. Thus, Eqs. (3), (8) and (9)  A rewrite of (20) gives us: An integration of (21) gives us: The integration constant ε vp0 is determined by looking at the instant (t 0 ) when the damage begins. We have from Eq. (17): To ensure the continuity of the deformation to the left and to the right of the instant when the damage is initiated, we have, starting from Eqs. (22) and (23): From relation (24), we obtain the integration constant ε vp0 , which gives: By inserting the constancy ε vp0 obtained in (24) into Eq. (22), we finally have the expression of the (16)  Thus, the total deformation is given from Eq. (4) by: ε = ε e + ε ve + ε vp . This gives, by grouping the results (11), (12), (17) and (26): The system of Eq. (27) describes the dynamics of the rheological model that we propose for modeling the damage of wood material under creep stress.
This model that we propose requires a total of five independent parameters all physically interpretable and easily accessible ( E e ,E 1 ,η 1 ,η vp and β ), which is low, compared to the model proposed by Kaliske et al. [21] which has a total of six independent parameters. The model captures the rapid change in strain during tertiary creep while accounting for damage.

Determination of the model parameters
The parameters,E e , E 1 , η 1 and η vp can be determined through low stress creep tests. The parameter σ y representing the plasticity threshold is determined experimentally through creep-recovery tests under high stresses. Its empirical value proposed in the literature is 35 to 50% of the breaking stress of the material [1,3,5]. The parameter β conditioning the evolution of the damage and the plastic strain's rate can be determined from the strain rate of tertiary creep. It is given by the slope of the logarithm of the tertiary creep strain rate. Indeed, starting from Eq. (21) we have: Relation (28) gives the logarithm of the tertiary creep strain rate. As shown by this relation (28), the logarithm of the tertiary creep rate is a line whose slope is the product of the stress by the parameter β . From this slope, we therefore deduce the parameter β by: where tan α is the slope of the line representing the logarithm of the tertiary creep rate. The parameter t 0 depends on the stress σ and the plastic deformation ε D damage threshold and can be obtained through relations (10) and (27b), which allow us to write: The numerical resolution of Eq. (30) provides the parameter t 0 .
One of the essential points for determining the parameters of the model is the determination of the instant t 0 from which the damage begins. In general, in the models proposed for the description of the failure by creep of materials, the moment t 0 from which the damage begins is very often ignored, in favor of the moment t r of failure. Indeed, authors are generally more interested in the moment t r of failure than the moment of onset of damage. It seems important to us to be equally interested in the moment at which the damage begins in the modeling of failure by creep, because for certain materials, like brittle materials, the failure is so fast or sudden, that one can easily think that the moment of the beginning of the damage is almost the same one as that of rupture. This moment being linked to the threshold deformation ε D of damage by relation (30), it is necessary and sufficient to determine the threshold deformation of damage ε D .
When the damaging stresses ( σ ≥ σ y ) are kept constant, the creep strain can only increase. This growth necessarily begins with the value: ε 0 = σ E e which is the instantaneous elastic strain. The strain field of creep failure is therefore: The onset of damage being consecutive to the large plastic deformations, it is evident that the threshold damage deformation is included in the universe of deformations defined by relation (31). It is therefore logical to think that any increment �ε of strain starting from the value ε 0 = σ E e results in a total strain value which may be a good candidate to represent the threshold strain for damage. One of the possible increments of strain is: The choice of this increment is motivated by the need to propose a model whose parameters depend on the easily accessible mechanical characteristics of the materials (modulus of elasticity, breaking load). This increment leads to a total strain value given by: Relation (34) allows us to write: The total strain given by relation (34) is therefore a good candidate to represent the value of the threshold damage strain. For the validation of our model, we adopt the value of the deformation of the relation (34) as the value of the threshold deformation of damage. Thus, we have: The choice of this value is justified by the need to take into account the fact that the origin and development of damage strongly depend on the stress level and the (33) ε = ε 0 + �ε.
mechanical characteristics of the material (elastic modulus, breaking load).

Validation of the model
For the validation of the model, we use the experimental data obtained by Wang et al. [28]. The authors performed creep tests for up to a year on wood-based composites. Two batches of test specimens were requested under two different stress levels, the first batch loaded at 27 MPa and the second loaded at 33 MPa. The mechanical characteristics of the material subject of this study were determined and are presented in the following table (Table 1).
The parameter E e , modulus of elasticity is determined from the slope of the stress-strain relationship at the beginning of the loading. The parameters E 1 , η 1 and η vp , are determined from the least squares method identical to that used by Wang et al. [16], using the SOLVER procedure in Excel software [29].
The determination of the parameters of the model led to the results set out in Tables 2 and 3. The calculation of the relative error of the model for each set of parameters was made according to the relation: where N is the number of points, ε cal (t i ) and ε exp (t i ) are, respectively, the calculated value and the measured value of the strain.
The analysis of the mean and standard deviation for each parameter shows a low dispersion for some ( E e and (36)  E 1 ) and a high dispersion for others ( η 1 ,η vp and β ). This high dispersion is mainly due to the heterogeneity of the mechanical properties of the material. In addition, some samples may have initial defects in their breasts that are created either during sample preparation or during packaging. These initial defects accelerate the process of damage to the material during stress, cause very rapid rupture as observed by the authors for certain specimens, and thus justify, to a certain extent, the high dispersion of certain parameters of the model.  The following figures (Figs. 3, 4 and 5) show the comparison between the experimental creep curves and the curves calculated from the model. In Figs. 3, 4 and 5, one distinguishes the three phases of creep (primary, secondary and tertiary, as well for the experimental curves as for the numerical curves proposed by our model. One notices very particularly, as regards the tertiary creep, that the numerical curve is almost superimposed on the experimental curve. Moreover, the calculation of the relative error shows that for the two levels of stress, the relative error is less than 1%, which means that the model is faithful and is therefore able to describe the three phases of creep. In these figures, we also distinguish the curve proposed by the model of Wang et al. It is noted that the model proposed by the authors is well suited to model the primary and secondary creep, but it is difficult to describe the rapid growth of the strain during the tertiary creep. This difficulty has been underlined in their work [28], and can be corrected using the model that we propose.
The analysis of creep failure has mainly been focused on the determination of the final failure time t r . The final rupture time proposed by the model was determined and compared to the experimental rupture time.
The determination of the breaking time was made from the relation: where D c is the critical damage failure threshold, and t r is the final failure time. The theory of continuous damage mechanics [27] considers that failure occurs when the damage reaches a critical value D c = 1. However, it is observed experimentally [30] that the value of damage at failure for the materials is always lower than 1. Within the framework of our study, we consider that the rupture occurs for a value of the damage at least equal to 90% of the theoretical critical value which leads us to D c = 0.90. The results are shown in Tables 4 and Table 5. The calculation of the relative error relating to the estimation of the breaking time was made according to the relation: where t ex r is the rupture time obtained experimentally, and t m r is the rupture time obtained from the model. The analysis of the failure time shows that the time to failure obtained by the model is slightly less than the experimental rupture time. Furthermore, the calculation of the relative error related to the failure time shows that our model describes well the failure by creep.

Conclusion
A model has been proposed in order to describe the primary, secondary and tertiary creep as well as the failure of wood-based materials. The model is a rheological model in which the choice of components and assembly was made in accordance with experimental observations and in the interests of minimizing the number of parameters. The proposed model has been tested using the experimental results of the creep tests carried out by Wang [28] and it is shown to be very capable of describing the three phases of creep, with a relative error of less than 1%. Also, the analysis of the time proposed by the model showed that the selected critical damage D c = 0.90 considered as the breaking point damage led to a slightly lower breaking time (Err = 0.01) than the breaking time obtained experimentally. The model requires a total of five independent parameters all physically interpretable and easily accessible. The time to failure is easily accessible, thanks to the simplicity of our model, without necessarily going through heavy algorithms which are generally expensive in time and memory space. This represents a significant advantage of our model, which in sum offers both a more realistic way of describing the three phases of creep by fully accounting for the phenomenon of damage during the tertiary phase, and a simple and fast way to analyze the rupture time, compared to other models in the literature (Kaliske et al. [21], Wang et al. [16]). However, the model was developed in a one-dimensional case and assuming constant temperature and humidity conditions. It would be important to question the effects of humidity on temperature. A real challenge would also be that of taking into account non-constant stresses and the extension of the model to the threedimensional case. All these challenges open up perspectives for future work.