Notes on computational aspects of the fractional-order viscoelastic model

This paper deals with the computational aspect of the investigation of the relaxation properties of viscoelastic materials. The constitutive fractional Zener model is considered under continuous deformation with a jump at the origin. The analytical solution of this equation is obtained by the Laplace transform method. It is derived in a closed form in the terms of the Mittag-Leffler function. The method of numerical evaluation of the Mittag-Leffler function for arbitrary negative arguments which corresponds to physically meaningful interpretation is demonstrated. A numerical example is given to illustrate the effectiveness of this result.

The theory of viscoelasticity is crucial in describing materials, such as rubber or polymers, which exhibit timedependent stress-strain behaviour.In many industry applications, like tyres, shock absorbers or other rubber components, a reliable constitutive equation must be established in order to describe all the relevant effects in the material response under small harmonic vibrations.The classical model of finite viscoelasticity is based on a finite series of Maxwell elements in parallel, which leads to a relaxation function given by a sum of decreasing exponentials (Prony series) (see [12]).Experimental investigations into the time-dependent relaxation behaviour of such materials are characterised by a very fast decrease of the stress at the beginning of the relaxation process and an extremely slow decay for large times.It can even take a long time before a constant asymptotic value is reached.Therefore, a lot of Maxwell elements are required to describe relaxation behaviour with sufficient accuracy.This ends in a difficult optimisation problem in order to identify a large number of material parameters.On the other hand, over the years, the concept of fractional derivatives has been introduced to the theory of viscoelasticity [13][14][15][16].Among these models, the fractional Zener model was found to be very effective to predict the dynamic nature of rubber-like materials with only a small number of material parameters [17][18][19][20].The solution of the corresponding constitutive equation leads to a relaxation function of the Mittag-Leffler type.It is defined by the power series with negative arguments.This function represents all essential properties of relaxation process under the influence of an arbitrary and continuous signal with a jump at the origin.However, from the numerical point of view, it is inconvenient to calculate the numerical values of the Mittag-Leffler function for large negative arguments.Our motivation is to discuss and solve this computational difficulty against a background of physical and application interest.
The paper is organized in the following way.In the first part (Sect.2), we recall the main properties of the Mittag-Leffler function and introduce the linear operators of fractional differentiation.In the second part (Sect.3), we consider differential equation of fractional order, which model the constitutive relation for viscoelastic materials.For this equation, we derive its analytical solution by applying the technique of Laplace transforms.It is given in the terms of the Mittag-Leffler function.We present the method of calculating numerical values of Mittag-Leffler function for arbitrary negative arguments, which is based on its integral representation.Finally, a numerical simulation is given to illustrate the effectiveness of the proposed approach.Section 4 contains some concluding remarks.

Preliminaries and notations
We start our considerations by introducing several special functions, fractional operators and their properties that will be useful in the second part of this work.Some definitions in this section are very common, and we write them just for reference.Others are far from being common but are still quite often used in the literature [1,21,22].

Laplace transform
Definition 2.1 Let f : R + → R be a locally integrable function on the interval [0, ∞) and s be a complex number.If the improper integral exists, then we call it the Laplace transform of f at s.

Definition 2.2
Let f : R + → R be a real function and suppose that F(s) = L[ f (t)](s) is well defined for some s = γ + iθ .Then, we have and it is called the inverse Laplace transform of F.
Property 2.1 Let f and g be any functions and c 1 and c 2 be any constants for which the needed operations are defined.Then, we can write the following properties of the Laplace transform Here the convolution is defined by The Heaviside step function may be defined equivalently as the integral of the Dirac delta function, that is 4 The Gamma function is defined by the integral It can be proved by integrating by parts that the Gamma function satisfies the following functional equation for arbitrary x > 0.
123 Definition 2.5 Let α > 0 and β > 0 be two real but arbitrary numbers.The two-parameter Mittag-Leffler function E α,β is defined by whenever the series converges. Since is the well known exponential function, the Mittag-Leffler function is also known as the fractional exponential function. Theorem for all complex numbers s satisfying the condition |a/s α | < 1.

Fractional derivatives
In the following subsection, we shall show the most popular definitions of the fractional-order derivative and describe their properties.We also present some similarities, differences and relations between them (see [1,21,22]).Let us start with the most popular definition of fractional derivative being a result of Liouville's and Riemann's work.It is based on a generalization of the Cauchy formula for n-fold integration.To be more precise, let f be a sufficiently regular function on the interval for t ∈ [a, b] and n ∈ N. Replacing the factorial operator with the Gamma function, being its natural generalization, Liouville and Riemann have created an integral operator in the following form is finite, then we call it the Riemann-Liouville fractional derivative with order α of function f .
In 1967, an Italian mathematician Michele Caputo [23] gave a different definition of fractional derivative, which has turned to be very useful in industrial problems.

Definition 2.8
Let α be an arbitrary real number in the interval (0, 1) and which is called the Caputo fractional derivative with order α of function f .
Let us investigate some basic properties of these differential operators of fractional order.First of all, we notice that Definitions 2.7 and 2.8 are consistent with the classical integer order derivative, i.e. (2.17) Property 2.2 Let f, g be arbitrary functions and c 1 , c 2 , α, a be real constants for which the needed operations are defined.Then the linearity of the Riemann-Liouville and Caputo fractional derivatives can be expressed by and (2.20) Proposition 2.1 Let f be arbitrary functions and α be a real constant for which the needed operations are defined.
Then the Laplace transform of the Riemann-Liouville and Caputo fractional derivatives satisfy Practical applicability of the Laplace transform of the Riemann-Liouville fractional derivative is limited by the absence of the physical interpretation of the limit values of fractional derivative at t = 0. Since there is no such problem in the case of Caputo fractional derivative, it can be useful for solving applied problems leading to linear fractional differential equations.

Fractional viscoelastic constitutive model
In this section, we present the application of fractional derivatives in the constitutive model of viscoelasticity.To be more precise, we formulate the Caputo type fractional differential constitutive equation and solve it, applying the Laplace transformation technique for an arbitrary signal with a jump at t = 0. Furthermore, we discuss the inconvenience of numerical calculations of the relaxation function derived in terms of the Mittag-Leffler function.

Generalized Zener model of viscoelastic materials
The fractional Zener model is the generalization of the well-known standard linear solid model consisting of an ideal elastic spring and Maxwell element connected in parallel.The qualitative behaviour of the conventional model is improved by replacing the integer order time derivatives of stress and with fractional-order derivatives.
The constitutive equation for the fractional Zener model is where σ is the stress, ε is the, A 0 , B 0 , B 1 and α ∈ (0, 1) are material constants,-order and the α-th-order fractional derivative is defined by the Caputo operator (2.16).This model has been proved to be efficient in describing the stress-strain relation of real viscoelastic materials.

Analytical solution of the fractional-order differential constitutive equation
Let us consider the viscoelastic material subjected to an imposed arbitrary continuous deformation ε(t) with a jump at t = 0. We assume that the stress is zero before the is loaded, i.e. the material is in equilibrium.According to the Boltzmann superposition principle in viscoelasticity, we can express the signal as the sum of constant ε 1 and continuous ε 2 functions (see Fig. 1).Then, we have where H denotes the Heaviside step function.Linearity of the constitutive equation (3.1) and linearity of the Caputo fractional operator (Property 2.2) allows us to measure the corresponding stress response which can be expressed as where σ 1 (t) and σ 2 (t) are the stress histories corresponding to ε 1 (t) and ε 2 (t) seperately.Let us start with the equation Inserting ε 1 (t) = ε 0 H (t), applying the Laplace transform, using Property 2.1 and doing simple transformations, we get Coming back to the time domain with the use of Proposition 2.1, we obtain Next we proceed with the second term ε 2 of the signal; thus, we have Keeping in mind that ε 2 vanishes at t = 0, we get now Applying the inverse Laplace transform and using Borel's theorem for the Laplace transform of a convolution leads to Since the only part of ε that is responsible for changes is ε 2 , we have ε (t) = ε 2 (t) for t > 0.Then, we get the final stress-strain relationship on the basis of (3.6) and (3.9) with the kernel k given in the form Using Eq. (2.10), we can rewrite it in a more compact way: Then, we have Fig. 1 The input signal ε(t) and its decomposition into constant part ε 1 (t) and continuous part ε 2 (t) where μ 0 = B 0 /A 0 and G is the so-called viscoelastic relaxation function defined as with parameter μ 1 = B 1 − μ 0 .The asymptotic behaviour of the Mittag-Leffler function for arbitrary negative argument with 0 < α < 1 can be described in the following way [12,21] for x → 0, and for x → +∞.Thus, it follows that the relaxation function G defined by (3.14) is able to represent physical properties of relaxation behaviour (i.e. a rapid decrease for small times and a very slow decay for large times) on the basis of only four material parameters A 0 , B 0 , B 1 and α, which define the so-called continuous relaxation spectrum.
The plots of G for different values of parameter α are illustrated in Fig. 2. The function G has a finite value and a vertical tangent at t = 0. Its decrease becomes extremely slow for large values of t if α is small.However, we have noticed computational difficulty in the numerical evaluation of the relaxation function (3.14) for large negative arguments.This problem arises as a result of summing series of terms with alternating signs using floating-point arithmetic.Therefore, we present another method of calculation for the Mittag-Leffler function, which is based on the continuous relaxation spectrum.

Integral representation of the Mittag-Leffler function
The concept of a general relation between the complex modulus of a linear viscoelastic system and its spectrum is well-known in the literature.We follow [24] and relate the relaxation function G to the continuous relaxation In practical applications, where the function G must be calculated numerically, it is more convenient to introduce the so-called cumulative relaxation spectrum: which is monotonically an increasing function such that h c (0) = 0 and h c (∞) < ∞.Then, integrating by parts (3.17), we obtain an alternative representation of the relaxation function: However, to take advantage of this formula, we have to first answer the question how to determine h c corresponding to the fractional Zener model and its analytical solution (3.13), which specifies the one-dimensional stress-strain relation.
In order to describe the relaxation properties of viscoelastic material, the dynamical tests are performed in which the harmonic oscillations are applied to a material specimen, and the resulting stress is measured until a steady-state value is reached.It is common practice in engineering to use complex variables to present the harmonic response.To this end, we assume the complex harmonic function where ω is the angular frequency, Δf is the amplitude and H is the Heaviside function.Applying the Caputo integral operator (2.16) and calculating the fractional derivative of the function f in the stationary case, i.e. for large values of t (for details see [1,4,22]), we have Without loss of generality, we assume ε 0 = 0 and consider the complex and stress as Then, by inserting a complex-valued deformation process ε into Eq.(3.13), we have with the quantity which defines the so-called complex dynamic modulus for a stationary stress response.Taking into account the representation (3.17) of the relaxation function G and carrying out the integration with respect to the variable τ yields The formula (3.25) shows that the complex modulus G * and relaxation spectrum h are related to each other.If a relaxation spectrum is given, the complex dynamic modulus is determined.On the other hand, we notice that the function (G * (iω) − μ 0 )/iω is the Stieltjes transformation of the relaxation spectrum h, i.e.

.26)
Assuming that h is an absolutely integrable function on [0, R], for every finite R > 0, one may calculate the inverse of the Stieltjes transformation using the formula where v is a positive value for which the function h is defined.It is convenient to express the complex arguments −v ± iξ in the form of the complex exponential functions (3.28) Then, calculating the limit ξ → 0 + , we obtain the general formula which can be used to compute the continuous relaxation spectrum if the complex dynamic modulus is known.Moreover, bearing in mind (3.21) and inserting the harmonic deformation ε in combination with the harmonic stress σ into (3.1),we transform the fractional Zener model into the frequency domain.Thus, we obtain again the expression (3.23), where the complex dynamic modulus has a form (3.30) The formula (3.29) in combination with (3.30) can be applied to calculate the relaxation spectrum h of the function G given by the Mittag-Leffler function (3.14).We obtain which finally leads to (3.32) For further calculations, it is useful to consider the cumulative relaxation spectrum h c given by (3.18).Then, analytical integration of function (3.32) yields the following formula This result provides an interesting formula for the efficient numerical calculation of the Mittag-Leffler function E α,1 for large negative arguments.Indeed, in combination with (3.19) and (3.33), we can formulate the following corollary.
From a practical and applicable point of view, we distinguish two cases for the numerical computation of E α,1 (x) with the prescribed accuracy ε > 0: (i) x < 0, |x| ≤ q < 1, (ii) x < 0, |x| > q, q < 1, where q is a fixed number.In case (i), we apply the series (3.35), whereas in case (ii), we use the integral representation (3.40).Since the number of terms in the sum of formula (3.35) in Theorem 3.1 depends on x, it is recommended to take a number q not very close to 1. Indeed, on the basis of (3.36), we can notice that the smaller value of q is taken, the smaller number of terms N is used in formula (3.35).
Using the above approach, we can quite easily evaluate the value of the relaxation function G given by the relation (3.47) Defining T = (q/A 0 ) 1/α for given q < 1, the Mittag-Leffler function is numerically calculated by formula (3.35) for 0 < t ≤ T and (3.40) for t > T , respectively.Our numerical results for arbitrary values of t, with parameters q = 0.8, A 0 = 0.5, μ 1 = 8, different values of α and accuracy ε = 2.2204 × 10 −16 (= machine epsilon), are [a, b] ⊂ R and I denote the integral operator over [a, b].The Cauchy integral formula is