Exact Mechanical Hierarchy of Non-Linear Fractional-Order Hereditariness

Abstract

Non-local time evolution of material stress/strain is often referred to as material hereditariness. In this paper, the widely used non-linear approach to single integral time non-local mechanics named quasi-linear approach is proposed in the context of fractional differential calculus. The non-linear model of the springpot is defined in terms of a single integral with separable kernel endowed with a non-linear transform of the state variable that allows for the use of Boltzmann superposition. The model represents a self-similar hierarchy that allows for a time-invariance as the result of the application of the conservation laws at any resolution scale. It is shown that the non-linear springpot possess an equivalent mechanical hierarchy in terms of a functionally-graded elastic column resting on viscous dashpots with power-law decay of the material properties. Some numerical applications are reported to show the capabilities of the proposed model.

1. Introduction

Time-dependent behavior of several complex materials is a crucial issue for a reliable mechanical description as well as virtual simulation. Indeed, a recent approach to the mechanics of several materials with “in-silico” simulations takes into account the non-local time behavior of materials in terms of material’s memory with respect to past strain/stress histories up to current time. This feature of materials, very often encountered in classical engineering materials such as concrete and wood, as well as in more advanced materials such as biomaterials and biological tissues, is called material hereditariness.

Material hereditariness has been initially formulated in terms of linear models that, taking full advantage of Boltzmann superposition, yield a time-dependent evolution ruled by convolution integrals involving the creep function $J\left(t\right)$ or relaxation function $G\left(t\right)$ [1,2].

In more details, creep is known as the the time-dependent strain evolution under constant stress; the relaxation is, instead, the stress decay under uniform strain [1,2].

The main assumption related to use of Boltzmann superposition is the linear dependence of creep (relaxation) on the real level of the strain (stress). Indeed, the creep function $J\left(t\right)$ is defined as the strain evolution $\epsilon \left(t\right)$ under unitary stress $\sigma \left(t\right)=1$, whereas stress relaxation $\sigma \left(t\right)$ is the evolution under constant unitary strain $\epsilon \left(t\right)=1$. In such circumstances, a widely used mathematical representation of creep and relaxation of polymers and rubbers is the power-law with real exponent as $G\left(t\right)={\displaystyle \frac{G_0}{\mathsf{\Gamma }(1-\beta )}}{\left({\displaystyle \frac{t}{\tau }}\right)}^{-\beta }$ and $J\left(t\right)={\displaystyle \frac{1}{G_0\mathsf{\Gamma }(1+\beta )}}{\left({\displaystyle \frac{t}{\tau }}\right)}^{\beta }$ with $\beta \in [0,1]$, $\mathsf{\Gamma }(·)$ the Euler-Gamma function, $\left[G_0\right]=F/L^2$ the elastic modulus of the material, and $\left[{\tau }_0\right]=T$ a material characteristic time [3,4].

The use of power-laws for creep and relaxation allow introducing constitutive stress–strain relations in terms of time non-local convolution integrals with power-law kernels yielding constitutive equations in terms of fractional differential calculus [5,6,7,8]. Among them, mechanical hierarchy has been introduced to capture the power-law evolution [9,10,11] and provide a mechanical description of real-order integro-differential operators [8]. Fractional-order operators are nowadays recognized as an interesting mathematical tool to model non-local problems in spatial coordinates [12,13,14,15,16,17] as well as in time parameter [18,19,20,21]. Fractional-order operators involve linear convolution of the state variables with power-law kernels and few contributions have been provided to deal with non-linear problems, mainly due to non-linear measures of stress and strains [22]. Time non-local formulations in presence of material non-linearity have been faced in the scientific literature since the middle of the last century based on the principle of fading memory allowing, after some manipulations, to express the constitutive equations as a sum of multiple integrals involving several material functions [23]. Neglecting multiple integrals with respect to the single integral term leads to the widely used Quasi-Linear Viscoelasticity (QLV) [24,25,26]. A comprehensive review of linear and non-linear material hereditariness has been reported in some studies [27,28,29,30,31,32,33,34]. QLV has not, however, been used in the context of fractional-order calculus to take full advantages of fractional-order formalism. Indeed, the formulation of QLV is usually presented in terms of material relaxation that is not readily obtained from experimental data for the inertia of loading equipment. Creep functions are, instead, more easily obtained but no closed-form relations among creep and relaxations for the QLV have been reported so far in the scientific literature to the best of the authors’ knowledge. Moreover, no mechanical justification for the use of QLV models have been introduced, providing a severe limitation in the use of the quasi-linear formulations.

Fractional-order QLV theory for measures and estimates of material parameters have been recently reported in the scientific literature for tendons of the knee with a specific experimental protocol [35,36].

Non-linear models involving fractional-order operators have been discussed in the context of non-linear dynamics [37] as well as in the context of large vibrations of plane structures with fractional-order damping [38].

In this paper, the authors aim to show that the Nutting non-linear relation [39] corresponds to a Fractional-order description of QLV, dubbed Quasi-Fractional Hereditary Material (Q-FHM) models. The main advantages of such approach relies in a closed-form relations among creep and relaxation parameters. The Q-FHM mechanical element is shown to be exactly equivalent to a mechanical hierarchy completely analogous to the well-established hierarchy for fractional hereditary materials. The proposed mechanical model constitutes a self-similar hierarchy that, as a result of the proper scaling of material coefficients and fulfilling conservation of linear momentum, yields a time-invariance at any resolution scale. Some numerical applications are reported in the paper to show the equivalence among the power-laws modeling and the results of the application of the proposed non-linear hierarchy.

2. Linear (LV) and Quasi-Linear (QLV) Viscoelasticity of Materials

Material viscoelasticity (often referred hereditariness) is experienced each time a long-standing controlled load (controlled displacement) experimental test shows time evolution of the measured displacements (measured load). In this section, we assume a 1D load–displacement relation and we switch to engineering measure of stress, namely $\sigma \left(t\right)$, and of strain, namely $\epsilon \left(t\right)$, without loss of generality. Under these circumstances, the constitutive behavior involves material function for strain evolution, namely ${\varphi }_c\left(\sigma ,t\right)$, which provides the strain evolution under constant stress, as well as a different material function yielding the stress decay under constant strain, namely ${\varphi }_r\left(\epsilon ,t\right)$. In passing, we observe that the material functions ${\varphi }_c$ and ${\varphi }_r$ depend, in general, on the applied stress and strain, respectively.

In the framework of linear hereditariness, the creep function satisfies the linearity conditions:

$${\varphi }_c\left(\lambda \sigma ,t\right)=\lambda {\varphi }_c\left(\sigma ,t\right)\forall \lambda \in R;{\varphi }_c\left({\sigma }_1+{\sigma }_2,t\right)={\varphi }_c\left({\sigma }_1,t\right)+{\varphi }_c\left({\sigma }_2,t\right)$$

(1)

A similar consideration holds true for the relaxation function that satisfies the linearity conditions in Equation (1):

$${\varphi }_r\left(\lambda \epsilon ,t\right)=\lambda {\varphi }_r\left(\epsilon ,t\right)\forall \lambda \in R;{\varphi }_r\left({\epsilon }_1+{\epsilon }_2,t\right)={\varphi }_r\left({\epsilon }_1,t\right)+{\varphi }_r\left({\epsilon }_2,t\right)$$

(2)

The linearity assumptions for the creep and relaxation functions allow introducing material hereditariness for unitary value of applied stress and strain, namely $\sigma =1$ and $\epsilon =1$, resulting in stress and strain independent material hereditary functions as ${\varphi }_c\left(1,t\right)=\sigma J\left(t\right)=1J\left(t\right)$ and ${\varphi }_r\left(1,t\right)=\epsilon G\left(t\right)=1G\left(t\right)$, respectively. Time-varying functions $\left[G\left(t\right)\right]=F/L^2$ and $\left[J\left(t\right)\right]=L^2/F$ are the well-known relaxation and creep functions, respectively.

In the following, the linearity conditions are extensively used to introduce the linear mathematical description of material hereditariness as well as to provide a rheological description of the experimental linear behavior observed for several conventional materials.

2.1. Mathematical Modeling of LV

The knowledge of the material functions $J\left(t\right)$ and $G\left(t\right)$, creep and relaxation functions, respectively, allows for the use of Boltzmann superposition principle [40], yielding the stress and strain at a generic time instant t due to an arbitrary stress $\sigma \left(\tau \right)$ or strain $\epsilon \left(\tau \right)$ history as:

$$\begin{array}{c}\hfill \epsilon \left(t\right)={\int }_0^{t}J\left(t-\tau \right)\dot{\sigma }\left(\tau \right)d\tau \end{array}$$

(3a)

$$\begin{array}{c}\hfill \sigma \left(t\right)={\int }_0^{t}G\left(t-\tau \right)\dot{\epsilon }\left(\tau \right)d\tau \end{array}$$

(3b)

The equalities in Equation (3) are the well-known integral constitutive relations in 1D linear hereditariness and it is well-known, after some straightforward manipulations, that creep and relaxation functions must satisfy the fundamental equation of linear hereditariness $\widehat{G}\left(s\right)\widehat{J}\left(s\right)={\displaystyle \frac{1}{s^2}}$, with $\left[s\right]={\displaystyle \frac{1}{T}}$ the Laplace parameter and $\widehat{\left[•\right]}$ the Laplace transform.

The specific functional class of creep and relaxation functions reported in Equation (3) may be guessed from experimental data collected in the course of experimental campaigns and they are very often expressed as single or linear combinations of exponential functions by means of Prony representation theorem [41] as:

$$\begin{array}{cc}\hfill & J\left(t\right)=\sum _{r=1}^M{J}_r\left(1-exp\left(-\frac{t}{{\tau }_j^{\left(c\right)}}\right)\right)\hfill \end{array}$$

(4a)

$$\begin{array}{cc}\hfill & G\left(t\right)=\sum _{r=1}^N{G}_{r}exp\left(-{\displaystyle \frac{t}{{\tau }_j^{\left(r\right)}}}\right)\hfill \end{array}$$

(4b)

where the coefficients of the expansions have physical measures $\left[G_r\right]={\displaystyle \frac{F}{L^2}}$ and $\left[J_r\right]={\displaystyle \frac{L^2}{F}}$ and the material characteristic times in creep and relaxations, namely ${\tau }_j^{\left(c\right)}$ and ${\tau }_j^{\left(r\right)}$, are additional material parameters that may be estimated by best fitting procedures together with the expansion coefficients. The integer numbers in the expansions, namely M and N, are the order of the Prony series used for creep and relaxation, respectively.

The expressions for creep and relaxation functions reported in Equation (4) cannot, however, satisfy the fundamental relation of linear hereditariness, and, henceforth, they must be used separately in stress- and strain-based constitutive relations reported in Equation (3). Some attempts to introduce analogous formulations joint in creep and relaxation led to unphysical negative values of the material relaxation times in the Prony expansion [41].

The functional stress–strain relations reported in Equation (4) possess an equivalent differential formulation in terms of elastic (Hookean) and viscous (Newtonian) elements. In more detail, we report in Figure 1a the differential formulation, named rheological representation of the Prony series expansion of the creep function $J\left(t\right)$. Similarly, the mechanical arrangements springs and dashpots reported in Figure 1b correspond to the rheological representation of the relaxation function $G\left(t\right)$ reported in Equation (4b).

Direct comparisons of Figure 1a,b show that the mechanics beyond the creep and relaxation functions described by Prony series expansion is quite different, as shown by the series and parallel arrangements of springs and dashpots that correspond to the prescribed analytical expression in Equation (3). Such a consideration is a direct consequence on the lack of mathematical consistency of creep and relaxation functions expressed in terms of Prony series expansions.

In passing, we observe that, as long as $N=M=1$, the well-known Maxwell elements representing relaxation and Kelvin–Voigt element for creep are obtained. Indeed, as $M=1$, the governing equation of the rheological model for the relaxation function reads:

$$\begin{array}{c}{\tau }_r\dot{G}\left(t\right)+G\left(t\right)={\eta }_1\dot{\epsilon }\left(t\right)={\eta }_1\delta \left(t\right)\hfill \end{array}$$

(5a)

$$\begin{array}{c}G\left(0\right)=E_1\hfill \end{array}$$

(5b)

where the relaxation time is defined as ${\tau }_r={\eta }_1/E_1$; $E_1$ and ${\eta }_1$ are the elastic modulus of the Hookean and Newtonian elements, respectively; and $\delta \left(•\right)$ is the Dirac’s delta function.

Solving the differential problem in Equation (5) yields the non-dimensional relaxation $\overline{G}\left(t\right)$ in the form:

$$G\left(t\right)=E_{1}exp\left(-{\displaystyle \frac{t}{{\tau }_r}}\right)$$

(6)

that correspond to the exponential form often used in linear hereditariness.

In a similar fashion, as $N=1$, the rheological model for the creep function $J\left(t\right)$ is ruled by the differential equation:

$$\begin{array}{c}{\tau }_c\dot{J}\left(t\right)+J\left(t\right)={\displaystyle \frac{\sigma \left(t\right)}{E_1}}={\displaystyle \frac{H\left(t\right)}{E_1}}\hfill \end{array}$$

(7a)

$$\begin{array}{c}J\left(0\right)={\displaystyle \frac{1}{E_1}}\hfill \end{array}$$

(7b)

with ${\tau }_c={\tau }_r={\displaystyle \frac{{\eta }_1}{E_1}}$, $H\left(t\right)$ the Heaveside function yielding the solution as:

$$J\left(t\right)=\left[1-exp\left(-{\displaystyle \frac{t}{{\tau }_c}}\right)\right]$$

(8)

2.2. Mathematical Modeling of QLV

Linear model introduced in previous studies does not account for the non-proportional evolution of creep as well as of relaxation often observed in polymeric and rubberlike materials.

This aspect, often observed at the beginning of the 1990s, has been modeled by means of a multiplicative decomposition of the material functions as:

$${\varphi }_r\left(\epsilon ,t\right)=s_n\left(t\right)G\left(t\right);{\varphi }_c\left(\sigma ,t\right)=e_n\left(\sigma \right)J\left(t\right)$$

(9)

with $\left[s_n\left(\epsilon \right)\right]=\left[e_n\left(\sigma \right)\right]$ non-dimensional, non-linear functions of strain and stress level at time instant t, respectively. Functions $G\left(t\right)$ and $J\left(t\right)$ are, instead, relaxation and creep functions observed in linear viscoelasticity reported in Equations (4a) and (4b). Under this condition, the assumption of effect superposition yields:

$$\begin{array}{c}\epsilon \left(t\right)={\int }_0^{t}J\left(t-\tau \right)d{s}_n\left(\tau \right)={\int }_0^{t}J\left(t-\tau \right){\dot{s}}_n\left(\tau \right)d\tau \hfill \end{array}$$

(10a)

$$\begin{array}{c}\sigma \left(t\right)={\int }_0^{t}G\left(t-\tau \right)d{e}_n\left(\tau \right)={\int }_0^{t}G\left(t-\tau \right){\dot{e}}_n\left(\tau \right)d\tau \hfill \end{array}$$

(10b)

The choice for $s_n\left(\epsilon \right)=\partial {\psi }_H/\partial \epsilon$ and $e_n\left(\sigma \right)=\partial {\psi }_G/\partial \sigma$, with ${\psi }_H\left(\epsilon \right)$ and ${\psi }_G\left(s_n\right)$ being the Helmoltz and the Gibbs free energy function of hyperelastic materials, yields the QL model of non-linear hereditariness as:

$$\begin{array}{c}\epsilon \left(t\right)={\int }_0^{t}J\left(t-\tau \right)\dot{\left({\displaystyle \frac{\partial {\psi }_H}{\partial \epsilon }}\right)}d\tau ={\int }_0^{t}J\left(t-\tau \right){\displaystyle \frac{{\partial }^2{\psi }_H\left(\epsilon \right)}{\partial {\epsilon }^2}}\dot{\epsilon }\left(\tau \right)d\tau \hfill \end{array}$$

(11a)

$$\begin{array}{c}\sigma \left(t\right)={\int }_0^{t}G\left(t-\tau \right){\displaystyle \frac{\partial {\psi }_G}{\partial \sigma }}\dot{\sigma }d\tau ={\int }_0^{t}G\left(t-\tau \right){\displaystyle \frac{{\partial }^2{\psi }_G\left(\sigma \right)}{\partial {\sigma }^2}}\dot{\sigma }\left(\tau \right)d\tau \hfill \end{array}$$

(11b)

The conventionally used functional class of time-varying functions $J\left(•\right)$ and $G\left(•\right)$ corresponds to the exponential combination expressed in Equations (4a) and (4b), respectively, for relaxation and creep. However, the mechanical models corresponding to Prony series expansion does not allow representing a rheological model similar to the linear case and, in this scenario, any attempt appears to be fundamentally flawed.

Moreover, the use of QL models do not allow for an estimate of the creep function $J\left(t\right)$ for measured relaxation and vice versa. Indeed, as long as the separable form of the material relaxation function is factorized in separable as ${\varphi }_r\left(\epsilon ,t\right)={\varphi }_e\left(\epsilon \right)G\left(t\right)$, the corresponding creep material function is ${\varphi }_c\left(\sigma ,t\right)\ne {\varphi }_e\left(\sigma \right)J\left(t\right)$ [42].

In the following section, we show that a different scenario is involved as we consider a power-law for creep and relaxation evolution as well-as a power-law for the stress (strain) dependence of the creep (relaxation) ${\varphi }_c$ (${\varphi }_r$) material functions.

3. Fractional Hereditary Materials (FHM) and Quasi-Fractional Hereditary Materials (Q-FHM)

In this section, we aim to show that power-laws expressed in terms of real powers of time for the material creep $t^{{\beta }_c}$ and relaxation functions $t^{-{\beta }_r}$ may be used for a consistent mathematical description of linear and non-linear models of material hereditariness.

Fractional-order operators have been shown to be very useful for describing the mechanical behavior of several engineering materials such as concrete [9,43], composites, polymers, and rubbers [2,44,45] under some restrictions.

Fractional-order calculus has also been applied in other fields of applied mechanics, such as heat transfer modeling [46], diffusive flow [47,48], wave propagation [49], and non-local elasticity [50,51,52,53]. For a comprehensive review, the readers may refer to [43]. Some stability mechanics problems involving non-conventional description of material external restraints have also been represented by means of fractional calculus [54,55].

3.1. Mathematical Modeling of FHM

In this section, we aim to provide a mathematical framework for the power-law time dependence of creep and relaxation functions. In such circumstances, the creep $J\left(t\right)$ and relaxation $G\left(t\right)$ are expressed as:

$$\begin{array}{c}J\left(t\right)={\displaystyle \frac{J_{\beta }}{\mathsf{\Gamma }\left(1+\beta \right)}}{t}^{\beta }={\displaystyle \frac{1}{G_{\beta }\mathsf{\Gamma }\left(1+\beta \right)}}{t}^{\beta }={\displaystyle \frac{1}{G_0\mathsf{\Gamma }\left(1+\beta \right)}}{\left({\displaystyle \frac{t}{{\tau }_0}}\right)}^{\beta }\hfill \end{array}$$

(12a)

$$\begin{array}{c}G\left(t\right)={\displaystyle \frac{G_{\beta }}{\mathsf{\Gamma }\left(1-\beta \right)}}{t}^{-\beta }={\displaystyle \frac{G_0}{\mathsf{\Gamma }\left(1-\beta \right)}}{\left({\displaystyle \frac{t}{{\tau }_0}}\right)}^{-\beta }\hfill \end{array}$$

(12b)

where the anomalous terms $\left[G_{\beta }\right]={\displaystyle \frac{F{T}^{\beta }}{L^2}}\ge 0$, $\left[J_{\beta }\right]={\displaystyle \frac{L^2}{F{T}^{\beta }}}\ge 0$ are material dependent coefficients, the exponent $0\le \beta \le 1$ is a material dependent decaying order, and $\mathsf{\Gamma }\left(•\right)$ is the Euler-Gamma function. The physical dimensions of the material coefficients allows representing the creep and relaxation functions introducing a two-term factorization $G_{\beta }=G_0{\tau }_0^{\beta }$ with ${\tau }_0$ a material dependent characteristic time and $G_0$ the conventional elastic modulus of the material observed in monotone tensile test.

Straightforward manipulations show that Equation (12) satisfies the fundamental relations of linear hereditariness and, by substituting it into Equation (2), the stress–strain constitutive equations of linear hereditariness read:

$$\begin{array}{c}\sigma \left(t\right)={\displaystyle \frac{G_0{\left({\tau }_0\right)}^{\beta }}{\mathsf{\Gamma }\left(1-\beta \right)}}{\int }_0^t{\left(t-\tau \right)}^{-\beta }\dot{\epsilon }\left(\tau \right)d\tau =G_0{\left({\tau }_0\right)}^{\beta }\left(D_0^{\beta }\epsilon \right)\left(t\right)\hfill \end{array}$$

(13a)

$$\begin{array}{c}\epsilon \left(t\right)={\displaystyle \frac{1}{G_0{\left({\tau }_0\right)}^{\beta }\mathsf{\Gamma }\left(\beta \right)}}{\int }_0^t{\left(t-\tau \right)}^{\beta -1}\sigma \left(\tau \right)d\tau ={\displaystyle \frac{1}{G_0{\left({\tau }_0\right)}^{\beta }}}\left(I_0^{\beta }\sigma \right)\left(t\right)\hfill \end{array}$$

(13b)

Notations introduced in the last equality at the right hand side of Equation (13), namely $\left(D_0^{\beta }f\right)\left(t\right)$ and $\left(I_0^{\beta }f\right)\left(t\right)$ denote, respectively, the Caputo fractional derivative and the Riemann–Liouville fractional integral of order $\beta$ of the generic function $f\left(t\right)$. Details about fractional-order operators are out of the aims of the paper and the reader may refer to more complete readings about the topic (e.g., [1,5]).

Close observation of Equation (12) reveals that the insurgence of fractional-order operators of Caputo-type is due to the specific choice of the time-decaying functions in terms of power-law types as $t^{\beta }$ or $t^{-\beta }$, respectively. This choice probably depends on the specific functional class of experimental data to be fit by the creep or relaxation functions. If alternative functional classes of the material functions are better suited to describe experimental data, e.g., a class of non-singular kernels as exponential-type functions, then alternative fractional-order operators are obtained in Equation (13) involving the so-called Caputo–Fabrizio fractional differential operator [56,57].

In the following, for simplicity, we refer to non-dimensional stress $s\left(t\right)=\sigma \left(t\right)/G_0$, yielding the constitutive equations for fractional-order non-linear hereditariness as:

$$\begin{array}{c}s\left(t\right)={\left({\tau }_0\right)}^{\beta }\left(D_0^{\beta }\epsilon \right)\left(t\right)\hfill \end{array}$$

(14a)

$$\begin{array}{c}\epsilon \left(t\right)={\left({\tau }_0\right)}^{-\beta }\left(I_0^{\beta }s\right)\left(t\right)\hfill \end{array}$$

(14b)

The linear model of fractional-order material hereditariness possesses a rheological element with intermediate behavior among spring and dashpot, dubbed springpot after [6], and it is modeled as a two-parameter element, namely the characteristic time ${\tau }_0$ and the order of the power-law $\beta$, as depicted in Figure 2.

3.2. Mathematical Modeling of Q-FHM

The presence of material non-linearity has often been observed in creep and relaxation material tests as coefficients of the power-law, i.e. the material characteristic time ${\tau }_0$ shows a strong dependence on the level of the applied stress or strain [36]. Under such circumstances, data analysis shows that material functions for constant stress, namely ${\varphi }_c\left(\sigma ,t\right)$, and for uniform strain, namely ${\varphi }_r\left(\epsilon ,t\right)$, may be expressed in a generalized, separable form as:

$$\begin{array}{cc}\hfill & {\varphi }_c\left(\sigma ,t\right)=J_e\left(\sigma \right)J\left(t\right)={\displaystyle \frac{{\sigma }_n}{G_0\mathsf{\Gamma }\left(1-{\beta }_c\right)}}{\left({\displaystyle \frac{t}{{\tau }_0}}\right)}^{{\beta }_c}\hfill \end{array}$$

(15a)

$$\begin{array}{cc}\hfill & {\varphi }_r\left(\epsilon ,t\right)=G_e\left(\epsilon \right)G\left(t\right)={\displaystyle \frac{G_0{\epsilon }_n}{\mathsf{\Gamma }\left({\beta }_r\right)}}{\left({\displaystyle \frac{t}{{\tau }_0}}\right)}^{-{\beta }_r}\hfill \end{array}$$

(15b)

where the specific functional class of $G_e\left(•\right)$ and $J_e\left(•\right)$ as well as the time-dependence expressed by the functions $G\left(t\right)$ and $J\left(t\right)$ may be estimated from experimental data, and they may be represented with power-law of time with averaged order ${\overline{\beta }}_c$ and ${\overline{\beta }}_r$, respectively, for creep and relaxation. Close observation of Equations (15a) and (15b) reveals that the separable form of the material function is assumed as the basis of QLV [26] where the non-linear auxiliary state variables are expressed in terms of the Helmoltz free energy $\mathsf{\Psi }\left(\epsilon \right)$ assumed to represent the elastic behavior in monotone tensile tests as ${\sigma }_n=G_e\left(\epsilon \right)={\displaystyle \frac{d\mathsf{\Psi }\left(\epsilon \right)}{d\epsilon }}$ and ${\epsilon }_n=J_e\left(\sigma \right)={\displaystyle \frac{d{\mathsf{\Psi }}^{(-1)}\left(\sigma \right)}{d\sigma }}$.

The use of the auxiliary variables allows introducing Boltzmann superposition principle, yielding, for non-dimensional variables $s\left(t\right)={\left(G_0\right)}^{-1}\sigma \left(t\right)$ and $s_n\left(t\right)={\left(G_0\right)}^{-1}{\sigma }_n\left(t\right)$:

$$\begin{array}{c}s\left(t\right)={\displaystyle \frac{{\left({\tau }_r\right)}^{{\beta }_r}}{\mathsf{\Gamma }\left(1+{\beta }_r\right)}}{\int }_0^t{\left(t-\tau \right)}^{-{\beta }_r}{\dot{\epsilon }}_n\left(\tau \right)d\tau ={\left({\tau }_r\right)}^{{\beta }_r}\left(D_0^{{\overline{\beta }}_r}{\epsilon }_n\right)\left(t\right)\hfill \end{array}$$

(16a)

$$\begin{array}{c}\epsilon \left(t\right)={\displaystyle \frac{{\left({\tau }_c\right)}^{-{\overline{\beta }}_c}}{\mathsf{\Gamma }\left(1-{\overline{\beta }}_c\right)}}{\int }_0^t{\left(t-\tau \right)}^{{\overline{\beta }}_c-1}{\dot{s}}_n\left(\tau \right)d\tau ={\left({\tau }_c\right)}^{-{\overline{\beta }}_c}\left(I_0^{{\overline{\beta }}_c}{s}_n\right)\left(t\right)\hfill \end{array}$$

(16b)

In the following, we assume that the non-linear functional class of the auxiliary variables are assumed as power-laws ${\epsilon }_n\left(t\right)={\left(\epsilon \right)}^{{\alpha }_r}$ and $s_n\left(t\right)={\left(s\right)}^{{\alpha }_c}$, yielding with the formalism of fractional-order calculus:

$$\begin{array}{c}s\left(t\right)={\left({\tau }_r\right)}^{{\beta }_r}\left(D_{0^{+}}^{{\beta }_r}\left[\epsilon {\left(t\right)}^{{\alpha }_r}\right]\right)\left(t\right)\hfill \end{array}$$

(17a)

$$\begin{array}{c}\epsilon \left(t\right)={\left[{\displaystyle \frac{1}{{\tau }_r^{{\beta }_r}}}\left(I_{0^{+}}^{{\beta }_r}s\right)\left(t\right)\right]}^{1/{\alpha }_r}\hfill \end{array}$$

(17b)

or involving the knowledge of the creep functions:

$$\begin{array}{c}\epsilon \left(t\right)={\left({\tau }_c\right)}^{-{\beta }_c}\left(I_{0^{+}}^{{\beta }_c}\left[s{\left(t\right)}^{{\alpha }_c}\right]\right)\left(t\right)\hfill \end{array}$$

(18a)

$$\begin{array}{c}s\left(t\right)={\left[{\left({\tau }_c\right)}^{{\beta }_c}\left(D_{0^{+}}^{{\beta }_c}\epsilon \right)\left(t\right)\right]}^{1/{\alpha }_c}\hfill \end{array}$$

(18b)

Observation of Equations (16a) and (16b) shows three main features:

(i) The constitutive equations for the non-dimensional stress involves a non-linear transform of the strain $\epsilon \left(t\right)\to {\epsilon }_n\left(t\right)$, a relaxation time ${\tau }_r$, and the time-decay order ${\beta }_r$.

(ii) The constitutive equation for the strain evolution involves the non-linear transform of the stress $s\left(t\right)\to s_n\left(t\right)$, a creep characteristic time ${\tau }_c$, and an evolution order ${\beta }_c$, with ${\beta }_c\ne {\beta }_r$.

(iii) Some specific relations among creep $\left({\alpha }_c,{\beta }_c,{\tau }_c\right)$ and relaxation $\left({\alpha }_r,{\beta }_r,{\tau }_r\right)$ parameters have been recently obtained as ${\alpha }_c=1/{\alpha }_r$, ${\beta }_c={\alpha }_c{\beta }_r$, ${\beta }_r={\alpha }_r{\beta }_c$, and

$${\tau }_r={\tau }_c\left[\mathsf{\Gamma }{\left({\beta }_r\right)}^{{\displaystyle \frac{1}{{\alpha }_r{\beta }_c}}}\mathsf{\Gamma }{\left({\beta }_c+1\right)}^{{\displaystyle \frac{1}{{\beta }_c}}}\right]$$

(19)

(iv) The proposed formulation for the constitutive models in terms of creep and relaxation functions may be framed in the context of thermodynamics of linear viscoelasticity [58,59], introducing the non-linear mapping in Equation (15) for the auxiliary state variables $s_n\left(t\right)$ and ${\epsilon }_n\left(t\right)$. In this context, the existent formulations for the free energy and the mechanical dissipation [60] associated with auxiliary strain and stress histories may be used for the formulation of the stress–strain relations.

It must be remarked that a similar formulation holds true also with other kinds of fractional-order operators with non-singular kernels that have recently been reported in the scientific literature [56].

In the next section, we aim to show that the quasi-linear model of fractional-order material hereditariness is exactly modeled by a rheological assembly of internal linear springs and internal linear dashpots depending on the level of the external agency.

4. Exact Mechanical Description of Fractional-Order Quasi-Linear Hereditariness

In this section, we aim to show that the proposed approach to fractional model of QLV possess an equivalent rheological model totally analogous to a mechanical hierarchy that has been proposed for the linear fractional-order viscoelasticity.

4.1. The Rheological Model of Fractional-Order Quasi-Linear Hereditary Materials (FQHM)

Rheological description of the non-linear dependence of the power-law in previous section may be captured as we introduce an unbounded linearly elastic bar that is externally restrained by a bed of viscous element along the column length, as shown in Figure 3.

The equilibrium of the generic element reads:

$$N_e(z+\Delta z)-N_e\left(z\right)=-N_v$$

(20)

where $N_v$ is the overall force of the dashpots along the column and it reads:

$$N_v\left(z\right)=-C\left(z,F_0\right)\Delta z\frac{\partial u}{\partial t}\left(z,t\right)$$

(21)

and the elastic axial stress along the column reads:

$$\begin{array}{c}{N}_e\left(z\right)=K(z,F_0)\left[u\left(z,t\right)-u\left(z-\Delta z,t\right)\right]\hfill \end{array}$$

(22a)

$$\begin{array}{c}{N}_e\left(z+\Delta z\right)=K(z+\Delta z,F_0)\left[u\left(z+\Delta z,t\right)-u\left(z,t\right)\right]\hfill \end{array}$$

(22b)

and with $u(z,t)$ the axial displacement along the column. Substitutions of Equations (21), (22a) and (22b) into Equation (20) yields the equilibrium equation in terms of the kinematic field $u(z,t)$ along the column:

$$K(z+\Delta z,F_0)\Delta u\left(z+\Delta z,t\right)-K(z,F_0)\Delta u\left(z,t\right)=C\left(z,F_0\right)\Delta z\frac{\partial u}{\partial t}\left(z,t\right)$$

(23)

which, after some algebraic manipulation, may be written as:

$$\frac{\Delta }{\Delta z}\left[K\left(z,F_0\right)\Delta u\right]=C\left(z,F_0\right)\frac{\partial u}{\partial t}$$

(24)

and introducing the relation $K\left(z,F_0\right)={\displaystyle \frac{K_a(z,F_0)}{\Delta z}}$ with $\left[K_a\right]=F$ and letting $\Delta z\to 0$, Equation (24) reads:

$$\frac{\partial }{\partial z}\left[K_a\left(z,F_0\right)\frac{\partial u}{\partial z}\right]=C\left(z,F_0\right)\frac{\partial u}{\partial t}$$

(25)

with the relevant boundary conditions:

$$\underset{z\to \infty }{lim}u\left(z,t\right)=0;\underset{z\to 0}{lim}{K}_a\left(z,t\right)\frac{\partial u}{\partial z}=F_0\forall t$$

(26)

and for the initial condition $u(z,0)=0$$\forall z$. The mechanical hierarchy that corresponds to the non-linear creep function in Equation (25) is obtained assuming that the applied force $F_0$ is constant, as in standard creep test. In the following, we assume that the axial stiffness and axial damping are expressed by the relation:

$$K_a\left(z,F_0\right)=K_{{\alpha }_c}\frac{F_0{z}^{{}^{-\delta }}}{\mathsf{\Gamma }\left(1-{\alpha }_c\right)}=K_0\frac{F_0^{\left(1-{\alpha }_c\right)}{z}^{-\delta }}{\mathsf{\Gamma }\left(1-\delta \right)}$$

(27)

with $\left[K_0\right]=L^{\delta }{F}^{(2+{\alpha }_c)}$ and $\mathsf{\Gamma }(·)$ is the Euler-Gamma function and $\delta \in [-1,1]$. The non-linear damping coefficient reads:

$$C\left(z,F_0\right)=C_0\frac{F_0^{\left(1-{\alpha }_c\right)}{z}^{{}^{-\delta }}}{\mathsf{\Gamma }\left(1-\delta \right)}$$

(28)

with $\left[C_0\right]=F^{(2+{\alpha }_c)}T{L}^{\delta }$. Substituting Equations (27) and (28) into Equation (25) yields:

$$\frac{\partial }{\partial z}\left[K_0{z}^{-\delta }\frac{\partial u}{\partial z}\right]=C_0{z}^{-\delta }\frac{\partial u}{\partial t}$$

(29)

that has the same mathematical structures as the governing equation of the linear hereditary hierarchy described in [9,10]. Laplace transformation $\mathcal{L}\left[u\right]=\widehat{u}\left(s\right)$ of Equation (29) yields, after some straightforward manipulations:

$$\frac{d}{dz}\left[z^{-\delta }\frac{d\widehat{u}}{dz}\right]=s{\tau }_0{z}^{-\delta }\widehat{u}$$

(30)

with relaxation time ${\tau }_0={\displaystyle \frac{C_0}{K_0}}$. The solution of Equation (30) may be obtained in terms of the first and second modified Bessel functions, $Y_{{\beta }_c}\left(z\sqrt{{\tau }_{0}s}\right)$ and $K_{{\beta }_c}\left(z\sqrt{{\tau }_{0}s}\right)$, with ${\beta }_c(\delta +1)/2$ as:

$$\widehat{u}\left(z,s\right)=z^{{\beta }_c}\left[B_1{Y}_{{\beta }_c}\left(z\sqrt{{\tau }_{0}s}\right)+B_2{K}_{{\beta }_c}\left(z\sqrt{{\tau }_{0}s}\right)\right]$$

(31)

yielding, after substitution into Equation (26), the condition $B_1=0$. The integration constant $B_2$ is obtained by the boundary condition in Equation (26), yielding:

$$B_2\underset{z\to 0}{lim}\left[z^{\left(2{\beta }_c-1\right)}\frac{F_0^{\left(1-{\alpha }_c\right)}{K}_0}{\mathsf{\Gamma }\left(1-\delta \right)}\frac{d}{dz}{z}_c^{\beta }{K}_{{\beta }_c}\left(z\sqrt{{\tau }_{0}s}\right)\right]=\frac{F_0}{s}$$

(32)

The constant $B_2$ is derived as:

$$B_2=\frac{F_0^{{\alpha }_c}{\left(s{\tau }_0\right)}^{-{\beta }_c/\phantom{{\beta }_{c}2}2}}{\mathsf{\Gamma }\left(2{\beta }_c\right)K_0{2}^{2{\beta }_c-1}}$$

(33)

The the axial displacement of the hierarchy at $z\to 0$, namely $u_0\left(t\right)$, corresponds to the material creep functions as:

$$\begin{array}{c}{u}_0\left(t\right)={\mathcal{L}}^{-1}\left[{\displaystyle \frac{F_0^{{\alpha }_c}{\left(s{\tau }_0\right)}^{-{\beta }_c/\phantom{{\beta }_{c}2}2}}{\mathsf{\Gamma }\left(2{\beta }_c\right)K_0{2}^{2{\beta }_c-{\alpha }_c}}}\underset{z\to 0}{lim}{z}^{{\beta }_c}{K}_{{\beta }_c}\left(z\sqrt{{\tau }_{0}s}\right)\right]=\hfill \\ ={\displaystyle \frac{F_0^{{\alpha }_c}{\tau }_0^{-{\beta }_c}}{K_0\mathsf{\Gamma }\left(2{\beta }_c\right)\mathsf{\Gamma }\left({\beta }_c\right)K_0{2}^{2{\beta }_c}}}{t}^{{\beta }_c}={\displaystyle \frac{F_0^{{\alpha }_c}}{K_a\mathsf{\Gamma }{\beta }_c}}{\left({\displaystyle \frac{t}{{\overline{\tau }}_c}}\right)}^{{\beta }_c}\hfill \end{array}$$

(34)

where we denote:

$${\tau }_c={\tau }_0^{{\beta }_c}\mathsf{\Gamma }\left(2{\beta }_c\right)\mathsf{\Gamma }\left(\beta \right)K_a{2}^{{\beta }_c}$$

(35)

that is completely equivalent to the non-linear creep function introduced to capture the non-linear behavior of the tendons observed in the mechanical tests. In passing, we observe that the expression of the creep function in Equation (35) does not allow for the use of time-superposition principle so that no convolution integrals may be defined for non-linear hereditariness provided by Equations (16a) and (16b) unless the non-linear transform $F_e\left(t\right)\to F{\left(t\right)}^{{\alpha }_c}$ is introduced. Under these circumstances, the effect superposition may be used and the fractional-order operators may be readily defined.

In the next section, some numerical applications are provided to show the reliability of the proposed formulation.

4.2. Numerical Analysis

In this section, the hierarchic mechanical model presented in the previous section is used to reproduce experimental results presented in Section 2. To find the response of the hierarchic model to a constant applied stress, we recast the equilibrium equations of the elements of the discretized model (Equation (24)) in matrix form as follows:

$$F_0^{1-{\alpha }_c}\left(p_V\mathbf{B}\dot{\mathbf{u}}+q_V\mathbf{A}\mathbf{u}\right)=\mathbf{v}{F}_0$$

(36)

where $\mathbf{v}$ is an influence vector and

$$p_V=\frac{C_0\Delta z^{1-\delta }}{\mathsf{\Gamma }(1+\delta )}$$

(37a)

$$q_V=\frac{K_0\Delta z^{-1-\delta }}{\mathsf{\Gamma }(1-\delta )}$$

(37b)

$$\mathbf{A}=\left[\begin{array}{ccccccc}{1}^{-\delta }& -1^{-\delta }& 0& ·& ·& ·& 0\\ -1^{-\delta }& 1^{-\delta }+2^{-\delta }& -2^{-\delta }& ·& ·& ·& 0\\ 0& -2^{-\delta }& 2^{-\delta }+3^{-\delta }& ·& ·& ·& 0\\ ·& ·& ·& ·& & & ·\\ ·& ·& ·& & ·& & ·\\ ·& ·& ·& & & ·& -{(N-1)}^{-\delta }\\ 0& 0& 0& ·& ·& ·& {(N-1)}^{\delta }+N^{-\delta }\end{array}\right]$$

(38a)

$$\mathbf{B}=\left[\begin{array}{ccccccc}{1}^{-\delta }& 0& 0& ·& ·& ·& 0\\ 0& 2^{-\delta }& 0& ·& ·& ·& 0\\ 0& 0& 3^{-\delta }& ·& ·& ·& 0\\ ·& ·& ·& ·& & & ·\\ ·& ·& ·& & ·& & ·\\ ·& ·& ·& & & ·& ·\\ 0& 0& 0& ·& ·& ·& N^{-\delta }\end{array}\right]$$

(38b)

and N is the number of laminae in which the hierarchical model has been discretized. To evaluate the solution set of ordinary differential equations in Equation (36), we premultiply the system by ${\mathbf{B}}^{-1/2}$ and introduce the coordinates transformation as $\mathbf{x}={\mathbf{B}}^{1/2}\mathbf{u}$; as a consequence, we can write Equation (36) as

$$p_V\dot{\mathbf{x}}+q_V\mathbf{D}\mathbf{x}=\tilde{\mathbf{v}}{F}_0^{{\alpha }_c}$$

(39)

being $\mathbf{D}={\mathbf{B}}^{-1/2}\mathbf{A}{\mathbf{B}}^{-1/2}$ and $\tilde{\mathbf{v}}={\mathbf{B}}^{-1/2}\mathbf{v}$. Moreover, in Equation (39), the load terms are simplified and the load appears only on the right-hand side of the system; as a consequence, for constant applied load history, the system is linear. For this reason, the matrix $\mathbf{D}$ can be easily diagonalized by the standard method of modal analysis [10]. To this purpose, let $\mathsf{\Phi }$ be the modal matrix whose columns are the eigenvectors of the matrix D. Then,

$${\mathsf{\Phi }}^T\mathbf{D}\mathsf{\Phi }=\mathsf{\Lambda }{\mathsf{\Phi }}^T\mathsf{\Phi }=\mathbf{I}$$

(40)

being $\mathbf{I}$ the identity matrix and $\mathsf{\Lambda }$ a diagonal matrix containing the eigenvectors of $\mathbf{D}$. By introducing the modal coordinates in Equation (39) as

$$\mathbf{x}\left(t\right)=\mathsf{\Phi }\mathbf{y}\left(t\right)\mathbf{y}\left(t\right)={\mathsf{\Phi }}^T\mathbf{x}\left(t\right)$$

(41)

the following system of uncoupled linear differential equations is obtained

$$p_V\dot{\mathbf{y}}+q_V\mathsf{\Lambda }\mathbf{x}=\overline{\mathbf{v}}{F}_0^{{\alpha }_c}$$

(42)

The generic jth equation reads

$${\delta }_j{\dot{y}}_j+y_j=\frac{{\varphi }_{1,j}}{q_V{\lambda }_j}{F}_0^{{\alpha }_c}$$

(43)

where ${\delta }_j=p_V/q_V{\lambda }_j>0$ and ${\varphi }_{1,j}$ is the first component of the jth eigenvector. Equation (43) represents the governing equation of the well-known viscoelastic Kelvin–Voigt model with unitary stiffness coefficient and damping coefficient equal to ${\delta }_j$. Since the applied force on the right-hand side of Equation (43) is constant, the solution is the creep function of the Kelvin–Voigt model amplified by the value of the force in the y domain, that is

$$y_j\left(t\right)=\frac{{\varphi }_{1,j}}{q_V{\lambda }_j}{F}_0^{{\alpha }_c}\left[1-exp(-t/{\delta }_j)\right]$$

(44)

The solution in the original domain related to the upper lamina, the one we are interested in, can easily be obtained by the following

$$u_1\left(t\right)={\mathbf{v}}^T{\mathbf{B}}^{-1/2}\mathsf{\Phi }\mathbf{y}\left(t\right)$$

(45)

Equation (45) can reproduce the nonlinear creep behavior described by the Nutting law. To this purpose, Figure 4 shows the creep function obtained with the aid of Equation (45) for different values of the applied force $F_0$. It is evident that, as the applied force increases, the creep curves do not scale linearly.

5. Conclusions

In this paper, the authors discuss in detail the non-linear hereditariness of materials, described by time non-local models, with a single-integral approach. As long as time-varying power-law is considered as integral kernel, a fractional-order single integral non-linear model is obtained with a three parameters generalized springpot: (i) non-linear exponent ${\alpha }_c$; (ii) decaying exponent ${\beta }_c$; and (iii) the characteristic time ${\tau }_c$. The model is completely equivalent to the Nutting relations observed for rubbers and polymers. The single integral model allows for the use of effect superposition of auxiliary state variables, in terms of stress and strain, that account for the position of some equivalence relations among creep and relaxations parameters: (i) non-linear exponent ${\alpha }_r$; (ii) decaying exponent ${\beta }_r$; and (iii) characteristic time ${\tau }_r$. They were observed and verified with a large experimental campaign on several materials. The generalized springpot possesses an equivalent mechanical hierarchy, totally analogous to the well-known mechanical hierarchy already established for the linear springpot, with non-linear internal springs and dashpots. Some numerical examples reporting the behavior of the non-linear hierarchy are reported in the paper.