Creep Modelling of a Material by Non-Linear Modified Schapery’s Viscoelastic Model

Priska Kuida Atchounga, Grégoire Kamdjo, Emmanuel Foadieng, Pierre Kisito Talla

Laboratory of Mechanics and Modelling of Physical Systems (L2MSP), Department of Physics, University of Dschang, Dschang, Cameroon.

DOI: 10.4236/wjet.2017.54063                

Abstract

This research work aims at modeling the creep behavior of a material by a non-linear schapery’s viscoelastic model. We started with analytical part where three powerful methods of creep modeling have been developed and compared. That is the Heaviside, the Nordin and Varna and lastly our own proposed methods. From this preliminary study, it came out that our method is different to the two others because we took into account the loading time at the creep beginning. Besides we studied several loading programs and retained a five order non-linear polynomial which is the program that gave us satisfactory results. The other loading functions led to divergent results and wasn’t present here as consequence. In the second part of this work, we devoted ourselves to the determination of non-linear parameters in the schapery’s viscoelasticity equation, through a well developed and illustrated methodology. From this study, it is straight forward that non-linear parameters are stress dependent; confirming the results of several authors that preceded us in this studying field.

Keywords

Non-Linear Viscoelasticity, Creep, Strain, Stress, Schapery

1. Introduction

Creep is a physical phenomenon affecting many materials like woods, iron etc. in engineering structures like buildings. Before choosing a given material in engineering works, one must know very well its creep behavior in order to appreciate the lifespan of the structure. In such structures a deformation occurs when a material undergoes certain load. When it comes to study this phenomenon in the laboratory, we usually choose a normalized test material and with the help of a test machine we submit this material to a certain stress and follow back the deformation that occurs over the time. The ability to carry out reliable creep tests in a reasonable time at low stress levels allows a designer to have much more confidence in the data for creep-rupture behavior for materials and allows confident prediction of structural lifetimes. The inconvenience of this experimental method is that it can’t permit to follow the behavior of a material over a large period of time, so it is limited. In order to solve this problem, it is necessary to develop theoretical methods that can allow to model and to predict the creep behavior over a very large period of time in terms of years or even century.

Many authors [1] [2] [3] [4] [5] devoted their works to the creep modeling through experimental and theoretical methods. In this paper, we develop a theoretical method from the well known schapery’s equation for viscoelasticity. The main task consists of determining the non-linear parameters. We started by presenting the method where the loading and the unloading of the material are described by Heaviside step function. Because this method presents shortcomings [6] , it has been followed by the Nordin and Varna method and completed lastly by our analytical method.

2. Non-Linear Viscoelastic Material Model

The non-linear viscoelastic Schapery model [7] [8] is given by

$\epsilon \left(t\right)=g_0\left(\sigma \right)D_0\sigma \left(t\right)+g_1\left(\sigma \right){\displaystyle {\int }_0^t\Delta D\left(\psi -{\psi }^{\prime }\right)}\frac{\text{d}\left(g_2\left(\sigma \right)\sigma \left(t\right)\right)}{\text{d}\tau }\text{d}\tau$ (1)

where the reduced times are given by

$\psi ={\displaystyle {\int }_0^t\frac{\text{d}{t}^{\prime }}{a_{\sigma }\left(\sigma \left(t^{\prime }\right)\right)}}$ (2)

and

${\psi }^{\prime }={\displaystyle {\int }_0^{\tau }\frac{\text{d}{t}^{\prime }}{a_{\sigma }\left(\sigma \left(t^{\prime }\right)\right)}}$ (3)

where a_σ is a shift factor. The parameters g₀, g₁, g₂ and a_σ are functions of strain. D₀ = D(0) is the initial value of creep compliance and $\Delta D=D\left(t\right)-D_0$ is the transient component of the creep compliance. When $a_{\sigma }=g_0=g_1=g_2=1$ , Equation (1) reduces to

$\epsilon \left(t\right)=D_0\sigma \left(t\right)+{\displaystyle {\int }_0^t\Delta D\left(t-\tau \right)}\frac{\text{d}\sigma \left(t\right)}{\text{d}\tau }\text{d}\tau$ (4)

Or

$\epsilon \left(t\right)={\displaystyle {\int }_0^{t}D\left(t-\tau \right)}\frac{\text{d}\sigma \left(t\right)}{\text{d}\tau }\text{d}\tau$ (5)

which is the Boltzmann’s superposition integral for linear viscoelasticity. In order to ensure linear viscoelastic behavior at small stresses, the following initial values must hold:

$a_{\sigma }\left(0\right)=g_0\left(0\right)=g_1\left(0\right)=g_2\left(0\right)=1$ (6)

Many methods have been developed to determine the material parameters in Equation (1), see [6] - [11] .

3. Methods of Analysis

3.1. Step-Stress Hypothesis

Under step-stress hypothesis, that is, $\sigma \left(t\right)={\sigma }_{0}H\left(t\right)$ , where H(t) is the Heaviside step function, Equation (1) takes the form

$\epsilon \left(t\right)=g_0\left({\sigma }_0\right)D_0{\sigma }_0+g_1\left({\sigma }_0\right){\displaystyle {\int }_0^t\Delta D\left(\psi -{\psi }^{\prime }\right)}\frac{\text{d}\left(g_2\left({\sigma }_0\right){\sigma }_{0}H\left(\tau \right)\right)}{\text{d}\tau }\text{d}\tau$ (7)

$\epsilon \left(t\right)=g_0\left({\sigma }_0\right)D_0{\sigma }_0+g_1\left({\sigma }_0\right){\displaystyle {\int }_0^t\Delta D\left(\psi -{\psi }^{\prime }\right)g_2\left({\sigma }_0\right){\sigma }_0\delta \left(\tau \right)\text{d}\tau }$ (8)

$\epsilon \left(t\right)=g_0\left({\sigma }_0\right)D_0{\sigma }_0+g_1\left({\sigma }_0\right)g_2\left({\sigma }_0\right){\sigma }_0\Delta D\left(\frac{t}{a_{\sigma }\left({\sigma }_0\right)}\right)$ (9)

where δ is the Dirac delta function. Equation (9) is the material response when the stress is applied by the Heaviside step function.

3.2. Method by Nordin and Varna

Let the stress be given by

$\sigma \left(t\right)=\frac{{\sigma }_0}{2}\left[H\left(t\right)+H\left(t-t_1\right)\right]$ (10)

Then

$g_2\left(\sigma \right)=g_2\left(\frac{{\sigma }_0}{2}\right)H\left(t\right)+\left[g_2\left({\sigma }_0\right)-g_2\left(\frac{{\sigma }_0}{2}\right)\right]H\left(t-t_1\right)$ (11)

and

$g_2\left(\sigma \right)\sigma \left(t\right)=\frac{{\sigma }_0}{2}\left[g_2\left(\frac{{\sigma }_0}{2}\right)H\left(t\right)+\left(2g_2\left({\sigma }_0\right)-g_2\left(\frac{{\sigma }_0}{2}\right)\right)H\left(t-t_1\right)\right]$ (12)

Differentiating Equation (12) with respect to time gives

$\frac{\partial \left(g_2\left(\sigma \right)\sigma \left(t\right)\right)}{\partial t}=\frac{{\sigma }_0}{2}\left[g_2\left(\frac{{\sigma }_0}{2}\right)\delta \left(t\right)+\left(2g_2\left({\sigma }_0\right)-g_2\left(\frac{{\sigma }_0}{2}\right)\right)\delta \left(t-t_1\right)\right]$ (13)

Substituting Equation (13) in Equation (1) when $t\ge t_1$ gives

$\begin{array}{l}\epsilon \left(t\right)=g_0\left({\sigma }_0\right)D_0{\sigma }_0+\frac{1}{2}{g}_1\left({\sigma }_0\right){\sigma }_0{g}_2\left(\frac{{\sigma }_0}{2}\right){\displaystyle {\int }_0^t\Delta D\left(\psi -{\psi }^{\prime }\right)}\\ +\frac{1}{2}{g}_1\left({\sigma }_0\right){\sigma }_0\left(2g_2\left({\sigma }_0\right)-g_2\left(\frac{{\sigma }_0}{2}\right)\right){\displaystyle {\int }_0^t\Delta D\left(\psi -{\psi }^{\prime }\right)\delta \left(t-t_1\right)\text{d}\tau }\end{array}$ (14)

After integrating with mathematical formulae Equation (14) yields

$\begin{array}{l}\epsilon \left(t\right)=g_0\left({\sigma }_0\right)D_0{\sigma }_0+\frac{1}{2}{g}_1\left({\sigma }_0\right){\sigma }_0{g}_2\left(\frac{{\sigma }_0}{2}\right)\Delta D\left(\frac{t_1}{a_{\sigma }\left({\sigma }_0/2\right)}+\frac{t-t_1}{a_{\sigma }\left({\sigma }_0\right)}\right)\\ \text{+}\frac{\text{1}}{\text{2}}{g}_1\left({\sigma }_0\right){\sigma }_0\left[2g_2\left({\sigma }_0\right)-g_2\left(\frac{{\sigma }_0}{2}\right)\right]\Delta D\left(\frac{t-t_1}{a_{\sigma }\left({\sigma }_0\right)}\right)\end{array}$ (15)

which represents the material response under a tress defined by Equation (10).

3.3. Proposed Method

We now consider case where the stress is given by

$\sigma \left(t\right)=\{\begin{array}{l}f\left(t\right),t\prec t_1\\ {\sigma }_0,t\ge t_1\end{array}$ (16)

where $f\left(0\right)=0$ and $f\left(t_1\right)={\sigma }_0$ . Then the strain at time $t\ge t_1$ is given by

$\epsilon \left(t\right)=g_0\left({\sigma }_0\right)D_0{\sigma }_0+g_1\left({\sigma }_0\right){\displaystyle {\int }_0^{t_1}\Delta D\left(\psi -{\psi }^{\prime }\right)\frac{\text{d}\left(g_2\left(\sigma \right)\sigma \right)}{\text{d}\tau }\text{d}\tau }$ (17)

where

$\begin{array}{c}\psi -{\psi }^{\prime }={\displaystyle {\int }_0^t\frac{\text{d}{t}^{\prime }}{a_{\sigma }\left({\sigma }_0\right)}}-{\displaystyle {\int }_0^{\tau }\frac{\text{d}{t}^{\prime }}{a_{\sigma }\left({\sigma }_0\right)}}\\ ={\displaystyle {\int }_0^{t_1}\frac{\text{d}{t}^{\prime }}{a_{\sigma }\left({\sigma }_0\right)}}+{\displaystyle {\int }_{t_1}^t\frac{\text{d}{t}^{\prime }}{a_{\sigma }\left({\sigma }_0\right)}-\left({\displaystyle {\int }_0^{t_1}\frac{\text{d}{t}^{\prime }}{a_{\sigma }\left({\sigma }_0\right)}+{\displaystyle {\int }_{t_1}^{\tau }\frac{\text{d}{t}^{\prime }}{a_{\sigma }\left({\sigma }_0\right)}}}\right)}\end{array}$ (18)

$\psi -{\psi }^{\prime }=\frac{t-t_1}{a_{\sigma }\left({\sigma }_0\right)}+{\displaystyle {\int }_{\tau }^{t_1}\frac{\text{d}{t}^{\prime }}{a_{\sigma }\left({\sigma }_0\right)}}$ (19)

By combining Equation (19) and Equation (17), we have

$\begin{array}{l}\epsilon \left(t\right)=g_0\left({\sigma }_0\right)D_0{\sigma }_0\\ +g_1\left({\sigma }_0\right)t_1\Delta D\left(\frac{t-t_1}{a_{\sigma }\left({\sigma }_0\right)}+{\displaystyle {\int }_{\tau }^{t_1}\frac{\text{d}{t}^{\prime }}{a_{\sigma }\left(\sigma \left(t^{\prime }\right)\right)}}\right)\frac{\text{d}\left(g_2\left(\sigma \left(\tau \right)\right)\sigma \left(\tau \right)\right)}{\text{d}\tau }\end{array}$ (20)

Using midpoint rule, which is third-order accurate with respect to t₁ [12] , and with τ = t₁/2 then it follows that

$\begin{array}{c}{\displaystyle {\int }_{\tau }^{t_1}\frac{\text{d}{t}^{\prime }}{a_{\sigma }\left(\sigma \left(t^{\prime }\right)\right)}}={\displaystyle {\int }_{t_1/2}^{t_1}\frac{\text{d}{t}^{\prime }}{a_{\sigma }\left(\sigma \left(t^{\prime }\right)\right)}}=\frac{\left(t_1-t_1/2\right)}{a_{\sigma }\left(\sigma \left(\frac{t_1+t_1/2}{2}\right)\right)}\\ =\frac{t_1/2}{a_{\sigma }\left(\sigma \left(\frac{3t_1}{4}\right)\right)}=\frac{t_1/2}{a_{\sigma }\left(f\left(\frac{3t_1}{4}\right)\right)}\end{array}$ (21)

If the ramp loading is approximated to be linear, then $f\left(3t_1/4\right)=3{\sigma }_0/4$ . So we have

${\displaystyle {\int }_{\tau }^{t_1}\frac{\text{d}{t}^{\prime }}{a_{\sigma }\left(\sigma \left(t^{\prime }\right)\right)}=\frac{t_1/2}{a_{\sigma }\left(3{\sigma }_0/4\right)}}$ (22)

Let us apply second-order [13] accurate numerical differentiation formula to the function in Equation (20), it comes that

${\left[\frac{\text{d}\left(g_2\left(\sigma \left(\tau \right)\right)\sigma \left(\tau \right)\right)}{\text{d}\tau }\right]}_{\tau =t_1/2}\approx \frac{g_2\left(\sigma \left(t_1\right)\right)\sigma \left(t_1\right)-g_2\left(\sigma \left(0\right)\right)\sigma \left(0\right)}{t_1}=\frac{g_2\left({\sigma }_0\right){\sigma }_0}{t_1}$ (23)

Now substituting Equation (22) and Equation (23) in Equation (20), it comes out that the response of the material when the applied stress is defined by Equation (16) is

$\epsilon \left(t\right)=g_0\left(\sigma {}_0\right)D_0{\sigma }_0+g_1\left({\sigma }_0\right)g_2\left({\sigma }_0\right){\sigma }_0\Delta D\left(\frac{t}{a_{\sigma }\left({\sigma }_0\right)}+\Omega \right)$ (24)

where

$\Omega =t_1\left(\frac{1}{2a_{\sigma }\left(3{\sigma }_0/4\right)}-\frac{1}{a_{\sigma }\left({\sigma }_0\right)}\right)$ (25)

Only a rather moderate modification is made in the proposed method as compared to the step loading case, where Ω = 0. Furthermore, in linear case the proposed correction method yields

$\epsilon \left(t\right)={\sigma }_{0}D\left(t-\frac{t_1}{2}\right)$ (26)

which is the Zapas-Phillips [11] [14] correction method for linear viscoelastic systems.

4. Numerical Studies

In the following subsections we compare the accuracy of the Heaviside step loading, Nordin-Varna and the proposed method.

4.1. Linear Case

In linear case non-linear parameters in Equation (1) equals to one. Then the Nordin-Varna method for the creep formulation at time t ≥ t₁ is given by

$\epsilon \left(t\right)=D_0{\sigma }_0+\frac{{\sigma }_0}{2}\Delta D\left(t\right)+\frac{{\sigma }_0}{2}\Delta D\left(t-t_1\right)$ (27)

$\epsilon \left(t\right)=\frac{{\sigma }_0}{2}\left[D\left(t\right)+D\left(t-t_1\right)\right]$ (28)

and with the proposed method it is given by

$\epsilon \left(t\right)=D_0{\sigma }_0+{\sigma }_0\Delta D\left(t-\frac{t_1}{2}\right)={\sigma }_{0}D\left(t-\frac{t_1}{2}\right)$ (29)

The exact value for the strain at time t ≥ t₁ is

$\epsilon \left(t\right)={\displaystyle {\int }_0^{t_1}D\left(t-\tau \right)\stackrel{\dot{}}{\sigma }\left(\tau \right)\text{d}\tau }$ (30)

If the loading is carried out with a constant stress rate, then

$\epsilon \left(t\right)=\stackrel{\dot{}}{\sigma }{\displaystyle {\int }_0^{t_1}D\left(t-\tau \right)\text{d}\tau }$ (31)

where $\stackrel{\dot{}}{\sigma }={\sigma }_0/t_1$ . Now it can be seen that by applying numerical integration rule, to Equation (31) we get the Nordin-Varna method and by applying midpoint rule we get the proposed method, respectively. Error estimate for the trapezoidal rule is [13]

$e=-\frac{t_1^3}{12}{D}^{″}$ (32)

and for the midpoint rule [13]

$e=\frac{t_1^3}{24}{D}^{″}$ (33)

To summarize, in linear case error in the proposed method is half of error in the Nordin-Varna method.

4.2. Non-Linear Case, Creep Test with Non-Linear Ramp

The creep test that we study has the following form

$\epsilon \left(t\right)=\{\begin{array}{l}f\left(t\right),t\prec t_1\\ {\sigma }_0,t\ge t_1\end{array}$ (34)

where σ₀ = 42.26 MPa, t₁ = 4 s and

$f\left(t\right)=-\frac{1}{6}{\left(\frac{2t-t_1}{t_1}\right)}^5{\sigma }_0+\frac{3}{2}\frac{t}{t_1}{\sigma }_0+\frac{1}{4}{\sigma }_0$ (35)

is the loading function, see Figure 1 below.

The transient component of the creep compliance is taken to be

$\Delta D=\alpha t^n{\text{MPa}}^{-1},\text{0}\prec n\prec 1$ (36)

5. Parameter Identification

Material parameters are determined as follows:

1) With our method the strain at time t ≥ t₁ is now given by

$\epsilon \left(t\right)=g_0\left({\sigma }_0\right)D_0{\sigma }_0+g_1\left({\sigma }_0\right)g_2\left({\sigma }_0\right)\alpha {\left(\frac{t}{a_{\sigma }\left({\sigma }_0\right)}+\Omega \right)}^n{\sigma }_0$ (37)

$\epsilon \left(t\right)=g_0\left({\sigma }_0\right)D_0{\sigma }_0+\frac{g_1\left({\sigma }_0\right)g_2\left({\sigma }_0\right)}{a_{\sigma }^n\left({\sigma }_0\right)}\alpha {\left[a_{\sigma }t+\Omega a_{\sigma }\right]}^n{\sigma }_0$ (38)

where

Figure 1. Non-linear ramp loading.

$\Omega =t_1\left[\frac{1}{2a_{\sigma }\left(3{\sigma }_0/4\right)}-\frac{1}{a_{\sigma }\left({\sigma }_0\right)}\right]$ (39)

2) Creep test data with different constant stress levels is fitted to

$f\left(t\right)=A+B{\left[a_{\sigma }\left({\sigma }_0\right)t+\Omega a_{\sigma }\left({\sigma }_0\right)\right]}^n$ (40)

where

$\{\begin{array}{l}f\left(t\right)=\epsilon \left(t\right)/{\sigma }_0\\ A=g_0\left({\sigma }_0\right)D_0\\ B=\frac{g_1\left({\sigma }_0\right)g_2\left({\sigma }_0\right)}{a_{\sigma }^n\left({\sigma }_0\right)}\alpha \end{array}$ (41)

3) Now the values of A(σ₀) and B(σ₀) are known for all constant stress levels.

4) Parameter A is fitted to some proper function A(σ). We have used five order polynomial to approximate A(σ), in this study. Then the parameter D₀ can be determined using initial condition $A\left(0\right)=g_0\left(0\right)D_0=D_0$ . Since $g_0\left(\sigma \right)=A\left(\sigma \right)/D_0$ the parameter g₀ can be also determined [15] .

5) Parameter B is fitted to some proper function B(σ). In this study, we have used five order polynomial to approximate B(σ). Then the parameter α can be determined using initial condition $B\left(0\right)=\frac{g_1\left(0\right)g_2\left(0\right)}{a_{\sigma }^n\left(0\right)}\alpha =\alpha$ .

6) We have $\frac{g_1\left({\sigma }_0\right)g_2\left({\sigma }_0\right)}{a_{\sigma }\left({\sigma }_0\right)}=\frac{B\left({\sigma }_0\right)}{\alpha }$ . The value of B(σ₀)/α is known for all

constant stress levels, this value is denoted C. Parameter C is fitted to some proper function C [16] . We have used five order polynomial to approximate C(σ). Then, parameters g₁, g₂ and a_σ can be determined using initial conditions $g_1\left(0\right)=g_2\left(0\right)=a_{\sigma }\left(0\right)=1$ .

7) Since $a_{\sigma }^n\left({\sigma }_0\right)=\frac{g_1\left({\sigma }_0\right)g_0\left({\sigma }_0\right)}{\alpha B\left({\sigma }_0\right)}$ then we have

$n=\frac{\ln \left(\left(g_1\left({\sigma }_0\right)g_2\left({\sigma }_0\right)\right)/\left(\alpha B\left({\sigma }_0\right)\right)\right)}{\ln \left(a_{\sigma }\left({\sigma }_0\right)\right)}$

The preceding methodology leads us to the following values:

$D_0=160\times {10}^{-6}{\text{MPa}}^{-1};\alpha =1.03\times {10}^{-6}{\text{MPa}}^{-1}\text{and}n=0.24$ (42)

Simulated non-linear parameters are

$\{\begin{array}{l}{g}_0\left(\sigma \right)=3.2\times {10}^{-7}{\sigma }^5-5.7\times {10}^{-5}{\sigma }^4+3.9\times {10}^{-3}{\sigma }^3-0.13{\sigma }^2+2\sigma -11\\ g_1\left(\sigma \right)=8.5\times {10}^{-8}{\sigma }^5-1.3\times {10}^{-5}{\sigma }^4+8.1\times {10}^{-4}{\sigma }^3-2.4\times {10}^{-2}{\sigma }^2+0.34\sigma -0.83\\ g_2\left(\sigma \right)=4.9\times {10}^{-7}{\sigma }^5-8.5\times {10}^{-5}{\sigma }^4+5.7\times {10}^{-3}{\sigma }^3-0.19{\sigma }^2+2.9\sigma -17\\ a_{\sigma }\left(\sigma \right)=5.7\times {10}^{-7}{\sigma }^5-9.7\times {10}^{-5}{\sigma }^4+6.4\times {10}^{-3}{\sigma }^3-0.21{\sigma }^2+3.2\sigma -18\end{array}$ (43)

6. Relative Errors and Creep Curves

We computed the relative error with the following formula

$e\left(t\right)=\frac{\Vert \epsilon \left(t\right)-\bar{\epsilon }\left(t\right)\Vert }{\Vert \epsilon \left(t\right)\Vert }$ (44)

where $\epsilon \left(t\right)$ the true value of the strain in the creep test and $\bar{\epsilon }\left(t\right)$ is the numerically approximated value of strain. The Table 1 below depicts relative errors:

Error estimates shows that the proposed method produces the smallest error. The computed strain in creep test with different methods is shown in Figure 2 below.

Non-linear parameters are shown in Figure 3.

Figure 3 depicts the predicted non-linear parameters (straight curve) with our proposed method and the true value of non-linear parameters (disconti- nuous curve with “+” sign). From these figures it is evident that the simulated non-linear parameters are in good agreement with the true value of the non-li- near parameters. We can notice from our results that stress highly influences the value of material non-linear parameters; which is a confirmation that parameters g₀, g₁, g₂ and a_σ in the Schapery’s viscoelasticity equation are stress dependent [7] [8] . These results are matching perfectly with those obtained by [1] [2] [3] [4] when they were dealing with physical and mechanical properties of some Cameroonians woods. Authors [17] [18] [19] when dealing with non-linear creep and relaxation obtained similar results in their scientific works.

Figure 2 is just depicting the advantage of predicting creep behavior of material with our method. It is clear that the predicted creep curve is so close to the true value of the creep to be distinguished. The lower value of the relative estimating errors is just reinforcing the method.

Figure 2. Strain as a function of time in creep test.

Table 1. Relative errors in creep test.

Figure 3. Non-linear parameters as a function of stress.

7. Conclusion

We have presented in this paper a powerful method which takes in to account the finite ramp time in the Schapery’s non-linear viscoelastic equation. It came out that the method is a good predicting tool of strain in creep test, because it is doing while minimizing the relative error. At the end material non-linearity parameters that have been simulated from the proposed method are in good agreement with those found in literature. The authors [20] , [21] and [16] applied different processes in their works to predict long term creep of composites and they came out with good results. In the future we can see how to apply our correction method to what they did in order to have different point of view.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Talla, P.K., Foadieng, E., Fouotsa, W.C.M., Fogue, M., Bishweka, S., Ngarguededjim K.E., Alabeweh, F.S. and Foudjet, A. (2015) A Contribution to the Study of Entandrophragma Cylindricum Sprague and Lovoa Trichilioides Harms Long Term behaviour. Revue scientifique et Technique Forêt et Environnement du Bassin du Congo, X, 10-21.
[2]	Talla, P.K., Mabekou, J.S., Fogue, M., Fomethe, A., Bawe, G.N., Foadieng, E. and Foudjet, A. (2010) Non-Linear Creep Behavior of Raphia vinifera L. Arecacea under Flexural Load. International Journal of Mechanics and Solids, 5, 151-172.
[3]	Talla, P.K. (2008) Contribution à l’analyse mécanique de Raphia vinifera L. Arecacea. Thèse de Doctorat (Ph.D), Université de Dschang, Faculté des Sciences, Cameroun.
[4]	Talla, P.K., Pelab, F.B., Fogue, M., Fomethe, A., Bawe, G.N., Foadieng, E. and Foudjet, A. (2007) Non-Linear Creep Behavior of Raphia vinifera L. Arecacea. International Journal of Mechanics and Solids, 2, 1-11.
[5]	Foadieng, E., Fogue, M. and Talla, P.K. (2012) Effect of the Span Length on the Deflection and the Creep Behavior of Raffia Bamboo Vinifera L. Arecacea Beam. International Journal of Material Science, 7, 153-167.
[6]	Nordin, L.O. and Varna, J. (2006) Methodology for Parameter Identification in Non-Linear Viscoelastic Material Model. Mech. Time Matter, 9, 259-280.
[7]	Schapery, R.A. (1969) On the Characterization of Non-Linear Viscoelastic Materials. Polymer Engineering & Science, 9, 295-310. https://doi.org/10.1002/pen.760090410
[8]	Lou, Y.C. and Schapery, R.A. (1971) Viscoelastic Characterization of a Non-Linear Fiber-Reinforced Plastic. Journal of Composite Materials, 5, 208-234. https://doi.org/10.1177/002199837100500206
[9]	Kshitish, A.P. (2009) Linear and Non-Linear Viscoelastic Characterization of Proton Exchange Membranes and Stress Modeling for Fuel Cell Application. Doctor of Philosophy thesis, Virginia Polytechnic Institute and State University, USA.
[10]	Chien, W.H., Rashid, K.A., Eyad, A.M., Dallas, N.L. and Gordon, D.A. (2011) Numerical Implementation and Validation of a Non-Linear Viscoelastic and Viscoplastic Model for Asphalt Mixes. International Journal of Pavement Engineering, 12, 433-447. https://doi.org/10.1080/10298436.2011.574137
[11]	Zapas, L.J. and Philips, J.C. (1971) Simple Shearing Flows in Polyisobutylene Solutions. Journal of Research of the National Bureau of Standards, 75A, 33-41. https://doi.org/10.6028/jres.075A.005
[12]	Bakhvalov, N.S. (1977) Numerical Methods: Analysis, Algebra, Ordinary Differential Equations. MIR, Moscow.
[13]	Rami, M.H.K. and Anastasia, H.M. (2003) Numerical Finite Element Formulation of the Schapery Non-Linear Viscoelastic Material Model. International Journal for Numerical Methods in Engineering, 59, 25-45.
[14]	Ahlberg, J.H., Nilson, E.N. and Walsh, J.L. (1967) The Theory of Splines and Their Applications. Academic Press, New York.
[15]	Benjamin, B. (2013) Interpolation Polynomiale. Agrégation de mathématiques, Option modélisation, Université de Rennes 1.
[16]	Ratchada, S. and Raffaella, D.V. (2011) A Mathematical Model for Creep, Relaxation and Strain Stiffening in Parallel-Fibered Collagenous Tissues. Journal of Medical Engineering and Physics, 33, 1056-1063. https://doi.org/10.1016/j.medengphy.2011.04.012
[17]	Calgcagno, B., Lopez, G.M., Kuhns, M. and Lakes, R.S. (2008) On the Non-Linear Creep and Recovery of Open Cell Earplug Foams. Journal of Cellular Polymers, 27, 165-178.
[18]	Kaouther, B.A.A. (2010) Relations entre propriétés rhéologiques et structure microscopique de dispersions de particules d’argile dans des solutions de polymères. Thèse de Doctorat (Ph.D.), Université de Haute Alsace, France.
[19]	Ashish, O., Ray, V.J. and Roderic, S.L. (2003) Interralation of Creep and Relaxation for Non-Linearly Viscoelastic Materials: Application to Ligament and Metal. Rheologica Acta, 42, 557-568. https://doi.org/10.1007/s00397-003-0312-0
[20]	Tuttle, M.E. and Brinson, H.F. (1985) Prediction of the Long-Term Creep Compliance of General Composite Laminates. Experimental Mechanics, 26, 89-102. https://doi.org/10.1007/BF02319961
[21]	Ever, J.B. (2007) Prediction of Long-Term Creep of Composites from Doubly-Shifted Polymer Creep Data. Journal of Composite Materials.