Simultaneous Characterization of Relaxation, Creep, Dissipation, and Hysteresis by Fractional-Order Constitutive Models

: We considered relaxation, creep, dissipation, and hysteresis resulting from a six-parameter fractional constitutive model and its particular cases. The storage modulus, loss modulus, and loss factor, as well as their characteristics based on the thermodynamic requirements, were investigated. It was proved that for the fractional Maxwell model, the storage modulus increases monotonically, while the loss modulus has symmetrical peaks for its curve against the logarithmic scale log ( ω ) , and for the fractional Zener model, the storage modulus monotonically increases while the loss modulus and the loss factor have symmetrical peaks for their curves against the logarithmic scale log ( ω ) . The peak values and corresponding stationary points were analytically given. The relaxation modulus and the creep compliance for the six-parameter fractional constitutive model were given in terms of the Mittag–Lefﬂer functions. Finally, the stress–strain hysteresis loops were simulated by making use of the derived creep compliance for the fractional Zener model. These results show that the fractional constitutive models could characterize the relaxation, creep, dissipation, and hysteresis phenomena of viscoelastic bodies, and fractional orders α and β could be used to model real-world physical properties well. the creep compliance for the six-parameter fractional constitutive model in terms of the Mittag–Lefﬂer functions. Then, the stress–strain hysteresis loops were simulated by making use of the derived creep compliance for the fractional Zener model. These results show that the fractional constitutive models could simulate the relaxation, creep, dissipation, and hysteresis phenomena of viscoelastic bodies, and the tuning of the fractional orders α and β could model the physical properties well.

The applications of fractional calculus to viscoelasticity were promoted by the wide use of viscoelastic materials in damping design and shock absorption. Scott-Blair [42,43] presented a fractional constitutive relation σ(t) = E (α) (t), where 0 < α < 1 and E and α are material constants for a viscoelastic material that is an intermediate between a pure elastic solid (Hooke model) and a pure viscous fluid (Newton model). In [9,44], this relation was called the Scott-Blair model. Macromolecule materials, such as butyl and polybutadiene, are typical viscoelastic bodies whose constitutive relations can be modeled by using fractional derivatives [45]. Fractional constitutive models, such as the fractional Maxwell, Kelvin-Voigt, and Zener models, were suggested by replacing the classical Newton element with the Scott-Blair element [9,23,[45][46][47][48][49][50], i.e., substituting the integerorder derivatives into the fractional-order derivatives. In [51], viscoelastic dampers were simulated by using the Maxwell and Kelvin-Voigt fractional models. In [52], a fractional integral form and differential form were compared for the the standard solid model.
In [53], a fractional five-parameter model was proposed, and it was compatible with the high-frequency metadata. In [54], a fractional six-parameter model was proposed to include all of the four types of viscoelasticity. In [55], energy storage and energy dissipation for the fractional equations were discussed. In [49,[56][57][58][59], distributed-order derivatives were introduced to model viscoelastic constitutive equations.
Suppose that f (t) is a piecewise continuous function defined on (0, +∞) and is integrable on each finite subinterval. Then, the Riemann-Liouville fractional integral of order λ is defined as for λ > 0, and is the Gamma function. The Riemann-Liouville fractional derivative of order α is defined, when it exists, as The functions that we considered are causal, viz., identically equal to zero when t < 0. We use the Laplace transform with a lower limit of integral 0 − , So, the Laplace transform of the Riemann-Liouville fractional derivative satisfies We note that the initial conditions at t = 0 + will not appear in the Laplace transform of the constitutive relation of the stress and strain on the assumption that the contributions from the initial conditions are canceled if the Laplace transform with the starting point t = 0 + is used [9].
The Mittag-Leffler functions are significant in fractional models. These types of functions have the definition [7,60,61] The following Laplace transform formula holds: Dissipation and hysteresis phenomena exist broadly in mechanical systems [9,49,[62][63][64]. Energy dissipation was considered in some fractional constitutive equations in [9,49,53,55,65]. In [36], fractional-order models were proposed to describe the broadband hysteresis of a piezoelectric actuator by looking into the relationship between the input voltage and the output displacement. Simulations and experiments were performed to validate the effectiveness of the fractional-order model.
The purpose of this work is to consider dissipation, creep, relaxation, and hysteresis for a six-parameter fractional constitutive model and its particular cases. In the next section, we consider the storage modulus, loss modulus, and loss factor, as well as their characteristics for different fractional constitutive models based on the thermodynamic requirements. For the fractional Maxwell model, we prove that the storage modulus increases monotonically, while the loss modulus has symmetrical peaks for its curve against the logarithmic scale log(ω). For the fractional Zener model, we demonstrate that the storage modulus monotonically increases, while the loss modulus and the loss factor have symmetrical peaks for their curves against the logarithmic scale log(ω). In Section 3, the relaxation, creep, and hysteresis for the six-parameter fractional constitutive model are considered by making use of the Mittag-Leffler functions.

Fractional Models and Thermodynamic Requirements
We consider the six-parameter fractional constitutive equation where the coefficients and the orders are subject to the constraints The creep compliance and relaxation modulus of a six-parameter model similar to that in Equation (7) were considered by using the complex inversion formula of the Laplace transform, and the results were expressed by infinite integrals. In this article, we focus on the thermodynamic requirements in this section and the Mittag-Leffler function expressions of the creep compliance and relaxation modulus in the next section for the six-parameter model (7). Based on these results, the storage modulus, loss modulus, and loss factor are investigated and hysteresis loops are simulated. Operating the Laplace transform for Equation (7) leads to The complex modulus is defined as The storage modulus G s and loss modulus G l are defined as the the real part and the imaginary part of the complex modulus G * (ω) = G s (ω) + iG l (ω).
The loss factor is defined as The loss factor is a measure of the damping ability of a linear viscoelastic material. We note that the loss factor is also called loss tangent [9]. The storage modulus G s and loss modulus G l are also denoted as G and G in [9,52,66].
For a constitutive equation of real viscoelastic materials, thermodynamic requirements must be satisfied. It is known from thermodynamics that the internal work and the dissipated energy must be positive. We can satisfy the thermodynamic restrictions by ensuring that both the storage modulus and the loss modulus are positive for all frequencies [65], i.e., G s (ω) ≥ 0 and G l (ω) ≥ 0 for all ω > 0.
Inserting the formula i λ = cos( πλ 2 ) + i sin( πλ 2 ) into Equation (11) and separating the real part from the imaginary part, we obtain the loss modulus and the storage modulus: where By using Equations (15) and (16), the following proposition can be directly checked. (8) and (9), Equation (7) satisfies the thermodynamic requirements.

Proposition 1. Under the assumptions in
The model (7) covers common special cases, including the classical integer-order equations. We firstly list five integer-order equations, as well as their loss and storage moduli and loss factors.
Kelvin-Voigt model (a = b 2 = 0, α = 1): Here G s increases monotonically from 0 to b 1 /a, G l has a symmetrical peak for plot with respect to the logarithmic scale log(ω).

Zener model
Here, G s increases monotonically from b 0 to b 1 /a, and G l and ζ have symmetrical peaks for their plots with respect to the logarithmic scale log(ω).
Next, we consider the fractional cases.

Fractional Scott-Blair model
The loss modulus and the storage modulus have the same power law with respect to the frequency ω, and the loss factor is independent of ω.
Fractional Kelvin-Voigt model (a = b 2 = 0): The two moduli G l and G s increase in their power law with respect to ω; ζ(ω) monotonically increases on [0, +∞) and satisfies Before considering new models, we give the following lemma.
has the unique stationary point and the peak value and f (ω) increases monotonically on [0, ω * ] and decreases monotonically on [ω * , +∞). For any positive real number p, the following equation holds: Proof. From the derivative we have the stationary point ω * , the peak value f (ω * ) in Equation (37), and the monotonicity on the intervals [0, ω * ] and [ω * , +∞). For Equation (38), first, by direct substitution, it is easy to verify that Next, we have By the relation (39), the Equation (38) is obtained.
The loss factor ζ(ω) monotonically decreases on [0, +∞) and satisfies For the storage and loss moduli, we give the following proposition.

Proposition 2.
For the fractional Maxwell model (40), the storage modulus G s (ω) increases monotonically on the interval [0, +∞) and satisfies the loss modulus G l (ω) has the unique stationary point and the peak value has symmetrical peaks for its curve against the logarithmic scale log(ω).

Proof. From the derivative
we know that the storage modulus G s (ω) increases monotonically on the interval [0, +∞). Equation (45) is direct from (42). The results for the loss modulus come from Lemma 1.

Fractional Zener model
First, we notice that when α = 0, the model becomes the Hook model, and there is no energy dissipation. For the special case ab 0 = b 1 , i.e., the coefficients satisfy the proportional relation 1/b 0 = a/b 1 , the model stores an invariant energy and has no energy dissipation like a Hook model, In the following, we always assume that 0 < α ≤ 1 and b 1 > ab 0 .

Proposition 3.
For the fractional Zener model (48), (i) the storage modulus G s (ω) monotonically increases on the interval [0, +∞) with the limits (ii) the loss modulus G l (ω) has the unique stationary point and the peak value (iii) the loss factor ζ has the unique stationary point and the peak value (iv) the graphs of the loss modulus and the loss factor against the logarithmic scale log(ω) have symmetrical peaks, the stationary points of the loss modulus and the loss factor satisfy the inequality the peak values of the loss modulus and loss factor increase monotonically with increasing α and satisfy the inequalities where either of the equal signs holds only if α = 1, and the two peak values satisfy the relation Proof. (i) It is easy to check that the derivative G s (ω) ≥ 0. This means that G s (ω) is monotonically increasing on the interval [0, +∞). The limits in Equation (53) are directly from (50). (ii) and (iii) result from Lemma 1. (iv) From Equation (38) in Lemma 1, we know that the graphs of the loss modulus and the loss factor against the logarithmic scale log(ω) have symmetrical peaks. By the hypothesis b 1 > ab 0 , we have b 0 /(ab 1 ) < 1/a 2 . This leads to the inequality ω f < ω l , referring to Equations (54) and (56). From Equations (55) and (57), the peak values G lmax and ζ max increase monotonically with increasing α and satisfy the inequalities in (59). ).
This is equivalent to by using (55). Finally, (60) is obtained by using the inequality b 1 > ab 0 .
In Figures 1 and 2, we plot the storage modulus, the loss modulus, and the loss factor against ω for two groups of different parameter values. In Figure 1, ζ max < G l max , while in Figure 2, ζ max > G l max . The effects in the order α on the storage modulus G s and the loss factor ζ are shown in Figures 3 and 4, respectively. The increase of order α makes the curves of G s ∼ ω become steeper, but α does not affect the limits G s (0) and G s (+∞). The increase in order α enhances the peaks of the curves of ζ ∼ ω. The plot of G l ∼ ω is similar to that of ζ ∼ ω. The following models involve the term of the fractional derivative of order β, and we assume that α < β.

Fractional four-parameter model
By combining their derivatives, the following assertions can be yielded: G l (ω) and G s (ω) increase monotonically on [0, +∞) with limits G l (0) = 0, G l (+∞) = +∞, (65) G s (0) = 0, G s (+∞) = +∞, (66) and ζ(ω) decreases monotonically on [0, +∞) with the limits Three monotonic curves of G s , G l , and ζ versus ω are plotted in Figure 5. Fractional five-parameter model (b 0 = 0): By examining their derivatives, we can conclude that G s (ω) increases monotonically on [0, +∞) with the limits G s (0) = 0, G s (+∞) = +∞; (72) G l (ω) satisfies G l (0) = 0 and increases monotonically to +∞ when ω is large enough; and ζ(ω) satisfies and decreases monotonically when ω is large enough. From Equation (73), the linear-logarithmic curve of ζ ∼ ω does not have symmetrical peaks if 2α = β. In Figures 6 and 7, curves of G s , G l , and ζ versus ω are plotted, where the orders α and β affect the limits of ζ as ω → 0 or +∞, but do not affect the limits of G s and G l . In order to further examine the effects of the orders α and β on the loss factor ζ, we plot the curves of ζ versus ω for different order values in Figures 8 and 9. The asymmetrical peaks of the loss factor ζ are tuned by the two orders α and β.  We noticed that the two fractional derivative terms on the right-hand side of the constitutive equation lead to an asymmetrical loss factor peak. In [53], a five-parameter fractional derivative model, which also includes two fractional derivative terms with respect to the strain, but different from those of (68), was presented to describe this phenomenon.

Relaxation, Creep, and Hysteresis
The relaxation modulus and creep compliance are two important material functions describing the mechanical characteristics of viscoelastic materials. The relaxation modulus R(t) is the behavior of stress relaxing with time under a suddenly applied unit strain (t) = H(t), where H(t) is the Heaviside unit-step function. Creep compliance C(t) is the strain response to an instantaneously applied stress σ(t) = H(t). Both the material functions are nonnegative. Furthermore, for 0 < t < +∞, R(t) is decreasing and C(t) is increasing.
Inserting¯ (s) = 1/s into Equation (10), we obtain the relaxation modulus in the Laplace domain:R If a = 0, i.e., the model (7) degenerates into the following five-parameter model, then the relaxation modulus is If a > 0, we rewrite Equation (76) as By Formula (6), the inverse Laplace transform yields the relaxation modulus in terms of the Mittag-Leffler functions The relaxation moduli R(t) in Equations (78) and (80) are plotted in Figures 12 and 13, respectively. The appearance of the fractional derivative term with respect to stress accelerates the relaxation of stress by comparing the two figures. In addition, the variation in the order tunes the relaxing rate.  In a similar manner, lettingσ(s) = 1/s in Equation (10), we obtain the creep compliance in the Laplace domain:C We give the creep compliance C(t) by discriminating the following three cases.
In this case, it is necessary that b 0 > 0 and a = 0 from the condition (8). The model degenerates into a Hook spring σ(t) = b 0 (t), and the creep compliance is In fact, there is no "creep".
Case II. b 2 = 0 and b 1 = 0 This case corresponds to the fractional Zener model (48). Rewriting Equation (81) as and applying the inverse Laplace transform using Formula (6), we derive the creep compliance in terms of the Mittag-Leffler functions: In Figure 14, the curves of the creep compliance C(t) in Equation (84) are plotted for different values of α. Apart from an initial period, the creep compliance C(t) increases with the rising α. Rewriting Equation (81) as and operating the inverse Laplace transform using Formula (6), we have the creep compliance: Here, the creep compliance is expressed as a series in which the Mittag-Leffler functions and their derivatives are involved.
In terms of the creep compliance C(t) and the relaxation modulus R(t), the strain and the stress have the expressions respectively, where * denotes the convolution.
Making use of the creep compliance C(t) and Formula (87), we may simulate the hysteresis phenomena. We consider the strain response to the sinusoidal stress loading σ(t) = sin(t). From (87), the strain response is For the numerical simulation here, we use the fractional Zener model and its creep compliance in Equation (84). Taking a = 0.5, b 0 = 1, b 1 = 3, and four different values of α (0.25, 0.5, 0.75, and 1), the stress-strain hysteresis loops are depicted in Figure 15, where t serves as the parameter varying from 0 to 8π. For comparison, the coordinate ranges are the same in the four subfigures. The unloading curve falls to the right side of the loading curve to form the hysteresis loops. In other words, the fluctuations of the strain fall behind those of the stress. With increasing α, the hysteresis loops grow fat and encircle a larger area, which means that the ability of dissipation of the model is raised.

Conclusions
We considered dissipation, creep, relaxation, and hysteresis resulting from a sixparameter fractional constitutive model and its particular cases. In Section 2, we considered the storage modulus, loss modulus, and loss factor, as well as their characteristics for different fractional constitutive models based on the thermodynamic requirements. We proved that for the fractional Maxwell model, the storage modulus G s increases monotonically, while the loss modulus G l has symmetrical peaks for its curve against the logarithmic scale log(ω), and for the fractional Zener model, the storage modulus G s monotonically increases, while the loss modulus G l and the loss factor ζ have symmetrical peaks for their curves against the logarithmic scale log(ω). The peak values and corresponding stationary points were accurately given. A series of results for the fractional Zener model were presented in Proposition 3. In Section 3, we derived the relaxation modulus and the creep compliance for the six-parameter fractional constitutive model in terms of the Mittag-Leffler functions. Then, the stress-strain hysteresis loops were simulated by making use of the derived creep compliance for the fractional Zener model. These results show that the fractional constitutive models could simulate the relaxation, creep, dissipation, and hysteresis phenomena of viscoelastic bodies, and the tuning of the fractional orders α and β could model the physical properties well.