A generalised fractional derivative model to represent elastoplastic behaviour of metals

International Journal of Mechanical Sciences 65 (2012) 12–17

Contents lists available at SciVerse ScienceDirect

International Journal of Mechanical Sciences journal homepage: www.elsevier.com/locate/ijmecsci

A generalised fractional derivative model to represent elastoplastic behaviour of metals Joseba Mendiguren a,n, Fernando Corte´s b,1, Lander Galdos a,2 a b

´n, Spain Department of Mechanical Engineering, Mondragon Unibertsitatea, Loramendi 4, 20500 Mondrago Deusto Institute of Technology (DeustoTech), Faculty of Engineering, University of Deusto, Avda. de las Universidades 24, 48007 Bilbao, Spain

a r t i c l e i n f o

a b s t r a c t

Article history: Received 16 January 2012 Received in revised form 8 July 2012 Accepted 21 August 2012 Available online 30 August 2012

In this paper, a one-dimensional elastoplastic model based on fractional calculus is presented. This model is conceived for metal deformation. First, a brief introduction about fractional calculus is exposed. Then, a fractional calculus based general model is proposed, build-up as a series of fractional derivatives. The particular cases of one, two and three terms are discussed. Finally, these three particular cases of the general model are fitted to experimental data and compared with the classical Hollomon and Ramberg–Osgood models. As a result, a model capable of representing the elastoplastic behaviour more precisely than the Hollomon and Ramberg–Osgood models is obtained, despite numerical efforts are enlarged. & 2012 Elsevier Ltd. All rights reserved.

Keywords: Fractional calculus Model behaviour Elastoplasticity Numerical methods Plasticity

1. Introduction Nowadays, automotive industry is becoming more and more competitive, car safety has to augment while vehicles must be lighter to reduce fuel consumption, reducing prices to be more competitive in the market. Due to this increment in the competitiveness of the market, the automotive engineering needs more accurate material models in order to improve the numerical simulations needed to reach competitive designs. Among all manufacturing processes used in automotive industry, metal plastic forming is one of the most important. This kind of process exploits the fact that when a metallic material is deformed enough to reach its elastoplastic region, the material remains with this new configuration. In fact, when a material is deformed, it has two regions: a first region of elastic behaviour until reach certain stress value, and when this value of stress is reached, a second elastoplastic region starts, in which permanent deformation occurs. Materials are usually characterised by tensile test. In this test, the material is stretched in one direction until rupture and elastic and plastic properties of the material are obtained [1]. The elastoplastic

n

Corresponding author. Tel.: þ34 943794700; fax: þ 34 943791536. E-mail addresses: [email protected] (J. Mendiguren), [email protected] (F. Corte´s), [email protected] (L. Galdos). 1 Tel.: þ34 944139003. 2 Tel.: þ34 943794700; fax: þ 34 943791536. 0020-7403/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijmecsci.2012.08.008

behaviour of the material can be represented from a two different point of view. On the one hand, the industrial representation where the objective is to represent the whole behaviour until rupture and even extended to a high strain ranges. On the other hand, from a computational plasticity point of view where the curve is split and the elastic behaviour is represented by means of the elastic modulus and the plastic behaviour is represented with the hardening curve. In the last years, the representation of the plastic behaviour from a computational plasticity point of view has been deeply developed. New combined hardening models are introduced and new recovering techniques are developed [2–5]. However, the correct representation of these two behaviours (elastic and plastic behaviours) together in a single model is still a challenge. A single model able to describe both behaviours more accurately than current models is a high research due to the exigences of the market. The use of this kind of representations is widely used in industry in order to compare the behaviour of two materials and to extend the tensile test data to get a piecewise data up to the tensile test rupture strain. That is why, the accuracy of the single model and specifically just before necking has become a research interest from an industrial point of view. Among the different models that can be found in the literature to describe elastoplasticity, both Hollomon and Ramberg–Osgood models are the most used. On the one hand, the Hollomon model satisfies

s ¼ K enH ,

ð1Þ

where K and nH are the model parameters, s and e being the stress and strain, respectively. On the other hand, the Ramberg–Osgood

J. Mendiguren et al. / International Journal of Mechanical Sciences 65 (2012) 12–17

model is given by s  s nRO e¼ þ , E H

ð2Þ

where E, H, and nRO are the model parameters, the E parameter represents the elastic modulus of the material. Both models are used to represent the whole strain–stress curve metals. For more details about these models, see [6,7,1]. Elastoplastic behaviour is composed, on the one hand, of elastic behaviour usually represented in metals by the Hooke law

s ¼ Ee,

ð3Þ

where E is the Young modulus of the material, and on the other hand, of plastic behaviour, which can be represented in many cases as a perfect plastic law, which satisfies

s ¼ s0 ,

ð4Þ

where s0 is the yield strength. Therefore, writing these behaviours as differential equations, the elastic behaviour exposed in Eq. (3) can be represented by 1

d s ¼ E, de1

ð5Þ

and the perfect plastic behaviour by 0

d s ¼ s0 : de0

ð6Þ

Consequently, the main hypothesis for the construction of the proposed model is that elastoplasticity can be modelled by means of a derivative order a comprised between zero and one a

d s ¼C dea

ð0 o a o 1Þ:

ð7Þ

The viability of using fractional calculus to represent the elastoplastic behaviour of metals stems from the application of fractional models to represent viscoelasticity (some of this applications can be seen in [8–11], among others). In fact, it can be taken into account that elastic behaviour is represented as the zero order derivative of the strain with respect to time, according to Hooke’s law (3) and viscosity is represented as the first order derivative with respect time t, according to Newton’s law

s ¼ ce_ ðtÞ,

ð8Þ

and then, viscoelasticity, which is a combination of elastic and viscous behaviours, may be represented by means of an a order derivative, satisfying 0 o a o 1 a

s¼a

d e : dt a

ð9Þ

In short, this paper leads with the capability of a fractional calculus based model, which is named fractional model hereafter, to represent the elastoplastic behaviour of metals

 First, a brief introduction about fractional calculus is presented.  Next, the new fractional model is proposed, analysing the 

behaviour of the one, two and three term particular cases, both analytically and numerically. Finally, particular cases of the fractional model are fitted to experimental data in order to compare them with the classical Hollomon and Ramberg–Osgood models.

Historically, fractional calculus has been relegated to theoretical mathematics, but in last decades, it has been successfully applied in different physical phenomena modelling, as viscoelastic behaviour. Miller and Ross [12] have published an introductory book to the fractional calculus. In the same manner, some years later Podlubny [13] has exposed in his book a brief review of fractional calculus applications, e.g. fractional capacitor theory, biology, control theory. More recently, different authors have published more complex fractional models [9,14] about different physical phenomena. First, definitions of fractional integral and fractional derivative proposed by Riemann–Liouville are presented. Next, the defini¨ tion of fractional derivative proposed by Grunwald–Letnikov is analysed. The definition of fractional a order integral of f(t) function defined by Riemann–Lioville is given by t0 Dat f ðtÞ ¼

1 GðaÞ

Z

t

ðttÞa1 f ðtÞ dt,

ð10Þ

t0

where t0 is the lower integration limit, a is comprised between zero and one 0 o a o 1, and GðaÞ represents the gamma function of the a argument, which can be calculated as Z 1 GðaÞ ¼ et ta1 dt: ð11Þ 0

Representing the m order derivative as Dm, m A R þ , the fractional a order derivative Riemann–Lioville’s definition yields m

a f ðtÞ, t0 Dat f ðtÞ ¼ d t0 Im t

ð12Þ

which implies t0 Dat f ðtÞ ¼

" # Z t m d 1 f ðtÞ d t , dt m Gðmt 0 Þ t0 ðttÞam þ 1

m1 o a o m, ð13Þ

whose Laplace transform satisfies Lft0 Dat f ðtÞg ¼ sa FðsÞ

m 1 X

sk f

ak1

ðt 0 Þ,

ð14Þ

k¼0 ðÞ

where f ðt 0 Þ is the value of the ðÞ order derivative at the integration lower limit t0. In the material behaviour modelling field [9,12–14] the derivative order is commonly comprised between zero and one, 0 o a o 1. When this particular case is given, the general fractional derivative definition (13) can be rewritten as Z t  1 d f ðtÞ ð15Þ t0 Dat f ðtÞ ¼ a dt , 0 o a o 1: Gð1aÞ dt t0 ðttÞ Some authors have developed numerical algorithms to evaluate these fractional derivatives, for example [15–19] among others. One of the most efficient is the G1 method, [15, p. 136], which is based on the fractional derivative definition given by ¨ Grunwald–Letnikov (see [13] or [15] for more detailed analysis about this definition). ¨ The fractional derivative definition given by Grunwald– Letnikov is based on the generalization of the classical derivative. 0 In fact, the first derivative f ðtÞ is given by 0

f ðtÞ ¼ lim

Dt-0

2. Some theoretical aspects about fractional calculus In this section, in order to establish a baseline of knowledge, a brief review of some aspects about fractional calculus is presented.

13

f ðtÞf ðtDtÞ , Dt

ð16Þ

the second order derivative yields 00

f ðtÞ ¼ lim

Dt-0

f ðtÞ2f ðtDtÞ þf ðt2DtÞ ðDtÞ2

,

ð17Þ

14

J. Mendiguren et al. / International Journal of Mechanical Sciences 65 (2012) 12–17

series of fractional derivatives

the third order derivative results in 000

f ðtÞ ¼ lim

f ðtÞ3f ðtDtÞ þ 3f ðt2DtÞf ðt3DtÞ ðDtÞ3

Dt-0

,

ð18Þ

j¼0

where t

Dt

,

aq Daq sðeÞ ¼ 1,

where Daq sðeÞ represents the aq order derivative of stress s with a respect to strain e, and which is equivalent to e0 De q sðeÞ when the initial condition sð0Þ ¼ 0 is considered. Three different particular cases are analysed in this work: a one-term model with qmax ¼ 1, a two-term model with qmax ¼ 2, and a three-term model with qmax ¼ 3.

ð20Þ

where ðnjÞ represents Newton’s binomial ! n n! : ¼ j j!ðnjÞ!

3.1. One-term particular case model ð21Þ

When a unique term is taken into account, the model is given by

This expression can be simplified using gamma function according to ! ! n jn1 GðjnÞ j , ð22Þ ð1Þ ¼ ¼ j j GðnÞGðj1Þ

a1 Da1 sðeÞ ¼ 1:

thus Eq. (19) can be rewritten as 0 1 N1 X GðjnÞ ðnÞ n @ f ðtÞ ¼ lim ðDtÞ f ðtjDtÞA: GðnÞGðj1Þ Dt-0 j¼0

sðeÞ ¼ ð23Þ

Therefore, the definition of fractional derivative proposed by ¨ Grunwald–Letnikov is obtained when the derivative order in Eq. (23) is replaced by a real value, and a derivative lower limit t0, which satisfies N¼

tt 0 , Dt

ð24Þ

is defined. Using the same nomenclature as used for Riemann– Liouville, the definition of fractional order derivative proposed by ¨ Grunwald–Letnikov can be expressed as 0 1 N 1 X G ðj a Þ a  a t0 Dt f ðtÞ ¼ lim @ðDtÞ ð25Þ f ðtjDtÞA: GðaÞGðj1Þ Dt-0 j¼0 The G1 numerical method leads from Eq. (25) with a finite Dt discrete time step, resulting in 1 1 NX a Aa,j þ 1 f nj , a Dt f n  a ðDtÞ j ¼ 0

ð26Þ

¨ where Aa,j þ 1 are the a order Grunwald coefficients, defined as

GðjaÞ : Aa,j þ 1 ¼ GðaÞGðj þ 1Þ

ð27Þ

In this work, the G1 numerical method is used due to its ¨ simplicity. This simplicity leads from the fact that the Grunwald coefficients can be represented by means of the recursion formula Aa,j þ 1 ¼

ð29Þ

q¼1

thus, it can be deduced that the n order derivative should satisfy 0 1 ! N 1 X n ðnÞ n j f ðtÞ ¼ lim @ðDtÞ ð1Þ ð19Þ f ðtjDtÞA, j Dt-0

N¼

qmax X

j1a Aa,j , j

ð28Þ

where Aa,1 ¼ 1, which simplifies the treatment of the algorithm. Eq. (28) shows the path-dependence property or non-locality of the fractional derivative, because to evaluate the fractional order derivative of a function f in a particular point fn, the precedent values of the function until the lower derivative limit have to be taken into account.

ð30Þ

According to the fractional order differential equation resolution techniques exposed in [13], and making use of Eqs. (10)–(15), the analytical solution of Eq. (30) results in 1 ea1 : a1 Gða1 þ1Þ

ð31Þ

From the analytical solution (31) of the one-term model, some statements can be concluded. First, if the derivative order is zero, a1 ¼ 0, the model represents a perfect plastic behaviour (4), the parameter model being related to the yield strength by a1 ¼ s1 0 . Also, if the derivative order is equal to one, a1 ¼ 1, the model represents a perfect elastic behaviour (3), and in this case, the model parameter is related to the Young modulus by a1 ¼ E1 . Finally, if the derivative order is comprised between zero and one, 0 o a1 o1, the fractional model becomes the Hollomon one (1), satisfying a1 ¼ nH and a1 ¼ ðK  GðnH þ1ÞÞ1 . 3.2. Two-term particular case model If two terms are taken into account in Eq. (29), qmax ¼ 2, the model results in a1 Da1 sðeÞ þ a2 Da2 sðeÞ ¼ 1:

ð32Þ

Using the same resolution techniques as employed to solve the one-term case, the analytical solution of the two-term model is given by

sðeÞ ¼

1 X

ð1Þk ak1

kþ1 Gð 2 ðk þ1Þk 1 þ 1Þ k ¼ 0 a2

a

a

ea2 ðk þ 1Þka1 ,

ð33Þ

under assumption that the first derivative order is lower than the second one, and both are between zero and one, 0 r a1 o a2 r 1. From the first derivative of Eq. (33), it can be concluded that if the highest derivative order a2 is equal to one, and the lowest derivative order is between zero and one, a2 ¼ 1 and 0 o a1 o1, the initial slope of the stress–strain curve is equal to 1=a2 , i.e. a2 ¼ E1 , where E represents the Young modulus. On the contrary, if the higher derivative value is not equal to one, the initial slope is infinity. In this context, for the special case in which a1 ¼ 0 and a2 ¼ 1, Eq. (33) can be reduced to

sðeÞ ¼

1 ð1eða1 =a2 Þe Þ, a1

ð34Þ

3. Analysis of the proposed model

because the simplification of Eq. (33) for this particular case is just McLaurin’s development of Eq. (34), which corresponds to the solution of the differential equation

In this section, a new fractional elastoplastic model is presented. The general form of the proposed model is defined as a

a1 sðeÞ þa2 D1 sðeÞ ¼ 1:

ð35Þ

J. Mendiguren et al. / International Journal of Mechanical Sciences 65 (2012) 12–17

Then, Eq. (34) represents a perfect elastoplastic behaviour, without hardening, in which a1 ¼ s1 0 , representative of some metals. If the solution of the stress given by Eq. (33) has to be implemented in an incremental algorithm in which a finite strain step De is considered, the numerical resolution at the n þ 1 step provided by the G1 numerical method (26) can be used, which satisfies a1 Pn a 2 Pn 1 Aa ,j þ 1 sn þ 1j j ¼ 1 Aa1 ,j þ 1 sn þ 1j  ðDeÞa1 ðDeÞa2 j ¼ 1 2 sn þ 1 ¼ : a1 a2 a1 þ a2 ðDeÞ ðDeÞ ð36Þ

3.3. Three-term particular case model The three-term fractional model satisfies a1 Da1 sðeÞ þ a2 Da2 sðeÞ þ a3 Da3 sðeÞ ¼ 1:

ð37Þ

As for the previous cases, the resolution of this fractional differential equation follows the techniques exposed in [13], and results in   1  k  1  1 X a2 j X a1 ðj þkÞ! sðeÞ ¼ ð1Þk a3 j ¼ 0 a3 k!j! a3 k¼0 

1

Gðjða3 a2 Þ þkða3 a1 Þ þ a3 þ 1Þ

ejða3 a2 Þ þ kða3 a1 Þ þ a3 :

ð38Þ

@1

n n a1 X a2 X Aa1 ,j þ 1 sn þ 1j  Aa ,j þ 1 sn þ 1j a1 a2 ðDeÞ j ¼ 1 ðDeÞ j ¼ 1 2

1 n a3 X Aa ,j þ 1 sn þ 1j A:  ðDeÞa3 j ¼ 1 3

ð39Þ

4. Example of application to experimental data In order to validate the presented model, the three previously analysed particular cases of the general fractional model are fitted to experimental data and compared with the results provided by the Hollomon and Ramberg–Osgood classical models. Two different materials are used to validate the model. On the one hand, tensile test of a 6082O aluminum alloy, performed under a 0.001 s  1 strain ratio in room temperature, taken from [20]. On the other hand, tensile test of the AHSS TRIP 700 steel, performed under a 0.001 s  1 strain ratio in room temperature, taken from [21]. The curve fitting procedure was carried out by the Nelder and Mead [22] minimisation method, which is implemented in the s fminsearch function of Matlab . The objective function f ðxÞ is defined as f ðxÞ ¼

imax X

9sm ðei Þsi 9,

ð40Þ

i

For numerical applications, the numerical solution of Eq. (37) can be obtained using the G1 numerical method, the stress satisfying

sn þ 1 ¼

0

15

1 a1 a2 a3 þ þ ðDeÞa1 ðDeÞa2 ðDeÞa3

where the vector x represents the model parameters, i is the index of summation, imax is the total number of experimental data, 9  9 denotes the absolute value, sm ðei Þ and si are the stress values of the fitted model and the ith experimental stress for the ith experimental strain ei , respectively. Tables 1 and 2 show the obtained parameters for the five models of each material.

Table 1 Parameters for the fitted models for the 6082O aluminium alloy. Model

Parameters x

Hollomon, Eq. (1)

K (MPa) 235.77 E ðMPaÞ 70,000

nH 0.1812 nRO 5.7452

Ramberg–Osgood, Eq. (2) Fractional one-term, Eq. (30) Fractional two-term, Eq. (32) Fractional three-term, Eq. (37)

H ðMPaÞ 233.07

a1 ðMPa1 Þ

a1

4:5874  103

0.1812

a1 ðMPa1 Þ

a1

4:6411  103

0.1710

a2 ðMPa1 Þ 1=70,000

1

a1 ðMPa1 Þ

a1

a2 ðMPa1 Þ

a2

a3 ðMPa1 Þ

a3

4:4851  104

0

4:8967  103

0.0236

7:3084  106

1

a2

Table 2 Parameters for the fitted models for the AHSS TRIP 700 steel. Model

Parameters x

Hollomon, Eq. (1)

K ðMPaÞ 1253.9 E ðMPaÞ 203,000

nH 0.2202 nRO 4.8267

a1 ðMPa1 Þ

a1

8:7210  104

0.2202

a1 ðMPa1 Þ

a1

8:8637  104

0.2034

a2 ðMPa1 Þ 1=203,000

1

a1 ðMPa1 Þ

a1

a2 ðMPa1 Þ

a2

a3 ðMPa1 Þ

a3

2:4788  101

0

2:4716  101

0.0012

2:2962  106

1

Ramberg–Osgood, Eq. (2) Fractional one-term, Eq. (30) Fractional two-term, Eq. (32) Fractional three-term, Eq. (37)

H ðMPaÞ 1230.1

a2

16

J. Mendiguren et al. / International Journal of Mechanical Sciences 65 (2012) 12–17

Ramberg–Osgood model has as parameter the elastic modulus, so this parameter has fixed to 70 GPa for the aluminium alloy and 203 GPa for the TRIP steel. In the fractional two-term model the second derivative order has fixed to the unit value, a2 ¼ 1. Therefore, the constant associated to this derivative order has also fixed to the inverse of the initial slope, a2 ¼ 1=70,000 MPa1 in the case of the aluminium and a2 ¼ 1=203,000 MPa1 for the TRIP steel. From the results of Tables 1 and 2, the equivalence between the Hollomon and fractional one-term models can be verified, as was mentioned in the previous section. Figs. 1 and 2 represent the experimental stress–strain curve together with the provided by four fitted models for each material respectively. It can be taken into account that the one-term fractional model and the Hollomond model are equivalent, thus only the latter is represented. From Figs. 1 and 2, it should be remarked that the four model curves accurately represent the elastoplastic region, but differences can be found in the elastic region and in the elastic-to-elastoplastic transition, as well as in the large strain behaviour. In this context, Figs. 3 and 4 show a detail of the elastic and the elastic-to-elastoplastic transition, while Figs. 5 and 6 show the large strain behaviour.

Fig. 3. A detail of the stress–strain curve of the 6082O aluminium alloy: experimental and four fitted models, elastic zone.

Fig. 4. A detail of the stress–strain curve of the TRIP 700 steel: experimental and four fitted models, elastic zone.

Fig. 1. Stress–strain curve of the 6082O aluminium alloy: experimental and four fitted models.

Fig. 5. A detail of the stress–strain curve of the 6082O aluminium alloy: experimental and four fitted models, large strain.

Fig. 2. Stress–strain curve of the TRIP 700 steel: experimental and four fitted models.

In Fig. 3, it can be pointed out that each model gives different precisions, where the fractional three-term model is the most accurate. The fitting of the TRIP 700 steel confirms this statement.

J. Mendiguren et al. / International Journal of Mechanical Sciences 65 (2012) 12–17

17

 The proposed model has been validated by means of curve fitting to experimental measurements for the 6082O aluminium alloy and TRIP 700 AHSS steel. As a result, a new model capable of representing elastoplastic behaviour of metals better than current models has been proposed. However, the manipulation of this model implies a higher computational effort, thus the employment of this model is justified when accurate representation of the whole stress–strain curve is needed.

Acknowledgments The work presented in this paper has been carried out with the financial support of the Department of Education, Universities and Research of the Basque Government. Fig. 6. A detail of the stress–strain curve of the TRIP 700 steel: experimental and four fitted models, large strain.

Table 3 Minimized function value for each model after fitting. Model

6082O f ðxÞ (MPa)

TRIP 700 f ðxÞ (MPa)

Hollomon Ramberg–Osgood Fractional two-term Fractional three-term

264.7 238.3 359.2 121.6

161.8 156.6 203.1 75.1

On the other hand, in Figs. 5 and 6 it is shown how the fractional three-term model represents accurately the larger strain behaviour, therefore at larger strains the extrapolation using this model is expected to be better than using the other ones. These observations are confirmed taking into account the error indicator shown in Table 3. In Table 3 the value of the minimised objective function f ðxÞ is exposed for each model. This value can be taken as an error indicator to compare the model accuracy. Taking into account the results shown in Table 3, as previously mentioned, for this application example, the three-term fractional model gives the best representation of the experimental data, improving the Hollomond model in a 54% and the Ramberg– Osgood model in a 49% for the 6082O aluminium alloy and a 53% and 51% for the TRIP 700 steel.

5. Conclusions In this work, a new elastoplastic model has been developed, using fractional calculus techniques, in order to improve the elastoplastic behaviour modelling taking advantage of the nonlocality property of fractional operators

 Despite the correct accuracy of the Hollomon model to



represent a elastoplastic region, from the presented results it should be pointed out the limitations of Hollomon model to represent both elastic region and elastic-to-elastoplastic transition. The new three-term fractional model results more accurate than the classic Hollomon and Ramberg–Osgood models.

References [1] Dieter G. Mechanical metallurgy. McGraw Hill; 1986. [2] Eggertsen PA, Mattiasson K. On the identification of kinematic hardening material parameters for accurate springback predictions. Int J Mater Form 2011;4(2):103–20. [3] Lemoine X, Durrenberger L, Zhu H, Kergen R. Mixed hardening models: parameters identification on AHSS steels. In: IDDRG2011; 2011. [4] Stepanskiy L. Probabilistic model of strain hardening of Fe–Mn-based austenitic steels. Scr Mater 2009;61(10):947–50. [5] Li X, Yang Y, Wang Y, Bao J, Lo S. Effect of the material-hardening mode on the springback simulation accuracy of V-free bending. J Mater Process Technol 2002;123(2):209–11. [6] Mendelson A. Plasticity: theory and applicationMacmillan; 1968. [7] Kachanov L. Foundations of the theory of plasticity. Amsterdam: NorthHolland; 1971. [8] Corte´s F, Elejabarrieta MJ. Finite element formulations for transient dynamic analysis in structural systems with viscoelastic treatments containing fractional derivative models. Int J Numer Methods Eng 2007;69(10):2173–95. [9] Corte´s F, Elejabarrieta MJ. Homogenised finite element for transient dynamic analysis of unconstrained layer damping beams involving fractional derivative models. Comput Mech 2007;40(2):313–24. [10] Schmidt A, Gaul L. Fe implementation of viscoelastic constitutive stress– strain relations involving fractional time derivatives. Constitutive models for rubber, vol. 2; 2001. p. 79–92. [11] Bagley R, Torvik P. On the fractional calculus model of viscoelastic behavior. J Rheol 1986;30(1):133–55. [12] Miller K, Ross B. Introduction to the fractional calculus and fractional differential equations. New York: John Wiley & Sons; 1993. [13] Podlubny I. Fractional differential equations. New York: Academic Press; 1999. [14] Vinagre BM, Monje CA. Introduction to the fractional control. Rev Iberam Autom Inf Ind 2006;3(3):5–23. [15] Oldham K, Spanier J. The fractional calculus. New York: Academic Press; 1974. [16] Schmidt A, Gaul L. On the numerical evaluation of fractional derivatives in multi-degree-of-freedom systems. Signal Process 2006;86(10):2592–601. [17] Diethelm K, Ford NJ, Freed AD, Luchko Y. Algorithms for the fractional calculus: a selection of numerical methodsComput Methods Appl Mech Eng 2005;194(6–8):743–73. [18] Podlubny I, Chechkin A, Skovranek T, Chen Y, Vinagre Jara B. Matrix approach to discrete fractional calculus II: partial fractional differential equationsJ Computat Phys 2009;228(8):3137–53. [19] Podlubny I, Skovranek T, Vinagre Jara B. Matrix approach to discretization of fractional derivatives and to solution of fractional differential equations and their systems. In: 2009 IEEE 14th international conference on emerging technologies & Factory Automation ETFA 2009. p. 6. [20] Aginagalde A, Go´mez X, Galdos L, Garcı´a C. Heat treatment selection and forming strategies for 6082 aluminum alloy. J Eng Mater Technol 2009;131: 044501. [21] Perez R, Benito J, Prado J. Study of the inelastic response of TRIP steels after plastic deformation. ISIJ Int 2005;45(12):1925–33. [22] Nelder JA, Mead R. A simplex method for function minimization. Comput J 1965;7(4):308–13.