© The Authors, published by EDP Sciences, 2018

1 Introduction

Fractional calculus is a branch of mathematical analysis that generalizes the derivative and integral of a function to non-integer order [1]. Application of fractional calculus in classical and modern physics greatly contributed to the analysis and our understanding of physico-chemical and bio-physical complex dynamical systems, since it provides excellent instruments for the description of memory and properties of various materials and processes. During the last two decades fractional calculus has been increasingly applied to mathematical modelling in physics [2,3], engineering [4,5], and especially to rheology [6,7], where fractional calculus constitute a valuable mathematical tool to handle viscoelastic aspects of systems and materials mechanics. Models involving fractional derivatives and operators have been found to better describe some real phenomena than integer-order differential equations [1-5], whence there are many new exciting areas of fractional models and fractional calculus applications, such as the automatic control, the modeling of biological, medical and environmental. A historical review of applications can be found in [8]. The survey on fractional models from biology and biomedicine is presented in [9], see also other papers cited therein. In recent decades fractional derivatives were found quite flexible, especially in the description of viscoelastic polymer materials [10].

Viscoelastic materials present a behaviour that implies dissipation and storage of mechanical energy. Research studies conducted during the past few decades proved that these models are also an important tool for studying the behaviour of biological materials [11]: wood, fruits, vegetables, animals tissues, etc. Viscoelasticity of the materials manifests itself in different ways, such as gradual deformation of a sample of the material under constant stress (creep behaviour), and stress relaxation in the sample when it is subjected to a constant strain. In general, viscoelasticity is a phenomenon associated with time variations in a material’s response. In an attempt to describe some of the above effects mathematically several constitutive laws have been proposed which describe the stress–strain relations in terms of quantities like creep compliance, relaxation modulus, the storage and loss moduli and dynamic viscosity. Some of these constitutive laws have been developed with the aid of mechanical models consisting of combinations of springs and viscous dashpots. The classical Maxwell model, which is a viscoelastic body that stores energy like a linearized elastic spring and dissipates energy like a classical fluid dashpot is, perhaps, the most representative example of such models.

For over five decades classical exponential behaviour models have been widely applied to describe the viscoelastic properties of biological materials. Maxwell, Kel- vin-Voight and Zener models are used to mathematical modelling of stress relaxation and creep processes [11-13]. For these models the relationship between the stress and deformation of the material is approximated though an ordinary differential or integral equations.

However, relaxation or creep processes deviating from the exponential Debye decay behaviour are often encountered in the dynamics of biological complex materials [6, 13]. For such materials a stretched exponential decay KWW model (Kohlrausch-Williams-Watts) [12], hyperbolic type decay Peleg model [13] or power type behaviour models [14,15] are used to approximate experimentally obtained relaxation modulus or creep compliance data. Recently, the experimental results obtained by the authors show that the behaviour of some viscoelastic biological materials shows good agreements with that of the fractional models [9,10,16]. To this end, fractional rheological models, originally pioneered by Nutting [17] and Scott Blair [18], have proven to be a concise and elegant framework for predicting the response of complex fluids such as liquid foods using a small number of parameters [16]. By replacing the springs and dashpots of the classical viscoelastic models by the Scott Blair elements, several fractional models, including the fractional Maxwell, fractional Voigt and fractional Kelvin models, have been proposed [6, 19].

In this paper fractional Maxwell model (FMM) model is considered, which relates the stress to the strain in the material by means of using differential fractional equation [3,6,7]. The FMM admit the closed form of analytical solution in terms of the known Mittag-Leffler function [20]. While fitting data to exponential sum models, like classic Maxwell model, is a very old problem, which has been studied for a long time, until now there are only a few papers concerned mainly to finding fractional Maxwell model.

We often determine the parameters in a model by obtaining the „best-possible" fit to experimental data. Common choice of the model quality measure is the mean square approximation error, leading to a least-squares identification problem. Unfortunately, the original nonlinear least-squares task for FMM is very difficult by means of the infinite series representation of the Mit-tag-Leffler function [20]. Thus, in [7, 21] an approach is considered, which consists on the determination of the FMM parameters by two approximate Scott Blair models identification in two separate time intervals, for small and large times respectively.

The aim of the paper is to develop a complete procedure for FMM identification using relaxation modulus data from ramp stress relaxation test.

2 Material and Method

2.1 Material

We consider a linear viscoelastic material subjected to small deformations for which the uniaxial, nonaging and izotropic stress-strain equation can be represented by a Boltzmann superposition integral [22]:

(1)

where σ(t) and ε(t) denotes the stress and stain, respectively, and G(t) is the linear time-dependent relaxation modulus. The modulus G(t) is the stress, which is induced in the viscoelastic material described by equation (1) when the unit step strain ε(t) is imposed. By assumption, the exact mathematical description of the relaxation modulus G(t) is completely unknown, but the value of G(t) can be measured with a certain accuracy for any given value of the time t.

2.2 Fractional Maxwell model

Fractional Scott Blair model [7] is described by the fractional differential equation

(2)

where E and τ are the elastic modulus and relaxation time, respectively, and α is non-integer positive order of fractional derivative of the strain ε(t). Here, means the derivative operator in the sense of the Ca- puto’s fractional derivative of a function f(x) of noninteger order y with respect to variable t and with starting point at t = 0, which is defined by [3]

where Γ(n) is Euler’s gamma function, n – 1 < γ > n.

The Scott Blair model (2) is intermediate model between ideal spring (see Figure 1a) Hooke’s model

and Newton’s model

of ideal fluids, represented by means of an ideal dashpot (see Figure 1b). To illustrate the structure of fractional models a fractional element, in addition to the standard purely elastic and purely viscous element, must be introduced - see Figure 1c.

Assuming unit-step strain ε(t) the uniaxial stress response of elementary fractional element (2), i.e. the time-dependent relaxation modulus G(t), is given by [6,7]

(3)

Thus, the elementary fractional element (3) is uniquely described by three parameters (E, τ, σ), as shown in Figure 1c. This approach makes it possible to include a whole range of dissipation mechanisms in a single (three parameter) rheological element.

Classic viscoelastic Maxwell model is the arrangement of ideal spring in series with a dashpot (see Figure 2a) described by the first order differential equation

(4)

which for unit-step strain ε(t) has exponential type response

(5)

with the relaxation time τ = η/E.

By replacing the spring and dashpot of the classic Maxwell model by two elementary Scott Blair elements (E₁, τ₁, α) and (E₂, τ₂, β) - see Figure 2b - we obtain fractional Maxwell model described by the fractional differential equation [6]. Without loss of generality we assume α ≥ β. Since, in view of (2), for the same stress σ(t) the strains for these elements are described by

the sum strain ε(t), which is induced in the series connection of two elementary fractional elements is described by the fractional equation

which using the composition rule for fractional differentiations

can be rewritten in compact form as

(6)

where the parameters of FMM (6) are functions of the parameters (E₁, τ₁, α) and (E₂, τ₁, β) of the model components given by

For details of the fractional model construction see [6,7]. For the unit-step strain the solution σ(t) = G(t) of FMM (6) is known for an arbitrary 1 > σ > β > 0 and given by the formula [6,7] :

(7)

where E_{κ, μ} (x) is the generalized two-parameter Mittag- Leffler function defined by series representation, convergent in the whole z-complex plane [3,20]

(8)

The above result is obtained in [6] by applying Fourier and Mellin transforms, while in [23] the same result is obtained by an approach involving the Laplace-Mellin transform technique.

Note, that the fractional Maxwell model (7) is uniquely defined by four parameters (E, τ, α, β), while the classic Maxwell model (4) is defined by means of only two parameters (E, η), or equivalently (E, τ). A deep insight into the complex properties of FMM gives the infinite hierarchical structures of elementary fractional model (2) presented in Figure 3 composed of ideal spring and dashpot elements, which compose the classic Maxwell model.

The course of the relaxation modulus G(t) (7) for fixed β parameter and a few values of the non-integer derivative order α are shown in Figure 4, for fixed α and a few values of β the respective modulus G(t) are plotted in Figure 5. Note, that the curves of FMM relaxation modulus shown in Figures 4 and 5 has the characteristic shape of the relaxation modulus of viscoelastic materials obtained in experiment. The following rule holds: the greater the parameter α is, the shorter the relaxation times are for β being fixed. Similar rule holds for model order β, while α is fixed.

Fig. 1. From left to right: the Hooke (a), Newton (b) and elementary fractional element (c).

Fig. 2. Classic Maxwell model (a) followed by the fractional Maxwell model (b).

Fig. 3. Sequential ladder realization on the fractional Scott Blair element [6; Fig. 2].

Fig. 4. The FMM relaxation modulus for fixed β = 0.1.

Fig. 5. The FMM relaxation modulus for fixed α = 1.

2.3 Approximate fractional models

In the interpretation of the asymptotic behaviour of the FMM a fundamental role plays the asymptotic expansions of the Mittag-Leffler function. Based on (8), the relaxation modulus G(t) (7) can be rewritten as

(9)

and since for small t, in particular for 1, the infinite sum of the right hand side of (9) can be truncated, the asymptotic approximation results

(10)

here ≅ means 'approximately equal'. From the following asymptotic property [24]

(11)

which holds for t → ∞, putting φ = α – β and μ = 1 – β in (11) based on (7) yields

whence for large times, especially for we obtain the next asymptotic approximation of FMM in the form

(12)

The two approximate rules (10), (12) are given in [6] for short times such that t ≪ τ and for large times t ≫ τ, respectively, based on power series representations of H-function and fact that the Mitag-Leffler function is expressible through H-function. The same asymptotic approximate models are considered in [7] and [21], but in the last paper stringent applicability conditions, derived above, are imposed. Thus, both for short as well as for large times, relaxation modulus (7) decreases almost according to the time-power elementary fractional element. The approximate models (10), (12) and the exact FMM are summarized in Figure 6, where logarithmic scale is used both for relaxation modulus and time scales. Here, the asymptotic models (10), (12), whose graph is a straight line are linear functions, which slope coefficients uniquely determine the fractional derivatives orders - β for small times and α for large times, respectively.

On the basis of this observation in [7] the FMM parameters (E, τ, α, β) of polymers PMMA and PTFE has been determined based on approximate models (10) and (12) in the time intervals defined by t ≪ τ and t ≫ τ, respectively. Next, in [21] the PMMA and PTFE fractional models have been corrected by applying more restrictive conditions and . In result better fit to experimental data is obtained.

In this paper a complete two-interval identification scheme based on approximate models (10) and (12) and using the linear least-squares is presented and applied to FMM of biological material identification. The scheme is inspired by [7,21] papers. However, precise applicability conditions are formulated here for the first time.

Fig. 6. FMM and approximate Scott Blair models for parameters: α = 0.8, β = 0.1, E = 1 [Pa], τ = 0.9 [s].

3 Results and discussion

3.1 Identification

A classical manner of studying viscoelasticity is by two- phase stress relaxation test, where the strain increases during the loading time interval until a predetermined strain ε₀ is reached at given ramp-time, after which that strain ε₀ is maintained constant at that value [11, 12, 25].

Suppose, a certain stress relaxation test performed on the specimen of the material under investigation resulted in a set of measurements of the relaxation modulus at the sampling instants t_t > 0, i = 1,..., N, where z(t_t) is measurement noise. For computational methods of relaxation modulus determination see, for example [12,25].

In general, identification consists of selecting within the given class of models such a model, which ensures the best fit to the measurement results. Fitting data to the original FMM (7) is a very difficult problem of nonlinear optimization, numerically difficult and often ill- conditioned.

Here, the linear least-squares identification routine will be used to estimate FMM parameters based on the logarithmic transformation of the experimental data and equations (10) and (12) which, respectively, yields

where

(13)

14

are introduced for brevity. Denoting χ(t) = log G(t) and introducing the new independent variable ϑ = log t, equivalent linear models are obtained (the sign = is neglected for simplicity)

15

(16)

Now, in the set of experimental data , or equivalently in the set {(χ_t, ϑi)}_i^N=₁, where

(17)

are log-transformed 'measurement' data, two separable subsets and must be determined such that n₁ < n₂, probably n₁ ≪ n₂. The subsets N = {1,…,n₁} and N₂ = {n₂, …,N} of respective indices are chosen during the recurrent identification scheme. It is clear that the union of the sets N U N₂ does not create the set N = {1, …,N}. Now, classical linear least squares method can be applied to find optimal approximate models.

As a measure of the model (15) accuracy the mean sum of squares is taken

(18)

the respective identification index for the second set of experimental data and log-linearized model (16) is

(19)

Therefore, the least-squares identification of the log- linearized models consists of determining the model parameters minimizing the indices (18), (19) by solving the following standard optimization problems

(20)

(21)

Based on the commonly known results concerning the linear least-squares problem solution, it can be shown that the solutions to (20), (21) tasks there exist and are unique whenever the sampling instants are such that 0 < t₁ < t₂ < … < t_N. The model parameters optimal in the sense of (21) are given by the known formulas:

(22)

(23)

The formulas for optimal and are analogous. According to definitions (13) and (14) the positivity constraint can be neglected for k₂,k₂ in (20), (21) optimization tasks. However, both α as well as β must be positive. It can be shown, based on the Czebyszev equality, that the optimal parameters and are positive, if and only if, the following conditions are satisfied

for tasks (20) and (21), respectively.

Note, that the models (15) and (16) include two unknown parameters (k₁,β) and (k₂,α). In view of (14) and (13) the pairs of parameters are not independent. Therefore, from (14) and (13) we have

(24)

which, in view of (14), next yields

(25)

3.2 Algorithm

Taking into account the above, the calculation of the approximate values of FMM parameters involves the following steps.

• 1.

  Perform the stress relaxation test, record and store the relaxation modulus measurements corresponding to the chosen sampling instants t_t > 0.

• 2.

  Determine the set of log-transform measurement data (17).

• 3.

  Based on the course of log-log plot of the measurement data, select in the set of measurements two separable subsets choosing n₁ and n₂, n₁ < n₂.

• 4.

  Compute the estimates and according to formulas (22), (23) and, next, compute the estimates of the relaxation time and elastic modulus based on (24) and (25).

• 5.

  In order to ascertain if the models with parameters are a satisfactory approximation of measurement data in time intervals compute the optimal identification indices and examine if and for S, a preselected small positive error. If not, change n_i and/or n₂ and go to step 4. Otherwise, go to step 6.

• 6.

  In order to ascertain if the Scott Blair models (10), (12) with parameters and are a satisfactory approximation of the original FMM examine if (26) (27)

  If both the above conditions are satisfied, stop the procedure taking as the FMM parameters. If not, change n₁ and/or n₂ and execute steps 4 and 5.

The stopping rules from Step 5 guarantee the good quality of the log-linearized models in the chosen time intervals and correspond with those commonly used in the optimal identification techniques. The conditions (26), (27) are imposed to guarantee the applicability of Scott Blair models (10), (12) to approximate the original FMM. Both the pairs of conditions must be satisfied simultaneously in order to guarantee good quality of the resulted model. Note also that, in particular, from (27) the obvious necessary applicability condition of the scheme results

i.e. the experiment time must be significantly larger than the relaxation time of the FMM.

While and are determined independently in appropriate time intervals, the relaxation time and depend on the identification results in both time intervals, simultaneously.

Note, that (26) and (27) are a posteriori conditions, since the applicability of the identification procedure cannot be checked, earlier than after the experiment is performed.

The next example shows, how the identification scheme can be used for identification of FMM of real material.

3.3 Example - FMM of carrot root

A cylindrical sample of 10.9 mm diameter and 8 mm height was obtained from the root of carrot Nantejska variety [14]. During the two-phase stress relaxation test performed by Bohdziewicz [14], in the first loading phase of the time period (0,43.7) seconds the strain rate 0.000033 [m · s^-i] is used and the force induced in the specimen was recorded with the mean sampling period Δ t_mean = 0.11007 [s]. Next, during the second phase of the time period (43.7,1187) seconds, at constant strain the corresponding time-varying force was recorded in N = 341 measurement points with sampling period Δt_mean = 3.34s (for details see [13,14]). Next, the relaxation modulus measurements were computed for the time interval (0,1143.3) seconds using the fast trapezoidal method of approximate relaxation modulus identification presented in [25].

Based on the log-log plot of the experiment results, the time intervals determined by t_ni = 0.5 [s] and t_n2 = 300 [s] have been chosen and the estimates of Scott Blair models parameters, and next the estimates of the relaxation time and elastic modulus have been computed. Freely available [26] Matlab code for evaluating the Mittag- Leffler function has been used. The inaccurate fit of the resulting FMM to the experimental data is illustrated in Figure 7. Next, the recurrent scheme of the algorithm have been realized and successive approximations of FMM parameters have been determined.

The results are summarized in Table 1, where both the model parameters, as well as the factors ε_{n1, n2} (26) and Ω_{n1, n2, mean}, defined as the mean value of the left hand side Ω_{n1, n2} of (27), are given for a few FMM approximations. Repeated composition of the scheme operations has led to parameters listed in the last row of Table 1. For these model parameters the accurate fit to experiment data is achieved, as shown in the Figure 8. The presented results convincingly prove, that the FMM give satisfactory approximation of the relaxation modulus of the material.

The last column of the table is a value of the empirical mean square model error

(28)

where a = (E, τ, α, β) is the vector of FMM parameters. The model errors are not big (see Table 1), but for a few first FMM approximations the accuracy is insufficient, especially in long time region.

Next, the generalized discrete Maxwell model, which presents a relaxation of exponential type given by a finite Dirichlet-Prony series [11,12]

(29)

where E_j, V_j and E_∞ represent the elastic modulus, relaxation frequencies and equilibrium modulus (long-term modulus), respectively, has been applied to approximate the experiment data. Least-squares approximations are used, and two, three and four-parameter optimal models were determined. The parameters of the classic optimal Maxwell models and the relaxation times τ_j = 1/V_j are given in Table 2. The differences between the three Maxwell models are characterized by the ERR errors (28) summarized in Table 2, here the vector a = (E_j, V_j, E_n). Figure 9 presents experiment data and the optimal Maxwell models.

It is easy to observe, that the Maxwell model is inappropriate for description of this relaxation process. Thus, the classic Maxwell models do not fully characterize the true viscoelastic behaviour of the biological material. A better fit to experimental data can be obtained, if the FMM is used. Note, that also the range of the relaxation time estimates is lower for FMM than for classic Maxwell models, even if the conditions (26), (27) of the scheme applicability are not sharply fulfilled - compare Tables 1 and 2.

Fig. 7. The inaccurate fit of FMM to experiment data for parameters: α = 0.6092, β = 0.0159, E = 6.9425 [MPa] and and τ = 384.7876 [s]

Fig. 8. The time course of FMM and the experimental data, parameters: α = 0.6892, β = 0.0139, E = 7.0425 [MPa] and τ = 457.0883[s].

Fig. 9. The optimal Maxwell models fit to the experimental data.

4 Conclusions

In this work, a complete simple procedure is proposed for the fractional Maxwell model identification based on the discrete-time stress relaxation experimental data. The asymptotic approximate log-linearized Scott Blair models, composed of only two parameters, are used in separate time intervals. The linear least squares are applied to find the best Scott Blair models in two separable appropriately chosen time intervals. An applicability conditions of the approximate identification scheme are derived. It is shown that FMM can be used to describe the viscoelastic mechanical properties of biological materials. The effectiveness is compared with the classical two, three and four parameter Maxwell model approximations. While the idea of using the approximate Scott Blair models to find the FMM parameters is known and already appeared in [7,21], the complete identification scheme is presented in this paper for the first time.

The authoress wish to thank Jerzy Bohdziewicz of the Faculty of Life Sciences and Technology, Wroclaw University of Environmental and Life Sciences, Wroclaw, Poland, for providing measurement results of relaxation tests of carrot roots and helpful suggestions concerning the rheological experiments.

References

All Tables

All Figures

	Fig. 1. From left to right: the Hooke (a), Newton (b) and elementary fractional element (c).
In the text

	Fig. 2. Classic Maxwell model (a) followed by the fractional Maxwell model (b).
In the text

	Fig. 3. Sequential ladder realization on the fractional Scott Blair element [6; Fig. 2].
In the text

	Fig. 4. The FMM relaxation modulus for fixed β = 0.1.
In the text

	Fig. 5. The FMM relaxation modulus for fixed α = 1.
In the text

	Fig. 6. FMM and approximate Scott Blair models for parameters: α = 0.8, β = 0.1, E = 1 [Pa], τ = 0.9 [s].
In the text

	Fig. 7. The inaccurate fit of FMM to experiment data for parameters: α = 0.6092, β = 0.0159, E = 6.9425 [MPa] and and τ = 384.7876 [s]
In the text

	Fig. 8. The time course of FMM and the experimental data, parameters: α = 0.6892, β = 0.0139, E = 7.0425 [MPa] and τ = 457.0883[s].
In the text

	Fig. 9. The optimal Maxwell models fit to the experimental data.
In the text