Fractional viscoelastic models for power-law materials

Abstract

Soft materials often exhibit a distinctive power-law viscoelastic response arising from broad distribution of time-scales present in their complex internal structure. A promising tool to accurately describe the rheological behaviour of soft materials is fractional calculus. However, its use in the scientific community remains limited due to the unusual notation and non-trivial properties of fractional operators. This review aims to provide a clear and accessible description of fractional viscoelastic models for a broad audience and to demonstrate the ability of these models to deliver a unified approach for the characterisation of power-law materials. The use of a consistent framework for the analysis of rheological data would help classify the empirical behaviours of soft and biological materials, and better understand their response.

---

1 Introduction

Characterising and understanding how a material deforms when subjected to external forces is critical for many branches of industry, from food and material processing to product design. For example, in additive manufacturing, solid polymer filaments are melted and then extruded through a nozzle. The material flow behaviour during extrusion affects both the processing time required, and the strength of the final printed object.^1,2 In the food industry, the texture and mouthfeel of bread is largely dependent on its mechanical properties, which themselves are largely determined by the manufacturing process.^3,4 Furthermore, the mechanical response of polymers is often used to infer their composition or microstructure.^5,6 The scientific study of material deformation is known as rheology. To further its inquiry, rheology uses mathematical modelling to capture an approximate representation of the behaviour of materials. This facilitates the classification and comparison of materials, and enables predictions that can then inform engineering design choices and improve manufacturing processes.

Studying how living systems respond to mechanical stimuli is fundamental to gaining a comprehensive understanding of their functions and underlying cellular processes. It is increasingly recognized that the mechanical properties of cells have a great impact on human development and disease progression.^14–16 Alterations in the mechanical behaviour of soft tissues when subjected to external stimuli are often associated with diseases.^17–21 Being able to measure and quantify such changes can provide insights into disease evolution, which in turn guide advancements in diagnostic tools and treatments.^22–24 In tissue engineering, knowing the response of human tissue to external forces is key to developing replacements that restore the structure and functionality of damaged tissue by matching the properties of the surrounding environment, which promotes cell in-growth and reduces the risk of implant failure.^25–28

There are two common limits in the behaviour of materials. Elastic materials are solids whose mechanical stress is related to their deformation state with respect to a reference, stress free, geometry; the work done to deform them is entirely stored as recoverable elastic energy. Newtonian fluids are viscous materials whose rate of deformation is a function of their mechanical stress; the mechanical work required to deform the material is fully dissipated. Most materials however exist on a spectrum between the elastic and Newtonian limit behaviours; they are referred to as viscoelastic materials. A defining characteristic of viscoelastic materials is their time-dependent behaviour. When a constant deformation is applied to viscoelastic materials, their internal stress decreases with time – this process is known as relaxation. When viscoelastic materials are subjected to a constant load (stress), their deformation (strain) increases with time – this process is known as creep. In applications where long-term durability and shape stability are required, creep may not be a desirable material property. In contrast, creep may be a desired property in manufacturing processes where a material is extruded or flows across cavities – e.g. in additive manufacturing. Various empirical methods are available to quantify the response of viscoelastic materials. Relaxation tests (in which a constant strain is applied and stress is recorded) and creep tests (in which a constant stress is applied whilst strain is recorded) constitute two testing paradigms for viscoelastic materials. Another common testing methodology requires the application of an oscillatory load to the material. Typical instruments to collect such data are dynamic mechanical analysers (DMA) for materials with solid-like behaviours at long time-scales and rheometers for liquid-like materials that need to be studied in shear deformations only. Relaxation, creep and oscillatory tests are the three canonical viscoelastic testing paradigms; the most appropriate test can depend on several factors including the experimental hardware available, and the most relevant physical modes of deformation for the application or research question.

Mathematical models are commonly used to describe the creep, relaxation and oscillatory behaviour of viscoelastic materials. One of the primary aims of modelling is to extract model parameters and relate them to the underlying molecular or microstructural deformation mechanisms. In other contexts, such as the food industry, the parameters can be related to sensory perception. Nevertheless, the identification of a suitable model facilitates comparison between the behaviour of different materials, states of organisation or environmental conditions. The theory of linear viscoelasticity is a widely used approach to analyse experimental data that yield to a well-defined mathematical representation of the stress–strain–time relation – often referred to as a constitutive relationship. By employing those mathematical methods it is possible to extract parameters that uniquely describe a given material – often referred to as material properties – and allow the prediction of the evolution of stress and strain within the material under arbitrary loading scenarios.

Advances in mechanical testing hardware have revealed that many materials possessing a complex microstructure exhibit a characteristic power-law signature in their creep, relaxation and spectral behaviours. A selection of data from the literature is shown in Fig. 1. Further examples include biological materials – e.g. tissues,^29–31 single cells,^{13,29,32–36} intra and extra cellular components;^10,37,38 gels,^39,40 polymers,^9,41 concrete,^42,43 asphalt,^44–46 ice,⁹ and food – e.g. cheese,^11,47 dough.^8,48 This behaviour represents a challenge to the commonly used viscoelastic models in obtaining a unique mathematical description of the constitutive relation with a limited number of parameters. The power-law behaviour can be approximated with traditional viscoelastic models via a large number of model parameters, which greatly hinders physical interpretation. To circumvent this limitation, traditional models are often discarded in favor of empirical fitting functions: mathematical ansatzes derived from qualitative inspection of the data rather than a governing constitutive equation. This approach precludes comparison of parameters across research studies, which can lead to a multitude of different interpretations of similar phenomena.⁴⁹

Fig. 1 Stress relaxation response of several materials follows a power-law behaviour. Examples of power-law materials include: (a) common polysaccharide (xanthan gum) used as food additive,⁷ bread dough,⁸ synthetic polymers such as nylon of diameter 1.125 mm,⁹ single collagen fibrils;¹⁰ (b) semi-hard zero-fat cheese,¹¹ pure ice at −35 °C,⁹ asphalt sealants¹² and single MDCK epithelial cell.¹³ For clarity the data has been separated into two subplots and the stresses have been normalized by the initial value for each data set.

A promising solution relies on the use of fractional calculus, a branch of mathematics that extends integration and differentiation operators to non-integer order,^50,51 to enrich the classical linear viscoelastic framework.^52,53 This led to the development of a new formalism for the modelling of viscoelastic materials known as fractional viscoelasticity. Fractional viscoelasticity has been applied to complex geological and construction materials such as bitumen (asphalt),^44,45 concrete,^42,43 rock mass,^54–59 waxy crude oil,^60,61 as well as polymers and gels,^41,62–65 and food.¹¹ Numerous examples can also be found of fractional viscoelasticity applied to biological materials such as epithelial cells,⁶⁶ breast tissue cells,^67,68 lung parenchyma,⁶⁹ blood flow,^70,71 as well as red blood cell membranes.⁷² Considering the microscale, some observations of anomalous diffusion in complex fluids have been interpreted based on fractional viscoelastic behaviours,^73–75 with applications to, for instance, transport movements in cell cytosol.^76,77 Acoustic waves are also commonly used to probe or image materials. Fractional calculus appears to be a powerful tool to characterise wave propagation in complex materials exhibiting power law behaviours.^78,79 In spite of the examples reported above, when we consider the number of power-law responses observed in the literature (some of which are modelled empirically), fractional viscoelasticity remains significantly underused. By using a consistent formalism for the analysis of power-law materials, a direct comparison of parameters across studies can be achieved,⁶⁶ which greatly extends the use of available data.

The main aims of this review are to demonstrate the applicability of fractional viscoelasticity to a wide range of real problems, and to carefully illustrate its advantages over commonly used modelling approaches. Furthermore, we aim to present fractional viscoelasticity in a way that is accessible to researchers without specialised knowledge of fractional calculus. To achieve this, we first illustrate the limitations of spring-dashpot models in capturing power-law viscoelasticity. Following this, the mathematical foundations of fractional viscoelasticity are introduced. We will describe the behaviour of networks of fractional viscoelastic elements with a focus on their limit behaviours at short and long time-scales. A number of studies which originally used empirical or spring-dashpot based traditional models are then revisited in order to demonstrate the benefits of fractional viscoelastic models. A physical interpretation of the phenomena is discussed where possible.

2 Introduction to linear viscoelasticity

The linear theory of viscoelasticity provides a powerful mathematical framework to link stress, strain and time and provide predictions of stress and strain distribution during arbitrary loading conditions. A material is assumed linear viscoelastic if the stress function σ(t) and strain function ε(t) are linearly connected: (i) if the strain function ε(t) is multiplied by a constant factor, the resulting stress would be scaled by the same constant, and vice versa. (ii) The response to a strain (or stress) that is a linear combination of two arbitrary strain (or stress) functions is given by the same linear combination of their two individual responses.

Most materials exhibit linear, or quasi linear, behaviour for small deformations while their response becomes nonlinear at large deformations. The linear theory of viscoelasticity is a commonly used approximation for the study of materials’ behaviour as it leads to an accessible and manageable mathematical formulation of the stress, strain and time relationship. Furthermore, the theory of linear viscoelasticity represents a good starting point for the study of more complex nonlinear responses that often gives rise to more complex models.^80–83

A common test performed to characterise a viscoelastic material is to analyse the stress response of the material when subjected to a constant strain (relaxation test). The implication of the linearity assumption is that the resultant stress function scales linearly with the magnitude of the step in strain.⁸⁴ Thus, for a step in strain of amplitude ε₀ at time t = 0, the resulting stress σ(t) can be written

σ(t) = G(t)ε₀, (1)

where G(t) is the relaxation modulus, a monotonically decreasing function. Similarly, in the case of a creep test, the strain response ε(t) is proportional to the stress step amplitude σ₀ that is imposed at t = 0:

ε(t) = J(t)σ₀, (2)

where J(t) is the creep compliance, a monotonically increasing function of time.

A key consequence of the linearity assumption is that the response of the material at time t is given by the sum of the responses to the perturbations imposed at all previous times. Therefore, given an arbitrary stress or strain history, we can compute either the strain response or stress response respectively by integrating the entire history of infinitesimal ‘step’ excitations.^85,86 Assuming that the material does not age, i.e. that its mechanical behaviour does not change with time during the course of a measurement, the strain (or stress) response at time t of a step in stress (or strain) at time τ would be given by the relaxation function G(t − τ) (or creep function J(t − τ)). This results in the following convolution integrals:

Another consequence of linearity is that, when a viscoelastic material is subjected to a sinusoidal stress, the recorded strain has the same frequency, but with a phase difference that may depend on frequency. For purely elastic materials the phase difference is 0° and for purely viscous materials the phase difference is 90°. For viscoelastic materials, the phase difference lies between these two values. If we assume a linear material, by considering an oscillatory excitation ε(t) = e^iωt, we obtain a stress response⁸⁶σ(t) = G*e^iωt from which we can define the so-called complex modulus, or dynamic modulus:

G*(ω) = G′(ω) + iG′′(ω), (5)

where ω is the frequency. The real part of this complex stress response, G′(ω), is defined as the storage modulus (as the energy is stored in an ideal elastic material). The imaginary part of the response, G′′, is defined as the loss modulus (as energy is dissipated in a purely viscous material).

For a linear material, the three moduli, relaxation, creep and dynamic modulus, are directly related. Their relationship can be expressed in the Laplace domain. For instance, by combining the Laplace transforms of eqn (3) and (4), one obtains G̃(s)J̃(s) = s⁻², where G̃(s) and J̃(s) are the Laplace transforms of the relaxation modulus and creep compliance respectively.⁸⁰ Therefore, we can predict one type of behavioural mode (e.g. creep) from observing a different behavioural mode (e.g. relaxation), although such a relationship is difficult to exploit in practice due to the necessarily limited experimental time/frequency domain and the presence of noise on the data.

3 Characteristics of traditional linear viscoelastic models

The expressions above are generic and applicable to any material that behaves in a linear manner, being solid and liquid like. Specific expressions for the moduli – i.e. relaxation, creep and dynamic moduli – can be deduced in several ways. A common approach relies on the use of fundamental viscoelastic units that can be combined in series or parallel – this element-by-element model building process is analogous to electrical circuit models. The most commonly used viscoelastic units are the Hookean spring and the Newtonian dashpot, with governing equations of σ(t) = kε(t) and σ(t) = η[small epsi, Greek, dot above](t) respectively. The simplest models that can be formed from these two viscoelastic units are their series (Maxwell model) and parallel (Kelvin–Voigt model) combinations, both of which are shown in Fig. 2. The creep compliance and relaxation modulus can be derived from the corresponding differential equation by: (a) assuming that the stress is constant and solving for the differential equation in strain, leading to the creep compliance, and (b) assuming that the strain is constant and solving for the differential equation in stress, leading to the relaxation modulus. The complex modulus is derived from the constitutive differential equation, either by considering oscillatory forms of the stress and strain functions, or by taking the Fourier transform of the differential equation, leading to an algebraic equation which can be rearranged for G* = [small sigma, Greek, circumflex](ω)/[small epsilon, Greek, circumflex](ω), where the hat (^∧) symbol denotes a Fourier transformed function.

Fig. 2 Properties of three common viscoelastic material models.

The differential equations relating stress and strain of traditional spring-dashpot models take the form of a linear ordinary differential equation with constant coefficients, and both relaxation modulus and creep compliance involve exponential functions (see Fig. 2). As an example of this, consider the stress response of a Maxwell model to a step in strain:

where G(t) is the relaxation modulus. The ratio τ = η/k has the units of time and can be considered as the characteristic relaxation time of the material. This time-scale is a useful quantity; it represents the time required for the stress to fall to 1/e of its initial value. More importantly, the time-scale may describe a physical relaxation process more effectively than the values of η and k separately, such as rates of binding or unbinding of molecular structures. To capture more complex material behaviours, additional elastic springs and viscous dashpot elements may be used. These additional elements give rise to multiple time-scales. Some models, for example the generalised Maxwell model shown in Fig. 3(a), give rise to an arbitrarily large number of time-scales (see ref. 80, 87 and 88 for more details about spring-dashpot based linear viscoelastic models).

Fig. 3 Use of traditional viscoelastic models to analyze power-law responses. (a) Sketch of the Maxwell model and its generalised form with the correspondent relaxation modulus showing an exponential form. (b) Approximation of a power-law response σ(t) = 5·t^−0.5 (grey line) using the generalised Maxwell model with n = 1 (blue), 2 (green), 4 (dashed purple), where for n = 1 it corresponds to the Maxwell model. Fitted model parameters in Table 2 in Appendix A. (c) Relaxation spectra of the power-law response and the generalised Maxwell model with n = 1, 2, and 4.

Power-law viscoelastic materials correspond to systems with a broad range of relaxation/creep time-scales arising from dissipation/deformation modes occurring at different time and length scales. Mathematically, power-law behaviour consists not of a discrete number of time-scales, but a continuous distribution of time-scales. For this reason, power-law rheology can only be approximated by the spring-dashpot models that yield exponential terms with discrete time-scales.^89–91 The larger the number of exponential terms, the closer the approximation.^92–96 We demonstrate this in Fig. 3 by incrementally increasing the number of Maxwell arms in a generalised Maxwell model, each of which contributes a material time-scale. With four Maxwell arms, a good approximation of a power-law is obtained over the time domain of the data (Fig. 3(b)). Another informative perspective on the approximation process can be gained via a relaxation spectrum, which represents the relative behavioural dominance at different times. The relaxation spectrum found here for the power-law relaxation modulus was approximated as^86,97

where G(t) is the relaxation modulus. Fig. 3(c) shows the relaxation spectrum of the power-law material in logarithmic scale, which is a power-law distribution of time-scales with the same exponent as the original relaxation data. The optimal fit of the various tested Maxwell models can be seen to yield time-scales that are equally spaced in logarithmic scale. As more Maxwell arms are added, a broader coverage of the true power-law distribution is attained. Hence, the quality of the approximation of power-law responses with a generalized Maxwell model depends on both (i) the number of Maxwell components and (ii) the time-window over which the power-law function is fitted.

This qualitative assessment of the above relaxation responses and their corresponding spectra shows that spring-dashpot models might be sufficient to approximate power-law behaviour. However, the approach has a number of disadvantages. A large number of model parameters make the computational analysis far more expensive, and interpretation of their physical meaning becomes significantly more difficult. The characteristic times extracted from the fitting process mostly capture the time window of the data and do not represent intrinsic material properties. Moreover, the fitting would optimise the match with data in the region fitted, but provide poor predictions at time-scales shorter or longer than those provided in the original data.

To circumvent the above disadvantages inherent to spring-dashpot based approximations, empirical power-law expressions have been used as a simple way to capture experimental measurements.^{29,32,33,98,99} These expressions are essentially a mathematical ansatz, formulated to capture the main qualitative features of the experimental data. The end result of this formulation is usually a single modulus of interest that may or may not be analytically converted into the other moduli, necessitating numerical analysis for the prediction of the material response to other external loading patterns. Further, the use of such ad hoc models limits the scope of the measurements since the model parameters may not be easily comparable across studies.

4 Fractional viscoelastic models

With regards to the modelling of power-law viscoelastic materials, we have now outlined the advantages and disadvantages of the phenomenological spring-dashpot approach and the empirical ad hoc approach. In this section we introduce an additional viscoelastic element, the springpot, which is able to capture a power-law behaviour with a minimum number of parameters. We then demonstrate how the springpot can be combined with other elements in series and parallel configurations to capture complex power-law signatures in various contexts.

4.1 The springpot captures power-law behaviour

By inspecting the power-law relaxation data (Fig. 1), it can be seen that the relaxation modulus takes the form G(t) = At^−β, where A is a constant and β the power exponent. Substituting this into eqn (3), the relationship between stress and strain for a power-law material becomes:A branch of mathematics called fractional calculus provides the definition of the generalization of the differentiation operation to non-integer order valid for a function f(t) = 0 for t < 0 as¹⁰⁰where Γ(·) is the gamma function. Redefining the coefficient A = c_β/Γ(1 − β) in eqn (8) yields the following result:which provides a simple constitutive equation for power-law materials. A more detailed derivation of this relationship is presented in the next section.

The springpot's schematic symbol and relationships to the spring and dashpot elements are shown in Fig. 4. For dimensional consistency, the unit of the constant c_β must be Pa s^−β. Due to its unusual physical unit, the parameter c_β lacks a tangible physical meaning, but is often interpreted as the ‘firmness’ of a material.¹⁰¹ It should be noted that in the limiting cases where β = 0 or β = 1, the springpot reduces to a spring or dashpot respectively, and the constant c_β represents the elastic spring constant, k (Pa), or the dashpot viscosity, η (Pa s), respectively.

Fig. 4 Sketch of the fractional element-springpot. It behaves as a spring when β = 0 and as a dashpot when β = 1.

To overcome the use of a constant c_β whose physical meaning may be unclear, the constitutive equation may be written as σ(t) = λ₀τ^β₀d^βε(t)/dt^β, where λ₀ has units of Pa and τ₀ is a characteristic time that has units of s.^102,103 However, it is not possible to design experiments to separately measure λ₀ and τ₀. Therefore, in practical problems, it is more relevant to define the fractional element with two parameters c_β and β,¹⁰⁴ and this is what we use in the current review.

4.1.1 Derivation of the springpot's governing equation. Although not essential for springpot practitioners, it is illuminating to derive the springpot's governing equation. Whilst the correspondence between power-law viscoelasticity and the fractional derivative has been presented before,^105,106 here we attempt to summarize the main steps of the derivation in a pedagogically accessible way.

Fractional derivatives can be understood as generalised integrals. Integrating a function (here ε(t)) an integer n number of times, we get:

where the final simplification on the RHS made use of Cauchy's repeated integral formula. Eqn (11) can in fact be generalised to any positive real number α by recalling the extension of factorial numbers to non-integer values via the Gamma function Γ(α) = (α − 1)!If we substitute α = 1 − β into eqn (12) and , we can then re-write the expression for the total stress reported in eqn (8) as followsSince the fundamental relation (I^a+bf)(t) = I^a(I^bf)(t) holds,¹⁰⁰eqn (13) can be expressed aswhere (I^−βε)(t) is the definition of fractional derivative (D^βε)(t). Equivalent to the above, we can writewhich is the governing equation of the springpot.⁵³

It should be noted that in the literature, different definitions of the fractional derivative have been put forward. Here we use the Caputo derivative definition since it naturally arises from experimental observations and it has shown to have a better applicability to real problems, where initial conditions are known in terms of derivatives of integer index.¹⁰⁷ If we assume that the system is at rest for time t ≤ 0, the fractional derivative is given by¹⁰⁸

where 0 < β < 1 and C denotes a Caputo fractional derivative. So to obtain the value of the Caputo fractional derivative at time t, we must integrate from the initial time t = 0. In other words the fractional derivative is a non-local operator. This is in contrast to integer order derivatives, whose value is determined by the limiting behaviour of the function at the evaluated point. Although hysteresis is also present in the generic hereditary integral shown in eqn (3), it is a fundamental feature of the fractional derivative and consequently the springpot. Lastly, it is important to note that the fractional operator is a linear operator. Therefore, the springpot element lies within the framework of linear viscoelasticity.

4.1.2 Rheological behaviour of the springpot. Despite the simple constitutive equation of the fractional springpot element, it possesses rich behavioural diversity. It is a generalization of the classical viscoelastic elements, the spring and the dashpot, and exhibits a behaviour that is intermediate between the two (see Fig. 5). By substituting a step of deformation ε(t) = ε₀Θ(t) (where Θ(t) is the Heaviside step function) into eqn (15) and making use of the definition of the fractional derivative in eqn (16), we can extract the relaxation modulus G(t) of the springpot:¹⁰⁷The relaxation modulus is as expected a power-law function of the time whose exponent matches the order of the fractional derivative. As a result, the corresponding relaxation spectrum of the springpot is also a simple power-law of exponent −β, and therefore captures very well, by design, the behaviour of power-law materials, with only two parameters. The creep compliance J(t) can be deduced from the expression of G(t). Both are related in the Laplace domain by G̃(s)J̃(s) = s⁻². The Laplace transform of G(t) in eqn (17) is given by G̃(s) = c_βs^β−1. Therefore, J̃(s) = 1/(c_βs^1+β), whose inverse Laplace gives the creep compliance¹⁰⁷

Fig. 5 Responses of a springpot when subjected to (a) constant strain and (b) a step in stress. The red curves are the responses for a spring (β = 0), the blue curves for a dashpot (β = 1) and the grey curves are for increasing values of β from 0.1 to 0.9 (orange, yellow and green curves are respectively β = 0.3, 0.5 and 0.7).

Qualitatively, the linear elastic solid and the Newtonian fluid behaviours are the two limit behaviours of the springpot for β = 0 and 1 respectively. To understand the interpolation of the springpot between elastic and viscous behaviour, we plot the springpot response to a step in strain for different values of β (Fig. 5(a)). When a dashpot (β = 1) is subjected to a step in strain, the stress is initially infinite but immediately dissipated afterwards (Fig. 5(a), blue curve). By gradually decreasing the exponent β the time taken by the springpot to dissipate the stress increases until the limit of the linear elastic solid is reached. When a spring (β = 0) is subjected to a constant deformation, the stress reached is proportional to the spring constant k (Fig. 5(a), red curve).

To illustrate the creep behaviour of the springpot, we consider its response to a step in stress σ₀ imposed for a limited time t* – loading phase – and then returned to 0 – unloading phase. Fig. 5(b) shows plots of the system's response for different values of β during the loading and unloading phases. A spring generates a strain directly proportional to the stress, as expected, and it immediately returns to its original length upon unloading. A dashpot linearly deforms while the stress is imposed, and does not return to its initial state after unloading. The response of the springpot is more complex. By applying the superposition principle, we can write an expression for the springpot response (0 < β < 1) after the load is removed at time t = t* that is given by

noting that when t tends to infinity, the strain tends to zero. Therefore, when the stress is removed, the strain in a springpot is always completely recovered; the springpot has shape memory, although energy would be dissipated during deformation. From eqn (19) one can notice that the recovery time scale during creep experiments depends on both the power-law exponent β and the time for which the load has been imposed t*. The greater the power-law exponent of the springpot β, the larger the time required for total strain recovery, as shown in Fig. 5(b). For example, by looking at the response of a springpot with β = 0.7 (green curve in Fig. 5(b)) over a finite time interval, one might qualitatively interpret the slow recovery process as evidence for plastic behaviour – although this is erroneous because strain eventually converges to zero. This shows that long experiments are needed to fully capture the rheology of power-law materials. It might be useful to re-write eqn (19) as a function of the non-dimensional time τ = t/t*It can be observed that the greater the time for which the load is applied to the material, the longer is the time taken for the strain to be completely recovered. This dependency should be taken into account during the design of experiments to fully capture the rheological response of power-law materials.

It is also illuminating to consider the response of the springpot in the frequency domain. However, the Caputo definition for the fractional derivative is not appropriate here since the assumption that ε(t) = 0 for t < 0 does not hold anymore if we excite the springpot with a periodic deformation. To overcome this problem, we use the Generalised Liouville–Caputo formulation

where L denotes the generalised Liouville–Caputo fractional derivative. Note that this definition is equivalent to the Caputo fractional derivative for functions f such that f(t < 0) = 0.¹⁰⁹ By now considering the oscillatory load ε(t) = e^iωt with eqn (21) (where ^LD^βf(t) = σ(t) and f(t) = ε(t)) the complex modulus of a springpot can be foundwhere ω is the frequency. Separating the real and imaginary parts, we find the storage and loss moduli respectivelyThey follow a power-law behaviour with the same exponent β. When β = 0 the loss modulus is zero as in the case of the spring, whilst when β = 1 the storage modulus is zero as in the case of a dashpot (Fig. 6(a)). The storage and loss modulus exactly match when β = 0.5 (Fig. 6(a)). For β < 0.5 the storage modulus is always greater than the loss modulus, whilst the opposite is true for β > 0.5 (Fig. 6(b)). The phase angle δ between the excitation and the response is related to the storage and loss moduli by tan(δ) = G′′/G′, from which we can derive the retardation phase for a springpot as , constant for all frequencies.

Fig. 6 Storage (solid line) and loss (dash line) moduli for: (a) spring (red curve), dashpot (blue curve) and the special case of β = 0.5 (yellow curves) for which storage and loss moduli are exactly the same; (b) for β < 0.5 the storage modulus is always greater than the loss modulus (β = 0.3, orange curves), whilst the opposite is true for β > 0.5 (β = 0.7, green curves).

4.1.3 Examples of the use of the springpot in practical cases. The first material investigated using fractional calculus was bitumen in 1944, in a study conducted by Scott-Blair and Veinoglou.⁴⁴ Since then, use of the springpot has been somewhat erratic, possibly because of its non-trivial mathematical foundations. However, the springpot has found use in a diverse array of materials outside of the geological and construction-materials contexts; for example, in gels. An early example is the work of Winter and Chambon⁹¹ who derived a springpot-like modulus for crosslinking polymers at their gelation point, which was used for analysis of polydimethylsiloxane gel data. More recently, a springpot has been used in modelling asphalt mixtures⁴⁵ and to capture the oscillatory rheological response of agarose, a polysaccharide extracted from red algae.⁶²

More examples of springpot usage can be found in the context of biological materials more generally. One of the first applications of the springpot in tissue biomechanics was for the study of the viscoelastic properties of lung.¹¹⁰ More recently, it has been used to capture the relaxation behaviour of human arteries¹¹¹ and the dynamic (oscillatory) response of brain tissue for the understanding of neurological disorders.¹¹² There are also several biomechanical studies which make use of a power-law empirical function that corresponds to one of the three springpot moduli discussed above. For example, the empirical function used to analyze the creep response of single cells can be easily related to the springpot's creep compliance.^32,113,114 Similarly, an empirical function used to analyze the relaxation response of smooth muscle cells is equivalent to that of a single springpot.¹¹⁵

4.2 Generalised fractional viscoelastic models

Many materials exhibit power-law viscoelastic behaviour. However, the power-law regime is often limited in time; some materials show a power-law response at short time-scale that converges to a well-defined plateau (solid like behaviour),^29,72 whilst others may exhibit a power-law response followed by a continuous flow of the material (fluid like behaviour).^116,117 In order to capture the diversity of behaviours, ad hoc empirical moduli and large networks of springs and dashpots (i.e. exponential terms) have been often used. However, as previously discussed, these approaches each have disadvantages.

A third useful approach is to combine springpots with themselves and with traditional spring and dashpot elements. In this way, the many advantages of the spring-dashpot network approach can be extended to concisely capture power-law viscoelastic regimes. Such models are often referred to as generalised viscoelastic models and they have been mathematically characterized in the past.^{102,103,118–120} It is worth highlighting that the research field of fractional calculus is currently very active, with recent progress and ongoing work to fully characterise the viscoelastic response of fractional models in terms of novel fractional operators.^79,121,122 The two simplest configurations are two springpots in series or in parallel, which are the fractional analogues of the Maxwell and Kelvin–Voigt models respectively. Examples of the use of the generalized models in practical cases are reported in Table 1 (see ref. 123 for a further review of fractional viscoelasticity in geotechnical engineering). For a comprehensive introduction to fractional calculus and fractional viscoelastic models refer to ref. 78. In the following section, we demonstrate the behavioural diversity that these generalised models are capable of.

Table 1 Examples of practical uses of generalized models

Model	Applications
Fractional Kelvin–Voigt	Modelling of human prostate tissue to develop novel criteria for cancer detection.^125
	Modelling of the oscillatory response of canine liver.^126
	Modelling of the relaxation response of beast cells and tissue samples.^40
	Modelling of the tissue mimicking materials CF-11 and gelatin.^39
	Modelling of the viscoelastic response of potato starch gel.^127
	Modelling of a three-dimensional particulate gel.^63
	Studying of the applicability of the model to wellbore creep^128
Fractional Maxwell	Modelling of the creep behaviour of rocks.^55
	Modelling of the viscoelastic response of tight sandstone.^56
	Modelling of both colloidal^64 and carbopol^129 gel rheology.
	Modelling of arterial tissue.^130
	Modelling the properties of food gels.^11
	Modelling the mechanical properties of collagen gel.^131
Fractional standard linear solid model	Study of single red blood cells.^72
	Modelling of breast cancer cells to develop new diagnostic tool.^68
	Modelling of artery walls for the study of aneurysms.^70,71
	Quantification of changes in viscoelastic response of lung parenchyma due to trauma.^69
	Modelling of the oscillatory response of cancerous cells; the information was then used to selectively attack malignant cells during treatments.^132
	Modelling of the viscoelastic response of brain tissue.^133
	Modelling of the viscoelastic properties of the glassy amorphous polymer poly-methyl-methacrylate.^41
	Modelling of the mechanical properties of polymer networks with addition of fillers.^134
	Modelling of stress relaxation of filled and unfilled polymeric materials.^135
Fractional burgers	Viscoelastic analysis of waxy crude oil.^60,61

---

4.2.1 The fractional Maxwell and Kelvin–Voigt models: heuristic overview. The choice of which viscoelastic model to use is largely informed by the qualitative behaviour of the experimental data. A key advantage of schematically representing viscoelastic models as networks of elements is that it facilitates visual intuition for a model's behaviour. Despite the mathematical complexity of the springpot, limit behaviours can be obtained as commonly done for springs and dashpots models. To address this, we provide below insight on the qualitative behaviour of the two simplest generalised fractional viscoelastic models, which can be extrapolated to more complex models.

When two springpots of exponents α and β (with α > β) are placed in series, the element with the lower exponent (more “springy”) dominates the short time-scale response, while the springpot with the higher exponent determines the long time-scale behaviour (see Fig. 7 first column). This is true for both the relaxation and creep responses, thus for β < α < 1:

When two springpots are combined in parallel the opposite is true: both relaxation and creep responses are dominated by the springpot with higher exponent at short time scale, while the springpot with the lower exponent dominates at long time-scale (see Fig. 7 second column), thus for β < α < 1:Indeed, traditional viscoelastic models can be seen as special cases of the generalised fractional models obtained by specialising the springpots to either springs or dashpots. The qualitative behaviours reported above are consistent with those observed in conventional viscoelastic models, where the Maxwell model shows a liquid-like behaviour at long time scale (dominated by the dashpot, α = 1), while the Kelvin–Voigt model shows a solid-like behaviour at long time-scale (dominated by the spring, β = 0).

Fig. 7 Qualitative behaviour of two springpots in series and parallel. When two springpots are placed in series, the short time scale response is dominated by the springpot with lower exponent, while the long time-scale response is controlled by the springpot with higher exponent (top figures). The contrary is true for two springpots in parallel (bottom figures). In the graphs at the bottom row, the storage modulus is the solid line, while the loss modulus is the dashed line. The parameters of the models are c_α = 1, α = 0.8, c_β = 1, β = 0.2.

Although the above heuristics directly apply to the two-springpot models, the insights can be extrapolated to more complex models in a straightforward manner. For example, consider the standard linear solid model (Fig. 2) but with the spring k₁ replaced by a springpot. If we look at the global network, the k₀ spring in parallel dominates the long-time scale response given its lower power. However, if we focus only on the upper arm, the model will follow the series configuration behaviour. Thus, the short time-scale response is dominated by the springpot and the intermediate response by the dashpot. A more graphical description of this heuristic process has been recently presented by Bonfanti et al.⁶⁶ To build further intuition, readers may also refer to the Annex (ESI† attached to this document), where the relaxation and creep behaviours of more complicated models are briefly reported.

In the following sections we demonstrate the benefits of using fractional models in a wide range of practical applications. We first show that a number of empirical functions previously introduced in the literature actually possess an equivalent fractional model. We then demonstrate that fractional models often capture the rich mechanical behaviour of power-law materials more accurately than spring-dashpot viscoelastic models while using less parameters.

4.2.2 Examples of empirical functions equivalent to fractional viscoelastic models. In this section we show how fractional viscoelastic models can be related to a selection of widely used empirical functions. We limit the study below to functions involving time-invariant parameters. Fractional calculus also provides a solid alternative to a number of empirical models based on non-constant parameters, such as time-dependent viscosity, to account for the complex behaviour of power-law materials.¹³⁶
4.2.2.1 Structural damping model. A number of dynamic (oscillatory) viscoelastic tests on biological tissue were fitted with an empirically derived model known as the structural damping (or hysteretic damping) model^33,137–142 whose complex modulus is given by:where is commonly called the hysteresivity or structural damping coefficient, β is the power-law exponent, G₀ and ω₀ are two scaling factors for stiffness and frequency respectively, and μ the Newtonian viscosity. This modulus can be shown to be exactly equivalent to that of a semi-fractional Kelvin–Voigt model consisting of a springpot in parallel with a dashpot. In fact, eqn (26) can be re-written as:which is equivalent to the complex modulus of a dashpot in parallel to a springpot G*(ω) = c_α(iω)^α + c_β(iω)^β, when α = 1, c_α = μ and c_β = G₀/ω^β₀.
4.2.2.2 Two power-law regime. The frequency response of single cells often exhibits two power-law regimes.^38,124 A lower-power exponent at low frequencies, followed by a higher power-law exponent at high frequencies. Such behaviour has been frequently been analyzed by the empirically derived superposition of two power-laws

G*(ω) = A(iω)^α + B(iω)^β, (28)

where one exponent is often fixed to β = 3/4. The expression above is exactly equivalent to the complex modulus of the fractional Kelvin–Voigt model shown in the inset of Fig. 8 and consisting of two springpots in parallel. Examples of fitting the storage and loss moduli of single cells to the fractional Kelvin–Voigt model are shown in Fig. 8.

Fig. 8 Fitting of storage modulus and loss modulus showing two power-law regimes using fractional Kelvin–Voigt model. Dynamic response of (a) smooth muscle cells cytoskeleton³⁸ and (b) kidney epithelial cells ATP-depleted.¹²⁴ Fitted model parameter values in Table 3.

4.2.2.3 Power-law cut-off. Recently, Khalilgharibi et al.²⁹ studied the relaxation response of epithelial monolayers, a sheet of cells devoid of substrate, to uncover the subcellular components responsible for stress dissipation. The relaxation response consists of an initial power-law phase in the first 5 s, followed by an exponential phase that reaches a plateau at ∼60 s. This biphasic relaxation response was captured by extending the power-law empirical function to include an exponential regime

G(t) = At^−αe^−t/τ + B. (29)

More recently however, traditional viscoelastic elements have been combined with a springpot to characterize the viscoelastic response of epithelial monolayers at both short and long time-scale.⁶⁶ The novel configuration features a springpot that sustains all the load at short time-scale, after which the load is then slowly transferred to a dashpot placed in series which gives rise to the exponential behaviour. This dissipative branch is placed in parallel with a spring to obtain the final plateau (inset in Fig. 9(a)). The mathematical expression for the relaxation modulus of the fractional model introduced by Bonfanti et al.⁶⁶ iswhere c_β and β are the springpot parameters, η the dashpot viscosity, k the spring stiffness, and E_a,b(z) is the Mittag–Leffler function, a special function that often arises from the solution of fractional differential equations^144,145 (see Appendix B). By fitting the fractional model to the empirical functions we can demonstrate that they are approximately equivalent (Fig. 9(a)) and therefore the empirical function can be mapped into the fractional viscoelastic framework.

Fig. 9 Fitting of empirical functions with fractional models. (a) Biphasic empirical function describing the relaxation response of epithelial monolayers fitted with the solid fractional model for different model parameters (fitted models parameters in Table 4). (b) Modified power-law empirical model for the relaxation response fitted with the fractional standard linear solid model (fitted model parameters in Table 5).

Interestingly, the same empirical model was independently developed to capture the relaxation response of one of the most widely used elastomers, polydimethylsiloxane (PDMS),¹⁴³ and the rheological behaviour of alginate-based gels.^146,147 The study of the liquid–solid transition during gelation has attracted significant attention in the past.¹⁴⁸ Friedrich et al.¹⁴³ reported the rheological response of stoichiometrically balanced polydimethylsiloxane (PDMS) in which the cross-linking reaction was interrupted at different times, before and after the critical point (referred to as the gel point). As shown in Fig. 10, the four-parameter fractional solid model can also be applied here to accurately captures the behaviour of the time-evolving cross-linking reaction pre- and post-gelation.

Fig. 10 Storage (blue) and loss (red) moduli of cross-linked PDMS samples at different extent of reactions (t − t_c, where t is the current time and t_c the time of gelation). The x-axis is shifted by A decades to facilitate reading. The dots are the experimental data (from ref. 143), while the dashed lines represent the fitted fractional viscoelastic model shown in the inset. Fitted parameters reported in Table 7 in Appendix A.

Physically, during gelation a polymer network forms. The critical gel point is associated with an abrupt change in viscosity,¹⁴⁹ which is defined as the creation of the first percolation cluster that spans the sample. By observing the values of the fitted parameters reported in Table 7 in Appendix A, the fluid-like behaviour of the PDMS pre-gel is confirmed by the zero value of the stiffness of the spring in parallel. As expected from physical considerations, at the gel point we observe a rapid increase of the viscosity (Fig. 12 in Appendix A). Post-gel the PDMS behaves as a solid, which is confirmed by the rapid increase of the stiffness k. As discussed above, the storage and loss moduli of the springpot are exactly equal when β = 0.5, or physically, when the springpot is exactly intermediate between a liquid and a solid. The physical relevance of this parameter is confirmed from the fitted springpot coefficient shown in Table 7, where β at time t − t_c = 0 (gel point) is ∼0.5.

4.2.2.4 Modified power-law models. Another empirical model, recently used for the analysis of the relaxation response of various benign and malignant cell lines, is the modified power-law (MPL) model⁹⁵ which can be written aswhere E₀ is the instantaneous or ‘glassy’ modulus, E_∞ is the plateau or ‘rubbery’ modulus, τ is a time scaling, and α determines the power-law gradient of relaxation. In contrast to a regular power-law model, the MPL model has the convenient property of being well defined at short time scales. Its behaviour is notably similar to a fractional standard linear solid model (a springpot in series with a spring, all in parallel to a spring), as shown in Fig. 9(b). The relaxation modulus of the fractional standard linear solid model is defined aswhere k_β and k_γ are the two spring constants, c_α′ and α′ are the springpot parameters (with a prime added to α to avoid confusion with the MPL α parameter), and E is the Mittag–Leffler function.^144,145 The similarity between the two models was also discussed by Bagley,¹⁵⁰ who made comparisons between their relaxation spectra. Here we show their similarity by directly matching their boundary condition parameters at t = 0 and t → ∞ and fitting the fractional standard linear solid parameters c_α′ and α′ to the MPL (Fig. 9(b)). The fractional model fits the MPL well, especially at longer time scales where both the Mittag–Leffler function and the MPL asymptotically approach a simple power-law.¹⁵¹ Interestingly, the MPL model has also been used in several asphalt and asphalt-concrete studies^46,152,153 indicating that the fractional standard linear solid may have utility in this field.

4.2.3 Fractional models capture complex power-law behaviours with less parameters. As discussed, spring-dashpot viscoelastic models can be used to capture power-law viscoelasticity but this results in an abundance of model parameters. We have shown that the springpot provides an efficient way to describe the rheology of power-law materials with only two parameters. This ability relies on the fact that while spring-dashpot models define specific time scales, the fractional viscoelastic element directly parameterises the relaxation spectrum, and thus the distribution of characteristic time scales. To further illustrate the ability of fractional models to accurately capture complex power-law responses with a minimum number of model parameters, we now re-analyse the time-response of different power-law materials, originally modelled using spring-dashpot networks, by using fractional viscoelastic models.

We first analyse the rheological behaviour of single cells presented by Darling et al.³⁴ They investigated the relaxation response of zonal articular chondrocytes by the use of AFM (Atomic Force Microscopy). The standard linear solid model originally used in the paper struggles to capture the short-time scale response (see Fig. 11(a)). We propose the use of the fractional Kelvin–Voigt model (springpot in parallel to a spring) that involves the same number of parameters (see inset in Fig. 11(a)). A cursory look at the comparison between the fits of traditional and fractional viscoelastic models to the relaxation response of zonal articular chondrocytes shown in Fig. 11(a) reveals the ability of the fractional model to significantly improve the quality of the fitting by accurately capturing the fast relaxation response while using the same number of parameters as the spring-dashpot model used in the original paper.

Fig. 11 (a) Relaxation data from zonal articular chondrocytes³⁴ fitted with a standard linear solid and fractional Kelvin–Voigt specialized by one spring. (b) Relaxation data from tomato mesocarp cells¹⁵⁴ fitted with a 2 time-scale generalised Maxwell model and a fractional Maxwell model specialized by one dash-pot. (c) Relaxation data from PCL/bio-active glass¹⁵⁵ fitted with a 2 time-scale generalised Maxwell model and a fractional special model. (d) Relaxation data from collagen fibrils¹⁰ fitted with a 2 time-scale generalised Maxwell model and a fractional special model. All fitted parameters are reported in Table 6.

Recently, Zhang et al.⁶⁰ tested tomato mesocarp cells under high speed microcompression to understand how industrial handling may affect the integrity of fresh fruits. To successfully capture the distribution of time-scales that gives rise to the power-law regime at short time-scale, a 2 time-scale standard linear solid that requires five parameters was originally used. By using a fractional Maxwell model (springpot in series to a dashpot) we achieve the same excellent fit to the data whilst reducing the number of parameters to three (Fig. 11(b)).

The 2 time-scale standard linear solid model has recently been applied to the analysis of the rheological behaviour of other materials, such as polycaprolactone bioactive glass tested under compression¹⁵⁵ and single collagen fibrils from the extracellular matrix under tensile testing.¹⁰ Given that they present the same qualitative behaviour as the epithelial monolayers (power-law followed by exponential behaviour until a final steady-state value is reached), we have successfully applied the fractional model recently developed by Bonfanti et al.⁶⁶ to data from these other materials. The model accurately fits both short and long time-scale responses while using one less parameter (Fig. 11(c and d)).

5 Conclusions

This review has demonstrated the ability of fractional viscoelastic models to accurately capture the rheological responses of a broad range of materials while using less parameters than traditional viscoelastic models, identifying material parameters that account for the wide distribution of time-scales involved in power-law behaviour. Despite the limitation to the linear response, fractional models exhibit rich behaviours consistent with empirical data in both relaxation and creep experiments. Looking at large deformations and failure of complex materials remains however a challenge that cannot be tackled by linear rheological models and we anticipate significant efforts to tackle these questions in future.

The main impediments to the use and dissemination of fractional viscoelastic models appears to be their relatively intricate mathematical formalism and the difficulty of their numerical implementation. This document illustrated the qualitative response of the springpot element and its composition into simple networks in series and in parallel configurations, helping the reader to build intuition. An Annex to the review (ESI†) provides an exhaustive list of rheological models including up to three elements, alongside the analytical expression of their moduli and graphs illustrating their behaviour. We also developed a software library, RHEOS, to facilitate the fitting of experimental data and prediction of power-law behaviours.¹⁵⁶ The software package allows non-experts to fit their own data using a wide selection of viscoelastic models, accounting for complex loading patterns (all figures of this work have been created using RHEOS). Altogether, these elements significantly lower the barriers to using fractional models to analyse the rheology of power-law materials.

Although fractional models extend the range of behaviours that can be modelled with rheological elements, they do not provide on their own explanations for the underlying mechanisms giving rise to the observed macroscopic power-law behaviour. A deeper understanding of soft material mechanics would require a more systematic theoretical analysis of experimental data that would provide a physical underpinning for the emergence of power-law behaviours. However, fractional models provide systematic approaches to capture material parameters that can be compared across studies, and such comparison is likely to provide a better handle on the underlying physical significance of power-law behaviours.

Conflicts of interest

There are no conflicts to declare.

Appendix A: fitted parameters

Table 2 Fitted parameters generalised Maxwell models from Fig. 3(b)

Parameters	1 time-scales	2 time-scales	4 time-scales
k 1 (Pa)	3.25	11.50	7.50
η 1 (Pa s)	86.57	10.12	1.23
k 2 (Pa)	—	1.88	4.71
η 2 (Pa s)	—	103.64	5.55
k 3 (Pa)	—	—	1.96
η 3 (Pa s)	—	—	17.14
k 4 (Pa)	—	—	1.06
η 4 (Pa s)	—	—	132.97

---

Table 3 Fitted parameters fractional Kelvin–Voigt model from Fig. 8. Note that the fitted storage and loss moduli of TC7 kidney epithelial cells are normalized with respect to their values at ω = 10 Hz

Data	c α	α	c β	β
Bovine trachea smooth muscle cells^38	0.01 (Pa nm^−1 s^−α)	0.78	0.9 (Pa nm^−1 s^−β)	0.1
TC7 kidney epithelial cells ATP-depleted^124	0.02 (s^α)	0.6	0.7 (s^β)	0.07

---

Table 4 Parameters of empirical function and three-element fractional model from Fig. 9(a)

	Empirical model		Fractional model
Curve 1 (bottom)	A (Pa s^β)	960	c β (Pa s^β)	1125
	α	0.15	β	0.1
	τ (s)	8.0	η (Pa s)	6650
	B (Pa)	800	k (Pa s)	800
Curve 2	A (Pa s^β)	960	c β (Pa s^β)	1400
	α	0.3	β	0.21
	τ (s)	8.0	η (Pa s)	6100
	B (Pa)	1300	k (Pa)	1300
Curve 3	A (Pa s^β)	960	c β (Pa s^β)	6000
	α	0.5	β	0.4
	τ (s)	8.0	η (Pa s)	2080
	B (Pa)	1800	k (Pa)	1800
Curve 4 (top)	A (Pa s^β)	960	c β (Pa s^β)	6700
	α	0.7	β	0.54
	τ (s)	8.0	η (Pa s)	3040
	B (Pa)	2300	k (Pa)	2300

---

Table 5 Parameters of empirical function and three-element fractional model from Fig. 9(b)

	Empirical model		Fractional model
Curve 1 (bottom)	E ∞ (Pa)	0.5	c α′ (Pa s^α)	0.007
	E 0 (Pa)	1.0	α′	0.72
	τ (s)	1 × 10^−3	k 1 (Pa)	0.5
	α	0.8	k 0 (Pa)	0.5
Curve 2	E ∞ (Pa)	0.5	c α′ (Pa s^α)	0.017
	E 0 (Pa)	1.0	α′	0.59
	τ (s)	1 × 10^−3	k 1 (Pa)	0.5
	α	0.6	k 0 (Pa)	0.5
Curve 3	E ∞ (Pa)	0.5	c α′ (Pa s^α)	0.051
	E 0 (Pa)	1.0	α′	0.42
	τ (s)	1 × 10^−3	k 1 (Pa)	0.5
	α	0.4	k 0 (Pa)	0.5
Curve 4 (top)	E ∞ (Pa)	0.5	c α′ (Pa s^α)	0.096
	E 0 (Pa)	1.0	α′	0.33
	τ (s)	1 × 10^−3	k 1 (Pa)	0.5
	α	0.3	k 0 (Pa)	0.5

---

Table 6 Fitted parameters of the traditional viscoelastic models and the fractional model from Fig. 11. Note that the fitted stress relaxation responses of PCL/bio-active glass are normalised with respect to the maximum stress

Material	Traditional model		Fractional model
Zonal articular chondrocytes^34	k 0 (nN)	260	c β (nN s^β)	420
	η 1 (nN s)	172	β	0.82
	k 1 (nN)	260	k (nN)	253
Tomato mesocarp cells^154	k 0 (Pa)	8 × 10^5	η (Pa s)	6 × 10^6
	η 1 (Pa s)	1 × 10^6	c β (Pa s^β)	1.5 × 10^6
	k 1 (Pa)	6 × 10^5	β	0.16
	η 2 (Pa s)	1 × 10^6
	k 2 (Pa)	4 × 10^4
PCL/bio-active glass^155 bottom sample	k 0	0.8	η (s)	1.6
	k 1	0.1	c β (s^β)	0.15
	η 1 (s)	0.01	β	0.02
	k 2	0.15	k	0.76
	η 2 (s)	1.46
PCL/bio-active glass^155 top sample	k 0	0.8	η (s)	1.5
	k 1	0.1	c β (s^β)	0.13
	η 1 (s)	0.04	β	0.13
	k 2	0.1	k	0.8
	η 2 (s)	1.02
Collagen fibrils^10 top sample	k 0 (MPa)	473	η (MPa s)	4700
	k 1 (MPa)	93	c β (MPa s^β)	165
	η 1 (MPa s)	3800	β	0.12
	k 2 (MPa)	80	k (MPa)	470
	η 2 (MPa s)	300
Collagen fibrils^10 middle sample	k 0 (MPa)	414	η (MPa s)	5100
	k 1 (MPa)	78	c β (MPa s^β)	162
	η 1 (MPa s)	331	β	0.21
	k 2 (MPa)	62	k (MPa)	408
	η 2 (MPa s)	2700
Collagen fibrils^10 bottom sample	k 0 (MPa)	395	η (MPa s)	4400
	k 1 (MPa)	75	c β (MPa s^β)	120
	η 1 (MPa s)	3080	β	0.12
	k 2 (MPa)	41	k (MPa)	390
	η 2 (MPa s)	163

---

Table 7 Fitted parameters of the fractional model from Fig. 10

t − tc	η (Pa)	c β (Pa s^β)	β	k (Pa)
−6	3.9	9.4	0.71	0.0
−2	52.3	50.5	0.61	0.0
0	336.3	121.0	0.54	20.2
2	1122.7	225.2	0.46	101.3
6	2963.5	920.6	0.36	1923.0

---

Fig. 12 Fitted parameters of the fractional model from Fig. 10.

Appendix B: Mittag–Leffler function

The Mittag–Leffler function arises frequently in the study of fractional viscoelasticity (and fractional calculus more generally). This appendix aims to provide an overview of the function's behavior. There are one, two and three parameter versions of the Mittag–Leffler function;¹⁵⁷ only the first two are introduced here as they are most relevant. The function has a series definition and an integral definition.¹⁴⁵ The series definition is more amenable to intuition and will be used here. The one parameter series representation is:The first special case of note is E₀(z) = 1/(1 − z), which can be confirmed by Taylor expanding the fraction about z = 0 and comparing with the series definition in eqn (33), noting that Γ(1) = 1. The second special case of note is E₁(z) = e^z, which can be seen by comparing eqn (33) to the series definition of the exponential function and noting that Γ(k + 1) = k! for k ∈ [Doublestruck N]⁰. These two special cases illuminate the fact that the Mittag–Leffler function is capable of interpolating between power law and exponential behaviour. The two asymptotic representations are¹⁵¹For values of 0 < α < 1 the behaviour as t → 0 approaches a stretched exponential, while for t → ∞ the Mittag–Leffler function is asymptotically equivalent to a power-law (Fig. 13).

Fig. 13 Qualitative behaviour of the Mittag–Leffler function for different values of α = 0.25, 0.50, 0.75.

The one parameter Mittag–Leffler function is in fact a special case of the two parameter function, E_α,1 = E_α. The series representation of the two parameter Mittag–Leffler function is

As might be expected, the two parameter Mittag–Leffler function is capable of even greater behavioural diversity than the one parameter version. Although the manifold behaviours are not summarised here as they lie beyond the scope of this review, the asymptotic behaviour of the two parameter Mittag–Leffler function can be found in a recent paper.¹⁵⁸

Acknowledgements

The authors wish to acknowledge present and past members of the Kabla and Charras labs for stimulating discussions. The authors would like to thank Prof. Manuel Duarte Ortigueira, Prof. Sverre Holm, Prof. Jordan Hristov, Dr Arran Ferandez and the referees for their comments on the pre-print version of this work. AB, JLK, and AJK acknowledge the BBSRC grants BB/M002578/1, BB/K018175/1, and BB/P003184/1. JLK would like thank the George and Lillian Schiff Foundation for the PhD funding which facilitated this project. G. C. was supported by a consolidator grant from the European Research Council (MolCellTissMech, agreement 647186) and by the BBSRC grant BB/M003280.

Notes and references

1. M. E. Mackay, J. Rheol., 2018, 62, 1549–1561.
2. A. Corker, H. C.-H. Ng, R. J. Poole and E. Garca-Tuñón, Soft Matter, 2019, 15, 1444–1456.
3. R. I. Tanner, F. Qi and S.-C. Dai, J. Non-Newtonian Fluid Mech., 2008, 148, 33–40.
4. A. Lazaridou, D. Duta, M. Papageorgiou, N. Belc and C. G. Biliaderis, J. Food Eng., 2007, 79, 1033–1047.
5. N. Hata, A. Tobolsky and A. Bondi, J. Appl. Polym. Sci., 1968, 12, 2597–2613.
6. E. van Ruymbeke, S. Coppola, L. Balacca, S. Righi and D. Vlassopoulos, J. Rheol., 2010, 54, 507–538.
7. D. Gagnon, X. Shen and P. Arratia, EPL, 2013, 104, 14004.
8. J. Figueroa, C. Manuel, Z. Hernández-Estrada and B. Ramrez-Wong, Cereal Chem., 2012, 89, 211–216.
9. X. Xu and J. Hou, Mech. Time-Depend. Mater., 2011, 15, 29–39.
10. Z. L. Shen, H. Kahn, R. Ballarini and S. J. Eppell, Biophys. J., 2011, 100, 3008–3015.
11. T. Faber, A. Jaishankar and G. McKinley, Food Hydrocolloids, 2017, 62, 325–339.
12. S. Liu, L. Mo, K. Wang, Y. Xie and M. Woldekidan, Constr. Build. Mater., 2016, 105, 1–13.
13. N. Desprat, A. Guiroy and A. Asnacios, Rev. Sci. Instrum., 2006, 77, 055111.
14. Y. Park, C. A. Best, K. Badizadegan, R. R. Dasari, M. S. Feld, T. Kuriabova, M. L. Henle, A. J. Levine and G. Popescu, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 6731–6736.
15. Y. Sun, L. G. Villa-Diaz, R. H. Lam, W. Chen, P. H. Krebsbach and J. Fu, PLoS One, 2012, 7, e37178.
16. H. A. Messal, S. Alt, R. M. Ferreira, C. Gribben, V. M.-Y. Wang, C. G. Cotoi, G. Salbreux and A. Behrens, Nature, 2019, 566, 126.
17. K. S. Cunningham and A. I. Gotlieb, Lab. Invest., 2005, 85, 9.
18. S. Phipps, T. Yang, F. Habib, R. Reuben and S. McNeill, Urology, 2005, 66, 447–450.
19. D. A. Moulding, E. Moeendarbary, L. Valon, J. Record, G. T. Charras and A. J. Thrasher, Blood, 2012, 120(18), 3803–3811.
20. E. S. White, Ann. Am. Thorac. Soc., 2015, 12, S30–S33.
21. A. Fritsch, M. Höckel, T. Kiessling, K. D. Nnetu, F. Wetzel, M. Zink and J. A. Käs, Nat. Phys., 2010, 6, 730.
22. J. Palacio-Torralba, S. Hammer, D. W. Good, S. A. McNeill, G. D. Stewart, R. L. Reuben and Y. Chen, J. Mech. Behav. Biomed. Mater., 2015, 41, 149–160.
23. T. W. Remmerbach, F. Wottawah, J. Dietrich, B. Lincoln, C. Wittekind and J. Guck, Cancer Res., 2009, 69, 1728–1732.
24. F. Rahimov and L. M. Kunkel, J. Cell Biol., 2013, 201, 499–510.
25. B. Chan and K. Leong, Eur. Spine J., 2008, 17, 467–479.
26. L. Nimeskern, G. J. van Osch, R. Müller and K. S. Stok, Tissue Eng., Part B, 2013, 20, 17–27.
27. F. Guilak, D. L. Butler, S. A. Goldstein and F. P. Baaijens, J. Biomech., 2014, 47, 1933–1940.
28. M. M. Laronda, A. L. Rutz, S. Xiao, K. A. Whelan, F. E. Duncan, E. W. Roth, T. K. Woodruff and R. N. Shah, Nat. Commun., 2017, 8, 15261.
29. N. Khalilgharibi, J. Fouchard, N. Asadipour, R. Barrientos, M. Duda, A. Bonfanti, A. Yonis, A. Harris, P. Mosaffa and Y. Fujita, et al. , Nat. Phys., 2019, 15(8), 839–847.
30. F. Serwane, A. Mongera, P. Rowghanian, D. A. Kealhofer, A. A. Lucio, Z. M. Hockenbery and O. Campàs, Nat. Methods, 2017, 14, 181.
31. S. Nicolle, P. Vezin and J.-F. Palierne, J. Biomech., 2010, 43, 927–932.
32. N. Desprat, A. Richert, J. Simeon and A. Asnacios, Biophys. J., 2005, 88, 2224–2233.
33. B. Fabry, G. N. Maksym, J. P. Butler, M. Glogauer, D. Navajas and J. J. Fredberg, Phys. Rev. Lett., 2001, 87, 148102.
34. E. Darling, S. Zauscher and F. Guilak, Osteoarthr. Cartil., 2006, 14, 571–579.
35. X. Trepat, M. Grabulosa, F. Puig, G. N. Maksym, D. Navajas and R. Farré, Am. J. Physiol.: Lung Cell. Mol. Physiol., 2004, 287, L1025–L1034.
36. A. E. Ekpenyong, G. Whyte, K. Chalut, S. Pagliara, F. Lautenschläger, C. Fiddler, S. Paschke, U. F. Keyser, E. R. Chilvers and J. Guck, PLoS One, 2012, 7, e45237.
37. E. Fischer-Friedrich, Y. Toyoda, C. J. Cattin, D. J. Müller, A. A. Hyman and F. Jülicher, Biophys. J., 2016, 111, 589–600.
38. L. Deng, X. Trepat, J. P. Butler, E. Millet, K. G. Morgan, D. A. Weitz and J. J. Fredberg, Nat. Mater., 2006, 5, 636.
39. F. Meral, T. Royston and R. Magin, Commun. Nonlinear Sci. Numer. Simul., 2010, 15, 939–945.
40. H. Zhang, Q. Zhe Zhang, L. Ruan, J. Duan, M. Wan and M. F. Insana, Meas. Sci. Technol., 2018, 29, 035701.
41. M. Alcoutlabi and J. Martinez-Vega, Polymer, 1998, 39, 6269–6277.
42. Y. Bouras, D. Zorica, T. M. Atanacković and Z. Vrcelj, Appl. Math. Modell., 2018, 55, 551–568.
43. F. Barpi and S. Valente, Eng. Fract. Mech., 2002, 70, 611–623.
44. G. S. Blair and B. Veinoglou, J. Sci. Instrum., 1944, 21, 149.
45. G. D. Mino, G. Airey, M. Di Paola, F. P. Pinnola, G. D’Angelo and D. L. Presti, J. Civil Eng. Manag., 2016, 22, 882–889.
46. Y. Kim, Y. Lee and H. Lee, J. Mater. Civil Eng., 1995, 7, 59–68.
47. R. Subramanian, K. Muthukumarappan and S. Gunasekaran, Int. J. Food Prop., 2006, 9, 377–393.
48. J. Lefebvre and N. Mahmoudi, J. Cereal Sci., 2007, 45, 49–58.
49. P.-H. Wu, D. R.-B. Aroush, A. Asnacios, W.-C. Chen, M. E. Dokukin, B. L. Doss, P. Durand-Smet, A. Ekpenyong, J. Guck and N. V. Guz, et al. , Nat. Methods, 2018, 15, 491–498.
50. H. Rudolf, Applications of fractional calculus in physics, World Scientific, 2000.
51. R. L. Magin, Fractional calculus in bioengineering, Begell House Redding, 2006.
52. G. S. Blair, B. Veinoglou and J. Caffyn, Proc. R. Soc. London, Ser. A, 1947, 189, 69–87.
53. G. S. Blair, J. Colloid Sci., 1947, 2, 21–32.
54. H. Zhou, C. Wang, B. Han and Z. Duan, Int. J. Rock Mech. Min. Sci., 2011, 48, 116–121.
55. F. Wu, J. F. Liu and J. Wang, Environ. Earth Sci., 2015, 73, 6965–6971.
56. X. Ding, G. Zhang, B. Zhao and Y. Wang, Sci. Rep., 2017, 7, 11336.
57. M. Liao, Y. Lai, E. Liu and X. Wan, Acta Geotech., 2017, 12, 377–389.
58. C.-C. Zhang, H.-H. Zhu, B. Shi and B. Fatahi, Comput. Geotech., 2018, 94, 72–82.
59. C. Ma, H.-b. Zhan, W.-m. Yao and H.-z. Li, Water Sci. Eng., 2018, 11, 131–138.
60. Z. Zhang, L. Hou, X. Chen, Y. Zhou, M. Liu and W. Zhou, J. Dispersion Sci. Technol., 2016, 37, 326–332.
61. L. Hou, C. Song, W. Yan and Z. Zhang, Rheol. Acta, 2014, 53, 349–356.
62. Q. Chen, B. Suki and K.-N. An, J. Biomech. Eng., 2004, 126, 666–671.
63. M. Bouzid, B. Keshavarz, M. Geri, T. Divoux, E. Del Gado and G. H. McKinley, J. Rheol., 2018, 62, 1037–1050.
64. S. Aime, L. Cipelletti and L. Ramos, J. Rheol., 2018, 62, 1429–1441.
65. J. L. Kaplan, T. A. Torode, F. B. Daher and S. A. Braybrook, bioRxiv, 2019, 565614.
66. A. Bonfanti, J. Fouchard, N. Khalilgharibi, G. Charras and A. Kabla, R. Soc. Open Sci., 2020, 7, 190920.
67. C. Coussot, S. Kalyanam, R. Yapp and M. F. Insana, IEEE Trans. Ultrason. Eng., 2009, 56, 715–726.
68. B. Carmichael, H. Babahosseini, S. Mahmoodi and M. Agah, Phys. Biol., 2015, 12, 046001.
69. Z. Dai, Y. Peng, H. A. Mansy, R. H. Sandler and T. J. Royston, Med. Eng. Phys., 2015, 37, 752–758.
70. P. Perdikaris and G. E. Karniadakis, Ann. Biomed. Eng., 2014, 42, 1012–1023.
71. Y. Yu, P. Perdikaris and G. E. Karniadakis, J. Comput. Phys., 2016, 323, 219–242.
72. D. Craiem and R. L. Magin, Phys. Biol., 2010, 7, 013001.
73. I. Goychuk, Adv. Chem. Phys., 2012, 150, 187.
74. J.-H. Jeon and R. Metzler, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 021147.
75. J.-H. Jeon, N. Leijnse, L. B. Oddershede and R. Metzler, New J. Phys., 2013, 15, 045011.
76. G. Guigas, C. Kalla and M. Weiss, Biophys. J., 2007, 93, 316–323.
77. T. J. Lampo, S. Stylianidou, M. P. Backlund, P. A. Wiggins and A. J. Spakowitz, Biophys. J., 2017, 112, 532–542.
78. F. Mainardi, Fractional calculus and waves in linear viscoelasticity: an introduction to mathematical models, World Scientific, 2010.
79. S. Holm, Waves with Power-Law Attenuation, Springer, 2019.
80. W. N. Findley and F. A. Davis, Creep and relaxation of nonlinear viscoelastic materials, Courier Corporation, 2013.
81. B. Keshavarz, T. Divoux, S. Manneville and G. H. McKinley, ACS Macro Lett., 2017, 6, 663–667.
82. N. A. Bharadwaj, K. S. Schweizer and R. H. Ewoldt, J. Rheol., 2017, 61, 643–665.
83. S. Müller, M. Kästner, J. Brummund and V. Ulbricht, Comput. Mater. Sci., 2011, 50, 2938–2949.
84. T. Kailath, Linear systems, Prentice-Hall Englewood Cliffs, NJ, 1980, vol. 156.
85. D. Gutierrez-Lemini, Engineering Viscoelasticity, Springer, 2014, pp. 23–52.
86. R. Lakes and R. S. Lakes, Viscoelastic materials, Cambridge University Press, 2009.
87. D. Roylance, Department of Materials Science and Engineering-Massachusetts Institute of Technology, Cambridge, MA, 2001, vol. 2139, pp. 1–37.
88. D. R. Bland, The theory of linear viscoelasticity, Courier Dover Publications, 2016.
89. T. Svensson, IEEE Trans. Circuit Theory, 1973, 20, 142–144.
90. G. Beylkin and L. Monzón, Appl. Comput. Harmon. Anal., 2005, 19, 17–48.
91. H. H. Winter and F. Chambon, J. Rheol., 1986, 30, 367–382.
92. N. Heymans and J.-C. Bauwens, Rheol. Acta, 1994, 33, 210–219.
93. A. Halldin, M. Ander, M. Jacobsson and S. Hansson, World J. Mech., 2014, 4, 348.
94. A. Bembey, M. Oyen, A. Bushby and A. Boyde, Philos. Mag., 2006, 86, 5691–5703.
95. Y. M. Efremov, W.-H. Wang, S. D. Hardy, R. L. Geahlen and A. Raman, Sci. Rep., 2017, 7, 1541.
96. B. Wang, W. Wang, Y. Wang and L. Liu, Nano/Molecular Medicine and Engineering (NANOMED), 2016 IEEE 10th International Conference on, 2016, pp. 51–54.
97. R. Christensen, Theory of viscoelasticity: an introduction, Elsevier, 2012.
98. Y. Kobayashi, M. Tsukune, T. Miyashita and M. G. Fujie, Phys. Rev. E, 2017, 95, 022418.
99. N. Bonakdar, R. Gerum, M. Kuhn, M. Spörrer, A. Lippert, W. Schneider, K. E. Aifantis and B. Fabry, Nat. Mater., 2016, 15, 1090.
100. K. Oldham and J. Spanier, The fractional calculus theory and applications of differentiation and integration to arbitrary order, Elsevier, 1974, vol. 111.
101. G. S. Blair and F. Coppen, Am. J. Psychol., 1942, 215–229.
102. H. Schiessel, R. Metzler, A. Blumen and T. Nonnenmacher, J. Phys. A: Math. Gen., 1995, 28, 6567.
103. F. Mainardi and G. Spada, Eur. Phys. J.-Spec. Top., 2011, 193, 133–160.
104. A. Jaishankar and G. H. McKinley, Proc. R. Soc. London, Ser. A, 2013, 20120284.
105. I. Podlubny, Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications, Elsevier, 1998, vol. 198.
106. T. Surguladze, J. Math. Sci., 2002, 112, 4517–4557.
107. M. Di Paola, A. Pirrotta and A. Valenza, Mech. Mater., 2011, 43, 799–806.
108. R. Gorenflo and F. Mainardi, Fractals and Fractional Calculus in Continuum Mechanics, 2008, vol. 1, pp. 223–276.
109. M. Ortigueira and D. Valério, Fractional Signals and Systems, Walter de Gruyter GmbH, 2020.
110. B. Suki, A.-L. Barabasi and K. R. Lutchen, J. Appl. Physiol., 1994, 76, 2749–2759.
111. D. Craiem, F. J. Rojo, J. M. Atienza, R. L. Armentano and G. V. Guinea, Phys. Med. Biol., 2008, 53, 4543.
112. E. Nicolas, S. Calle, S. Nicolle, D. Mitton and J.-P. Remenieras, Ultrasonics, 2018, 84, 119–125.
113. E. Zhou, S. Quek and C. Lim, Biomech. Model. Mechanobiol., 2010, 9, 563–572.
114. P. A. Pullarkat, P. A. Fernández and A. Ott, Phys. Rep., 2007, 449, 29–53.
115. J. D. Hemmer, J. Nagatomi, S. T. Wood, A. A. Vertegel, D. Dean and M. LaBerge, J. Biomech. Eng., 2009, 131, 041001.
116. D. Tripathi, S. Pandey and S. Das, Appl. Math. Comput., 2010, 215, 3645–3654.
117. C. Celauro, C. Fecarotti, A. Pirrotta and A. Collop, Constr. Build. Mater., 2012, 36, 458–466.
118. W. G. Gloeckle and T. F. Nonnenmacher, Macromolecules, 1991, 24, 6426–6434.
119. W. Glöckle and T. F. Nonnenmacher, Rheol. Acta, 1994, 33, 337–343.
120. R. H. Pritchard and E. M. Terentjev, J. Rheol., 2017, 61, 187–203.
121. J. Y. Hristov, Front. Phys., 2018, 6, 135.
122. J. Hristov, Eur. Phys. J. Plus, 2019, 134, 283.
123. J. Lai, S. Mao, J. Qiu, H. Fan, Q. Zhang, Z. Hu and J. Chen, Math. Probl. Eng., 2016, 2016, 1–15.
124. B. D. Hoffman, G. Massiera, K. M. Van Citters and J. C. Crocker, Proc. Natl. Acad. Sci. U. S. A., 2006, 103, 10259–10264.
125. M. Zhang, P. Nigwekar, B. Castaneda, K. Hoyt, J. V. Joseph, A. di Sant’Agnese, E. M. Messing, J. G. Strang, D. J. Rubens and K. J. Parker, Ultrasound Med. Biol., 2008, 34, 1033–1042.
126. M. Z. Kiss, T. Varghese and T. J. Hall, Phys. Med. Biol., 2004, 49, 4207.
127. M. D. Choudhury, S. Chandra, S. Nag, S. Das and S. Tarafdar, Colloids Surf., A, 2012, 407, 64–70.
128. P. Yu, Z. Jinzhou and L. Yongming, Petrol. Explor. Dev., 2017, 44, 1038–1044.
129. P. Lidon, L. Villa and S. Manneville, Rheol. Acta, 2017, 56, 307–323.
130. J. J. Shen, C. G. Li, H. T. Wu and M. Kalantari, Korea-Aust. Rheol. J., 2013, 25, 87–93.
131. A. Holder, N. Badiei, K. Hawkins, C. Wright, P. Williams and D. Curtis, Soft Matter, 2018, 14, 574–580.
132. M. Fraldi, A. Cugno, L. Deseri, K. Dayal and N. Pugno, J. R. Soc., Interface, 2015, 12, 20150656.
133. M. Kohandel, S. Sivaloganathan, G. Tenti and K. Darvish, Phys. Med. Biol., 2005, 50, 2799.
134. R. Metzler, W. Schick, H.-G. Kilian and T. F. Nonnenmacher, J. Chem. Phys., 1995, 103, 7180–7186.
135. R. Metzler and T. F. Nonnenmacher, Int. J. Plast., 2003, 19, 941–959.
136. V. Pandey and S. Holm, Phys. Rev. E, 2016, 94, 032606.
137. J. J. Fredberg and D. Stamenovic, J. Appl. Physiol., 1989, 67, 2408–2419.
138. J. Rother, M. Büchsenschütz-Göbeler, H. Nöding, S. Steltenkamp, K. Samwer and A. Janshoff, J. R. Soc., Interface, 2015, 12, 20141057.
139. P. Cai, Y. Mizutani, M. Tsuchiya, J. M. Maloney, B. Fabry, K. J. Van Vliet and T. Okajima, Biophys. J., 2013, 105, 1093–1102.
140. B. Fabry, G. N. Maksym, J. P. Butler, M. Glogauer, D. Navajas, N. A. Taback, E. J. Millet and J. J. Fredberg, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 68, 041914.
141. B. A. Smith, B. Tolloczko, J. G. Martin and P. Grütter, Biophys. J., 2005, 88, 2994–3007.
142. P. Cai, R. Takahashi, K. Kuribayashi-Shigetomi, A. Subagyo, K. Sueoka, J. M. Maloney, K. J. Van Vliet and T. Okajima, Biophys. J., 2017, 113, 671–678.
143. C. Friedrich and L. Heymann, J. Rheol., 1988, 32, 235–241.
144. H. J. Haubold, A. M. Mathai and R. K. Saxena, J. Appl. Math., 2011, 2011, 1–51.
145. R. Gorenflo, J. Loutchko and Y. Luchko, Fract. Calc. Appl. Anal., 2002, 6(4), 491–518.
146. M. Nobile, V. Pirozzi, E. Somma, G. Gomez D'Ayala and P. Laurienzo, J. Polym. Sci., Part B: Polym. Phys., 2008, 46, 1167–1182.
147. M. Moresi and M. Bruno, J. Food Eng., 2007, 82, 298–309.
148. H. H. Winter and M. Mours, Neutron spin echo spectroscopy viscoelasticity rheology, Springer, 1997, pp. 165–234.
149. G. Odian, et al., Principles of polymerization, John Wiley & Sons, 2004.
150. R. L. Bagley, AIAA J., 1989, 27, 1412–1417.
151. F. Mainardi, Discrete Cont. Dyn. Syst., 2014, 19, 2267–2278.
152. H.-J. Lee and Y. R. Kim, J. Eng. Mech., 1998, 124, 32–40.
153. S. A. Forough, F. M. Nejad and A. Khodaii, Int. J. Pavement Eng., 2016, 17, 314–330.
154. Z. Li, Z. Zhang and C. Thomas, Innov. Food Sci. Emerg. Technol., 2016, 34, 44–50.
155. A. Shahin-Shamsabadi, A. Hashemi, M. Tahriri, F. Bastami, M. Salehi and F. M. Abbas, Mater. Sci. Eng., C, 2018, 90, 280–288.
156. J. Kaplan, A. Bonfanti and A. Kabla, J. Open Source Softw., 2019, 4, 1700.
157. R. Garra and R. Garrappa, Commun. Nonlinear Sci. Numer. Simul., 2018, 56, 314–329.
158. J. Wang, Y. Zhou and D. O'Regan, Integr. Transforms Spec. Funct., 2018, 29, 81–94.

---

Footnote

† Electronic supplementary information (ESI) available: Fractional viscoelastic models for power-law materials – Annex: summary of their constitutive equations, moduli & qualitative behaviours. See DOI: 10.1039/d0sm00354a

---