*Article* **A General Return-Mapping Framework for Fractional Visco-Elasto-Plasticity**

**Jorge L. Suzuki 1, Maryam Naghibolhosseini <sup>2</sup> and Mohsen Zayernouri 3,\***


**Abstract:** We develop a fractional return-mapping framework for power-law visco-elasto-plasticity. In our approach, the fractional viscoelasticity is accounted for through canonical combinations of Scott-Blair elements to construct a series of well-known fractional linear viscoelastic models, such as Kelvin–Voigt, Maxwell, Kelvin–Zener, and Poynting–Thomson. We also consider a fractional quasilinear version of Fung's model to account for stress/strain nonlinearity. The fractional viscoelastic models are combined with a fractional visco-plastic device, coupled with fractional viscoelastic models involving serial combinations of Scott-Blair elements. We then develop a general returnmapping procedure, which is fully implicit for linear viscoelastic models, and semi-implicit for the quasi-linear case. We find that, in the correction phase, the discrete stress projection and plastic slip have the same form for all the considered models, although with different property and timestep-dependent projection terms. A series of numerical experiments is carried out with analytical and reference solutions to demonstrate the convergence and computational cost of the proposed framework, which is shown to be at least first-order accurate for general loading conditions. Our numerical results demonstrate that the developed framework is more flexible and preserves the numerical accuracy of existing approaches while being more computationally tractable in the viscoplastic range due to a reduction of 50% in CPU time. Our formulation is especially suited for emerging applications of fractional calculus in bio-tissues that present the hallmark of multiple viscoelastic power-laws coupled with visco-plasticity.

**Keywords:** power-law visco-elasto-plasticity; time-fractional integration; fractional quasi-linear viscoelasticity

### **1. Introduction**

Power-law behavior has been observed in living cells [1,2] and bio-tissues [3–5]. This stems from the ubiquitous self-similar and scale-free nature of the tissue/cell microstructure, which can be physically and mathematically scaled up to continuum level, manifesting in the power-law behavior in the lumped sense. Such power-law relationships have been seen in auditory hair cells, positioned in the sensory organ of hearing, cochlea [6], vocal fold tissues [7], and bladder tissues [5]. Experimental evidence suggests that complex material behavior may possess more than a single power-law scaling in the viscoelastic regime, particularly in multi-fractal structures, such as in cells [8] and biological tissues [9], due to their complex, hierarchical, and heterogeneous microstructures. For such cases, a single fractional rheological element is not sufficient to capture the observed behavior even if the data suggest a linear viscoelastic behavior. Stamenovi´c et al. [8] measured the complex shear modulus of cultured human airway smooth muscle and observed two distinct power-law regimes, separated by an intermediate *plateau*. Kapnistos et al. [10]

**Citation:** Suzuki, J.L.;

Naghibolhosseini, M.; Zayernouri, M. A General Return-Mapping Framework for Fractional Visco-Elasto-Plasticity. *Fractal Fract.* **2022**, *6*, 715. https://doi.org/ 10.3390/fractalfract6120715

Academic Editors: Angelo B. Mingarelli, Leila Gholizadeh Zivlaei and Mohammad Dehghan

Received: 20 September 2022 Accepted: 8 November 2022 Published: 1 December 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

found an unexpected tempered power-law relaxation response of entangled polystyrene ring polymers compared with the usual relaxation plateau of linear chain polymers. Such behavior was interpreted through self-similar conformations of double-folded loops in the ring polymer, instead of the repetition observed in linear chains.

In addition to multiple viscoelastic power-law behaviors, there also exists evidence of bio-plasticity in soft media [11,12]. The creep behavior of human embryonic stem cells under differentiation was studied by Pajerowski et al. [11] through micro-aspiration experiments at different pressures. The cell nucleus demonstrated distinguished viscoelasto-plastic power-law scalings, with *α* = 0.2 for the plastic regime, independent of the applied pressure. It is discussed that such low power-law exponent arises due to the fractal arrangement of chromatin inside the cell nucleus. Studies on force-induced mechanical plasticity of mouse embryonic fibroblasts were performed by Bonadkar et al. [12]. They found that the viscoelastic relaxation and the permanent deformations followed a stochastic, normally distributed, power-law scaling *β*(*ω*) with values ranging from *β* ≈ 0 to *β* ≈ 0.6. The microstructural mechanism of plastic deformation in the cytoskeleton is due to the combination of permanent stretching and buckling of actin fibers.

Regarding existing modeling approaches of anomalous plasticity, several works employed fractional calculus to account for the visco-plastic regimes of different classes of materials [13]. Three of the main approaches include: time-fractional, space-fractional, and stress-fractional modeling. The time-fractional approaches focus on introducing memory effects into non-equilibrium viscous variables [14,15] and consequently modeling power-laws in both viscoelastic and visco-plastic regimes, which is applied for polymers, cells, and tissues. Suzuki et al. [14] developed a fractional visco-elasto-plastic model that provides a constitutive interpolation between rate-independent plasticity and Perzyna's visco-plasticity by introducing a Scott-Blair (SB) model acting the plastic regime. This model utilizes a rate-dependent-yield function, which was later proved to be thermodynamically consistent in a further extension of the model to account for continuum damage mechanics [16]. A three-dimensional space-fractional approach to elasto-plasticity was also developed by Sumelka [17] in order to consider the spatial nonlocalities. The model is based on rate-independent elasto-plasticity, and nonlocal effects are modeled using a fractional continuum mechanics approach, where the strains are defined through a space-fractional Riesz–Caputo derivative of the displacements. Finally, the stress-fractional models for plasticity have been found to be applicable for modeling the soil mechanics and geomaterials that follow a non-associated plastic flow [18,19], i.e.,in which the yield surface expansion in the stress space does not follow the usual normality rule and may be nonconvex. Sumelka [18] proposed a three-dimensional fractional visco-plastic model, where a fractional flow-rule with the order 0 < *α* < 1 in the stress domain naturally modeled the non-associative plasticity. This model recovers the classical Perzyna visco-plasticity as *α* → 1, and the effect of the fractional flow rule can be a compact descriptor of microstructure anisotropy. Later on, Sun and Sumelka [19] developed a similar stress-fractional model, which was successfully applied for soils under compression. We refer the reader to the Sun et al. review work on fractional calculus applications in plasticity [20].

In this work, we develop a generalized fractional visco-elasto-plastic model, where the visco-plastic device can be coupled with several existing fractional linear/nonlinear viscoelastic representations. More specifically, we utilize a fractional visco-plastic device developed in [14,16], which is then coupled with a series of linear fractional models, such as Scott-Blair (SB), Kelvin–Voigt (FKV), Maxwell (FM), Kelvin–Zener (FKZ), Poynting– Thomson (FPT); also a fractional quasi-linear viscoelastic (FQLV) model for large strains Figure 1. Then, a generalized fractional return-mapping algorithm is proposed, which overcomes existing difficulties in previous developments by first fully discretizing all fractional operators and then performing the predictor–corrector procedure. More specifically, the existing approaches are built on the notion of employing the predictor–corrector approach before the full discretization of fractional operators while treating trial states for stress and internal variables as continuous functions of time. This prevents models with serial combinations of SB elements to be incorporated in associated yield functions in a straightforward fashion. The main features of the proposed framework are:


We carry out a number of numerical experiments involving fabricated and reference solutions under monotone and general loading conditions and observe a global accuracy ranging from O(Δ*t*) to O(Δ*t* <sup>2</sup>−*β*), according to the regularity induced by the associated fractional differential equations (FDEs) and loading conditions.

This work is organized as follows. In Section 2, we present the mathematical definitions employed in this work. In Section 3, we describe the considered linear/quasi-linear fractional viscoelastic models, coupled with fractional visco-elasto-plasticity as explained in Section 4. All corresponding models are discretized and posed in a unified fractional return-mapping form in Section 5. Convergence analyses and computational performance evaluation of presented models and return-mapping algorithms are performed in Section 6, followed by the concluding remarks in Section 7.

**Figure 1.** Fractional linear viscoelastic models employed in this work, constructed from serial/parallel combinations of "building block" SB elements. The SB building blocks naturally account for an infinite fractal arrangement of Hookean/Newtonian elements. The employed fractional quasilinear model is not represented by a mechanical analogue although the time-dependent component of the relaxation function has an SB-like representation.

#### **2. Definitions of Fractional Calculus**

We start with some preliminary definitions of fractional calculus [21]. The left-sided Riemann–Liouville integral of order *β* ∈ (0, 1) is defined as

$$(\prescript{RL}{t\_{\mathbb{L}}} \mathcal{Z}\_t^{\mathbb{\beta}} f)(t) = \frac{1}{\Gamma(\beta)} \int\_{t\_{\mathbb{L}}}^t \frac{f(s)}{(t - s)^{1 - \beta}} \, ds, \quad t > t\_{L'} \tag{1}$$

where Γ represents the Euler gamma function, and *tL* denotes the lower integration limit. The corresponding inverse operator, the left-sided fractional derivative of order *β*, is then defined based on (1) as

$$(\prescript{RL}{t\_L}{\mathcal{D}}\_t^{\mathcal{S}} f)(t) = \frac{d}{dt}(\prescript{RL}{t\_L}{\mathcal{T}}\_t^{1-\mathcal{S}} f)(t) = \frac{1}{\Gamma(1-\beta)} \frac{d}{dt} \int\_{t\_L}^t \frac{f(s)}{(t-s)^{\beta}} \, ds, \quad t > t\_L. \tag{2}$$

The left-sided Caputo derivative of order *β* ∈ (0, 1) is obtained as

$$(\,^C\_{t\!\!L} \mathcal{D}\_t^{\beta} f)(t) = (\,^{RL}\_{t\!\!L} \mathcal{D}\_t^{1-\beta} \frac{df}{dt})(t) = \frac{1}{\Gamma(1-\beta)} \int\_{t\!\!L}^t \frac{f'(s)}{(t-s)^{\beta}} \, ds, \quad t > t\_L. \tag{3}$$

The definitions of Riemann–Liouville and Caputo derivatives are linked by the following relationship:

$$(\prescript{RL}{t\_{\rm L}} \mathcal{D}\_t^{\mathcal{S}} f)(t) = \frac{f(t\_{\rm L})}{\Gamma(1-\beta)(t+t\_{\rm L})^{\mathcal{S}}} + (\prescript{C}{t\_{\rm L}} \mathcal{D}\_t^{\mathcal{S}} f)(t),\tag{4}$$

which can be obtained through integration by parts followed by the application of the Leibniz rule on (2). We should note that the aforementioned derivatives coincide when dealing with homogeneous Dirichlet initial/boundary conditions. Finally, we define the two-parameter Mittag–Leffler function *Ea*,*b*(*z*) [22] as:

$$E\_{a,b}(z) = \sum\_{k=0}^{\infty} \frac{z^k}{\Gamma(ak+b)}, \quad \text{Re}(a) > 0, \quad b \in \mathbb{C}, \quad z \in \mathbb{C}.\tag{5}$$

#### **3. Fractional Viscoelasticity**

We present the linear and quasi-linear fractional viscoelastic models that we couple with the visco-plastic return-mapping procedure.

#### *3.1. Linear Viscoelasticity*

**Scott-Blair (SB) Model:** The rheological *building block* for our framework is the fractional SB viscoelastic element, which compactly represents an anomalous viscoelastic constitutive law relating the stresses and strains:

$$
\sigma(t) = \mathbb{E}^{\mathbb{C}}\_0 \mathcal{D}\_t^{\mathbb{B}} \varepsilon(t), \quad t > 0, \quad \varepsilon(0) = 0,\tag{6}
$$

with pseudo-constant <sup>E</sup><sup>1</sup> [*Pa*.*sβ*] <sup>≥</sup> 0, and constant fractional order 0 <sup>&</sup>lt; *<sup>β</sup>* <sup>&</sup>lt; 1, which provides a material interpolation between the Hookean (*β* → 0) and Newtonian (*β* → 1) elements. The pair (*β*,E) uniquely represents the SB constants, where the *pseudo-constant* E [*Pa*.*sβ*] compactly describes textural properties, such as the firmness of the material [23,24]. In this sense, E is interpreted as a snapshot of a non-equilibrium dynamic process instead of an equilibrium state. The corresponding rheological symbol for the SB model represents a fractal-like arrangement of springs and dashpots [25,26], which we interpret as a compact, upscaled representation of a fractal-like microstructure. Regarding the thermodynamic admissibility, we refer the reader to Lion [27] for the SB model and Suzuki et al. [16] for the combination of the SB element with more complex mechanisms of visco-plasticity and damage. The relaxation function *G*(*t*) [*Pa*] for the SB model is given by the following inverse power-law form:

$$G^{SB}(t) := \frac{\mathbb{E}}{\Gamma(1-\alpha)} t^{-\beta} \, \_{\prime} \tag{7}$$

which is the convolution kernel of the differ-integral form in (6).

**Fractional Kelvin–Voigt (FKV) Model:** Through a parallel combination of SB elements, we obtain the following stressd–strain relationship [25]:

$$
\sigma(t) = \mathbb{E}\_1 \, \prescript{\mathbb{C}}{0} \mathcal{D}\_t^{\beta\_1} \varepsilon(t) + \mathbb{E}\_2 \, \prescript{\mathbb{C}}{0} \mathcal{D}\_t^{\beta\_2} \varepsilon(t), \quad t > 0, \quad \varepsilon(0) = 0,\tag{8}
$$

with fractional orders 0 <sup>&</sup>lt; *<sup>β</sup>*1, *<sup>β</sup>*<sup>2</sup> <sup>&</sup>lt; 1, and associated pseudo-constants <sup>E</sup><sup>1</sup> [*Pa*.*sβ*<sup>1</sup> ] <sup>≥</sup> 0, and <sup>E</sup><sup>2</sup> [*Pa*.*sβ*<sup>2</sup> ] <sup>≥</sup> 0. The corresponding relaxation modulus *<sup>G</sup>*(*t*) [*Pa*] is also an additive form of two SB elements:

$$G^{FkV}(t) := \frac{\mathbb{E}\_1}{\Gamma(1-\beta\_1)} t^{-\beta\_1} + \frac{\mathbb{E}\_2}{\Gamma(1-\beta\_2)} t^{-\beta\_2},\tag{9}$$

which has a response characterized by two power-law regimes with a transition from faster to slower relaxation. Assuming *β*<sup>2</sup> > *β*1, the asymptotic responses for small and large time-scales are given by *<sup>G</sup>FKV* <sup>∼</sup> *<sup>t</sup>* <sup>−</sup>*β*<sup>2</sup> as *<sup>t</sup>* <sup>→</sup> 0 and *<sup>G</sup>FKV* <sup>∼</sup> *<sup>t</sup>* <sup>−</sup>*β*<sup>1</sup> as *<sup>t</sup>* <sup>→</sup> <sup>∞</sup>.

**Fractional Maxwell (FM) Model:** Through a serial combination of SB elements, we obtain the fractional Maxwell (FM) model [24], given by:

$$
\sigma(t) + \frac{\mathbb{E}\_2}{\mathbb{E}\_1} \, ^\mathbb{C}\_0 \mathcal{D}\_t^{\mathbb{A}\_2 - \beta\_1} \sigma(t) = \mathbb{E}\_2 \, ^\mathbb{C}\_0 \mathcal{D}\_t^{\mathbb{A}\_2} \varepsilon(t), \quad t > 0,\tag{10}
$$

with pseudo-constants <sup>E</sup><sup>1</sup> [*Pa*.*sβ*<sup>1</sup> ] <sup>&</sup>gt; 0 and <sup>E</sup>2[*Pa*.*sβ*<sup>1</sup> ] <sup>≥</sup> 0, fractional orders 0 <sup>&</sup>lt; *<sup>β</sup>*<sup>1</sup> <sup>&</sup>lt; *β*<sup>2</sup> < 1 with ,0 < *β*<sup>2</sup> − *β*<sup>1</sup> < 1 and two sets of initial conditions for strains *ε*(0) = 0 and stresses *σ*(0) = 0. We should note that in the case of non-homogeneous initial conditions, there needs to be a compatibility condition [22] between stresses and strains at *t* = 0. The corresponding relaxation function for this building block model assumes a more complex Miller–Ross form [24]:

$$G^{FM}(t) := \mathbb{E}\_1 t^{-\beta\_1} E\_{\beta\_2 - \beta\_1, 1 - \beta\_1} \left( -\frac{\mathbb{E}\_1}{\mathbb{E}\_2} t^{\beta\_2 - \beta\_1} \right). \tag{11}$$

The presence of a Mittag–Leffler function in (11) leads to a stretched exponential relaxation for smaller times and a power-law behavior for longer times. We also observe that the limit cases are given by *<sup>G</sup>FM* <sup>∼</sup> *<sup>t</sup>* <sup>−</sup>*β*<sup>1</sup> as *<sup>t</sup>* <sup>→</sup> 0 and *<sup>G</sup>FM* <sup>∼</sup> *<sup>t</sup>* <sup>−</sup>*β*<sup>2</sup> as *<sup>t</sup>* <sup>→</sup> <sup>∞</sup>, indicating that the FM model provides a behavior transitioning from slower-to-faster relaxation. We refer the reader to [5,24,28] for a number of applications of the aforementioned models. We should notice that both FKV and FM models are able to recover the SB element with a convenient set of pseudo-constants and *β*<sup>1</sup> = *β*2.

**Fractional Kelvin–Zener (FKZ) model:** The fractional generalization of the standard linear solid (SLS) model is given by an FM branch in parallel with a third SB element, given by the following FDE:

$$\left[1+\frac{\mathbb{E}\_2}{\mathbb{E}\_1}\,^{\mathbb{C}}\_0\mathcal{D}\_t^{\mathcal{S}\_2-\mathcal{S}\_1}\right]\sigma(t) = \left[\mathbb{E}\_2\,^{\mathbb{C}}\_0\mathcal{D}\_t^{\mathcal{S}\_2}+\mathbb{E}\_3\,^{\mathbb{C}}\_0\mathcal{D}\_t^{\mathcal{S}\_3}+\frac{\mathbb{E}\_2\mathbb{E}\_3}{\mathbb{E}\_1}\,^{\mathbb{C}}\_0\mathcal{D}\_t^{\mathcal{S}\_2+\mathcal{S}\_3-\mathcal{S}\_1}\right]\varepsilon(t),\tag{12}$$

with fractional orders 0 < *β*<sup>1</sup> < *β*<sup>2</sup> < 1 and conditions 0 < *β*<sup>2</sup> − *β*<sup>1</sup> < 1 and 0 < *<sup>β</sup>*<sup>2</sup> <sup>+</sup> *<sup>β</sup>*<sup>3</sup> <sup>−</sup> *<sup>β</sup>*<sup>1</sup> <sup>&</sup>lt; 1, pseudo-constants <sup>E</sup><sup>1</sup> [*Pa*.*sβ*<sup>1</sup> ] <sup>&</sup>gt; 0, <sup>E</sup><sup>2</sup> [*Pa*.*sβ*<sup>2</sup> ] <sup>≥</sup> 0 and <sup>E</sup><sup>3</sup> [*Pa*.*sβ*<sup>3</sup> ] <sup>≥</sup> 0, and the same initial conditions as in the FM model. We should note that the FM model is recovered when E<sup>3</sup> = 0, and the FKV model is recovered when setting E<sup>1</sup> = 0. The relaxation function is obtained in a straightforward fashion as the summation of the relaxation functions from the SB and FM models:

$$G^{FKZ}(t) := \mathbb{E}\_1 t^{-\beta\_1} E\_{\beta\_2 - \beta\_1, 1 - \beta\_1} \left( -\frac{\mathbb{E}\_1}{\mathbb{E}\_2} t^{\beta\_2 - \beta\_1} \right) + \frac{\mathbb{E}\_3}{\Gamma(1 - \beta\_3)} t^{-\beta\_3},\tag{13}$$

which leads to three inverse power-law regimes for short, intermediate, and long times with particular relationships between *β*1, *β*2, *β*<sup>3</sup> [29].

**Fractional Poynting–Thomson (FPT) Model:** Finally, we introduce our last fractional linear viscoelastic model, given by the serial combination between an FKV model and an SB element:

$$
\sigma \left[ 1 + \frac{\mathbb{E}\_1}{\mathbb{E}\_3} \mathbb{\hat{G}} \mathcal{D}\_t^{\mathfrak{E}\_1 - \mathfrak{E}\_3} + \frac{\mathbb{E}\_2}{\mathbb{E}\_3} \mathbb{\hat{G}} \mathcal{D}\_t^{\mathfrak{E}\_2 - \mathfrak{E}\_3} \right] \sigma(t) = \left[ \mathbb{E}\_1 \, \prescript{\mathbb{C}}{}\_0 \mathcal{D}\_t^{\mathfrak{E}\_1} + \mathbb{E}\_2 \, \prescript{\mathbb{C}}{}\_0 \mathcal{D}\_t^{\mathfrak{E}\_2} \right] \varepsilon(t), \tag{14}
$$

with 0 < *β*<sup>3</sup> < *β*<sup>1</sup> < 1 and 0 < *β*<sup>3</sup> < *β*<sup>2</sup> < 1, additional conditions 0 < *β*<sup>1</sup> − *β*<sup>3</sup> < 1 and <sup>0</sup> <sup>&</sup>lt; *<sup>β</sup>*<sup>2</sup> <sup>−</sup> *<sup>β</sup>*<sup>3</sup> <sup>&</sup>lt; 1, and pseudo-constants <sup>E</sup><sup>1</sup> [*Pa*.*sβ*<sup>1</sup> ] <sup>≥</sup> 0, <sup>E</sup><sup>2</sup> [*Pa*.*sβ*<sup>2</sup> ] <sup>≥</sup> 0, and <sup>E</sup><sup>3</sup> [*Pa*.*sβ*<sup>3</sup> ] <sup>&</sup>gt; <sup>0</sup> and homogeneous initial conditions *σ*(0) = 0 and *ε*(0) = 0. Similar to the FKZ model, we recover the FM model when setting either E<sup>1</sup> or E<sup>2</sup> to zero; although the FKV model cannot be recovered except for a trivial case when *σ*(*t*) = 0.

#### *3.2. Quasi-Linear Fractional Viscoelasticity*

Although fractional linear viscoelastic models provide suitable relaxation functions that describe the anomalous viscoelastic dynamics of a number of soft materials, at times, complex microstructural deformation mechanisms and large strains induce material nonlinearities; hence, the relaxation function itself depends on the applied strain levels. To incorporate this additional effect, we also consider the following FQLV model [30,31]:

$$
\sigma(\mathbf{t}, \varepsilon) = \int\_0^t \mathbf{G}(\mathbf{t} - \mathbf{s}) \frac{\partial \sigma^\varepsilon(\varepsilon)}{\partial \varepsilon} \mathbf{\dot{z}} \, d\mathbf{s}, \tag{15}
$$

where the convolution kernel is given by a multiplicative decomposition of a reduced relaxation function *G*(*t*) and an instantaneous nonlinear elastic tangent response with stress *σ<sup>e</sup>* . In the work by Craiem et al. [31], the reduced relaxation function has a fractional Kelvin–Voigt-like form with one of the SB replaced with a Hookean element. Here, we assume a simpler rheology and adopt an SB-like reduced relaxation in the form:

$$G(t) = Et^{-\alpha}/\Gamma(1-\alpha)\tag{16}$$

with the pseudo-constant *E* with units [*sα*]. We adopt the same two-parameter exponential nonlinear elastic part as in [31]:

$$
\sigma^{\varepsilon}(\varepsilon) = A \left( \varepsilon^{B\varepsilon} - 1 \right),
\tag{17}
$$

with *A* having units of [*Pa*]. Plugging in (16) and (17) into (15), we obtain:

$$
\sigma(t,\varepsilon) = \frac{EAB}{\Gamma(1-a)} \int\_0^t \frac{e^{R\varepsilon(s)}\dot{\varepsilon}(s)}{(t-s)^a} \, ds,\tag{18}
$$

which differs slightly from the linear SB model (6) in the sense that an additional exponential factor multiplies the function being convoluted.

#### **4. Fractional Visco-Elasto-Plasticity**

With all fractional viscoelastic models defined in Section 3, we can couple any of them, subject to a viscoelastic strain *εve*(*t*), to the fractional visco-plastic device, illustrated in Figure 2. The visco-plastic device is composed of a parallel combination of a Coulomb element with initial yield stress *σ<sup>Y</sup>* [*Pa*], an SB element with pseudo-constant K [*Pa*.*sβ<sup>K</sup>* ], and fractional order *βK*, and a Hookean spring with constant *H* [*Pa*]. The entire visco-plastic part is subject to a visco-plastic strain *<sup>ε</sup>vp*(*t*) : <sup>R</sup><sup>+</sup> <sup>→</sup> <sup>R</sup>. In order to obtain the kinematic equations for the internal variables, we start with an additive decomposition of the total logarithmic strain *<sup>ε</sup>*(*t*) : <sup>R</sup><sup>+</sup> <sup>→</sup> <sup>R</sup> acting on the visco-elasto-plastic device:

$$
\varepsilon(t) = \varepsilon^{\upsilon \varepsilon}(t) + \varepsilon^{\upsilon p}(t) \tag{19}
$$

The visco-plastic effects are accounted for through the definition of a memory- and rate-dependent-yield function *<sup>f</sup>*(*σ*, *<sup>α</sup>*) : <sup>R</sup> <sup>×</sup> <sup>R</sup><sup>+</sup> <sup>→</sup> <sup>R</sup><sup>−</sup> ∪ {0} in the following form [14]:

$$f(\sigma, \mathfrak{a}) := |\sigma| - \left[\sigma^Y + \mathbb{K}\_0^C \mathcal{D}\_t^{\mathcal{S}\_K}(\mathfrak{a}) + H\mathfrak{a}\right]. \tag{20}$$

where *<sup>α</sup>* <sup>∈</sup> <sup>R</sup><sup>+</sup> represents the internal hardening variable, and the above form accounts for the isotropic hardening. The set of admissible stresses lie in a closed convex space, where the associated boundary respects the yield condition of classical plasticity (see [16], Lemma 4.1, setting the damage as *D* = 0). From the defined yield function (20) and the principle of maximum plastic dissipation [32], the following properties hold: (i) associativity of the flow

rule, (ii) associativity in the hardening law, (iii) Kuhn–Tucker complimentary conditions, and (iv) convexity. The set of evolution equations for the internal variables *εvp* and *α* is obtained by:

$$\varepsilon^{\upsilon\upsilon} = \frac{\partial f}{\partial \sigma} \gamma\_{\prime} \quad \text{if} \quad = -\frac{\partial f}{\partial R} \dot{\gamma}\_{\prime} \tag{21}$$

where *<sup>γ</sup>*˙(*t*) : <sup>R</sup><sup>+</sup> <sup>→</sup> <sup>R</sup><sup>+</sup> denotes the plastic slip rate. Evaluating the above equations using (20), we obtain the evolution for visco-plastic strains and hardening [14]:

$$
\mathfrak{s}^{vp} = \text{sign}(\mathfrak{r})\gamma. \tag{22}
$$

$$
\dot{\mathfrak{a}} = \dot{\gamma}.\tag{23}
$$

**Proposition 1.** *The closure for the plastic slip rate <sup>γ</sup>*˙(*t*) <sup>∈</sup> <sup>R</sup><sup>+</sup> *with an SB viscoelastic part of constants* (E, *βE*)*,* (K, *βK*)*, and H (model M1 [14]) with homogeneous initial conditions for the internal variables and their respective rates, i.e., εvp*(0) = *α*(0) = *γ*(0) = 0*, γ*˙(0) = 0*, and ε*˙ *vp*(0) = *α*˙(0) = *γ*˙(0) = 0*, is given by the following fractional Cauchy problem:*

$$\mathbb{E}\,^{\mathbb{C}}\_{0}\mathcal{D}^{\mathbb{\beta}\_{E}}\_{t}\dot{\gamma}(t) + \mathbb{K}\,^{\mathbb{C}}\_{0}\mathcal{D}^{\mathbb{\beta}\_{K}}\_{t}\dot{\gamma}(t) + H\dot{\gamma}(t) = \text{sign}(\sigma)\mathbb{E}\left[\frac{\dot{\varepsilon}(0)t^{-\beta\_{E}}}{\Gamma(1-\beta\_{E})} + \,^{\mathbb{C}}\_{0}\mathcal{D}^{\mathbb{\beta}\_{E}}\_{t}\dot{\varepsilon}(t)\right] \tag{24}$$

**Proof.** See Appendix A.

**Figure 2.** Fractional visco-elasto-plastic model subject to a uniaxial stress *σ*. Here, any of the linear and quasi-linear (not illustrated) fractional viscoelastic models under strains *εve* can be chosen and then coupled with a fractional viscoplastic rheological device under strains *εvp*.

#### **5. A Class of Return-Mapping Algorithms for Fractional Visco-Elasto-Plasticity**

Given the presented viscoelastic and visco-plastic models, respectively, in Sections 3 and 4, we now demonstrate how to solve each resulting system of nonlinear equations according to the choice of viscoelastic models. The considered fractional return-mapping approach in this work is fully discrete, i.e., we first discretize all fractional derivatives using a finitedifference approach and then employ trial states for the internal variables in a predictor– corrector scheme.

We discretize the fractional Caputo derivatives in (6)–(10) through an implicit L1 finite-difference scheme [33]. Extensions to account for fast time-stepping approaches [34] are straightforward, since they mostly affect the history terms computation. Let Ω = (0, *T*] be decomposed into a uniform time grid with *N* time steps of size Δ*t*, such that *tn* = *n*Δ*t*, with *n* = 0, 1, ... , *N*. The time-fractional Caputo derivative of a real-valued function *<sup>u</sup>*(*t*) <sup>∈</sup> *<sup>C</sup>*2(Ω) at time *<sup>t</sup>* <sup>=</sup> *tn*+<sup>1</sup> is therefore discretized as [33]:

$$\left. \begin{matrix} {\_0^C} \mathcal{D}\_t^\beta u(t) \end{matrix} \right|\_{t=t\_{n+1}} = \frac{1}{\Delta t^\beta \Gamma(2-\beta)} [\boldsymbol{\mu}\_{n+1} - \boldsymbol{\mu}\_n + \mathcal{H}^\boldsymbol{u} \boldsymbol{u}] + \mathcal{O}(\Delta t^{2-\beta}),\tag{25}$$

with the history term <sup>H</sup>*ν<sup>u</sup>* given by the following form:

$$\mathcal{H}^{\mathcal{G}}u = \sum\_{j=1}^{n} b\_j^{\mathcal{G}} \left[ u\_{n+1-j} - u\_{n-j} \right] \tag{26}$$

with weights *b β <sup>j</sup>* := (*<sup>j</sup>* <sup>+</sup> <sup>1</sup>)1−*<sup>β</sup>* <sup>−</sup> *<sup>j</sup>* <sup>1</sup>−*β*.

#### *5.1. Time-Fractional Integration of Viscoelastic Models*

In the following, we present the discretized forms for each considered fractional viscoelastic model from Section 3, which are represented in a fully implicit fashion.

**Scott-Blair Model**: Evaluating both sides of (6) at *t* = *tn*+1, we obtain:

$$
\sigma\_{n+1} = \mathbb{E}\_1 \, \_0^C \mathcal{D}\_t^{\mathcal{B}\_1} \varepsilon(t) \big|\_{t=t\_{n+1}} \tag{27}
$$

in which by applying (25), we directly obtain:

$$\sigma\_{n+1} = \mathbb{C}\_1^{SB} \left[ \varepsilon\_{n+1} - \varepsilon\_n + \mathcal{H}^{\beta\_1} \varepsilon \right] \tag{28}$$

with the strain history <sup>H</sup>*β*<sup>1</sup> *<sup>ε</sup>* and constant *<sup>C</sup>SB* <sup>1</sup> , shown in Appendix B for the SB and the following model discretizations.

**Fractional Kelvin–Voigt Model**: Evaluating both sides of (8) at *t* = *tn*+1, we obtain:

$$
\sigma\_{n+1} = \mathbb{E}\_0^C \mathcal{D}\_t^{\beta\_1} \varepsilon(t) \big|\_{t=t\_{n+1}} + \mathbb{E}\_0^C \mathcal{D}\_t^{\beta\_2} \varepsilon(t) \big|\_{t=t\_{n+1}} \tag{29}
$$

which, applying (25) for the fractional derivatives of order *β*<sup>1</sup> and *β*2, leads to:

$$
\sigma\_{n+1} = \mathbb{C}\_1^{KV} \left[ \varepsilon\_{n+1} - \varepsilon\_n + \mathcal{H}^{\mathcal{S}\_1} \varepsilon \right] + \mathbb{C}\_2^{KV} \left[ \varepsilon\_{n+1} - \varepsilon\_n + \mathcal{H}^{\mathcal{S}\_2} \varepsilon \right]. \tag{30}
$$

**Fractional Maxwell Model**: Evaluating both sides of (10) at *t* = *tn*+1, we obtain:

$$
\sigma\_{n+1} + \frac{\mathbb{E}\_2}{\mathbb{E}\_1} \, ^C\_0 \mathcal{D}\_t^{\theta\_2 - \beta\_1} \sigma(t) \Big|\_{t=t\_{n+1}} = \mathbb{E}\_2 \, ^C\_0 \mathcal{D}\_t^{\theta\_2} \varepsilon(t) \Big|\_{t=t\_{n+1}} \tag{31}
$$

in which applying (25) for the fractional derivatives of strains and stresses leads to:

$$\sigma\_{n+1} = \frac{\mathbb{C}\_1^M \left[ \varepsilon\_{n+1} - \varepsilon\_n + \mathcal{H}^{\beta\_2} \varepsilon \right] + \mathbb{C}\_2^M \left[ \sigma\_n - \mathcal{H}^{\beta\_2 - \beta\_1} \sigma \right]}{1 + \mathbb{C}\_2^M} \,\_{\prime} \tag{32}$$

with the emergence of a stress history term <sup>H</sup>*β*2−*β*1*σ*.

**Fractional Kelvin–Zener Model**: Evaluating both sides of (12) at *t* = *tn*+1, we obtain:

$$\begin{split} \left. \sigma\_{n+1} + \frac{\mathbb{E}\_{2}}{\mathbb{E}\_{1}} \mathop{\mathbb{C}}\_{0} \mathcal{D}\_{t}^{\mathbb{f}\_{2}-\mathbb{f}\_{1}} \boldsymbol{\sigma}(t) \right|\_{t=t\_{n+1}} &= \mathbb{E}\_{2} \mathop{\mathbb{C}}\_{0} \mathcal{D}\_{t}^{\mathbb{f}\_{2}} \boldsymbol{\varepsilon}(t) \big|\_{t=t\_{n+1}} + \mathbb{E}\_{3} \mathop{\mathbb{C}}\_{0} \mathcal{D}\_{t}^{\mathbb{f}\_{3}} \boldsymbol{\varepsilon}(t) \big|\_{t=t\_{n+1}} \\ + \frac{\mathbb{E}\_{2} \mathbb{E}\_{3}}{\mathbb{E}\_{1}} \mathop{\mathbb{C}}\_{0} \mathcal{D}\_{t}^{\mathbb{f}\_{2}+\mathbb{f}\_{3}-\mathbb{f}\_{1}} \boldsymbol{\varepsilon}(t) \big|\_{t=t\_{n+1}} \end{split}$$

which after applying (25) for the fractional derivatives of strains and stresses leads to:

$$\begin{split} \sigma\_{n+1} &= (1 + \mathsf{C}\_{4}^{KZ})^{-1} \Big[ \mathsf{C}\_{1}^{KZ} \Big( \Delta \varepsilon\_{n+1} + \mathsf{\mathcal{H}}^{\beta\_{2}} \varepsilon \Big) + \mathsf{C}\_{2}^{KZ} \Big( \Delta \varepsilon\_{n+1} + \mathsf{\mathcal{H}}^{\beta\_{3}} \varepsilon \Big) \\ &+ \mathsf{C}\_{3}^{KZ} \Big( \Delta \varepsilon\_{n+1} + \mathsf{\mathcal{H}}^{\beta\_{2} + \beta\_{3} - \beta\_{1}} \varepsilon \Big) + \mathsf{C}\_{4}^{KZ} \Big( \sigma\_{n} - \mathsf{\mathcal{H}}^{\beta\_{2} - \beta\_{1}} \sigma \Big) \Big] \end{split} \tag{33}$$

with Δ*εn*+<sup>1</sup> = *εn*+<sup>1</sup> − *εn*.

**Fractional Poynting–Thomson Model**: Finally, we evaluate both sides of (14) and obtain:

$$\begin{split} \left. \sigma\_{n+1} + \frac{\mathbb{E}\_{1}}{\mathbb{E}\_{3}} \mathop{\mathbb{C}}\_{0} \mathcal{D}\_{t}^{\beta\_{1} - \beta\_{3}} \sigma(t) \right|\_{t = t\_{n+1}} + \frac{\mathbb{E}\_{2}}{\mathbb{E}\_{3}} \mathop{\mathbb{C}}\_{0} \mathcal{D}\_{t}^{\beta\_{2} - \beta\_{3}} \sigma(t) \big|\_{t = t\_{n+1}} &= \mathbb{E}\_{1} \mathop{\mathbb{C}}\_{0} \mathcal{D}\_{t}^{\beta\_{1}} \varepsilon(t) \big|\_{t = t\_{n+1}} \\ + \mathbb{E}\_{2} \mathop{\mathbb{C}}\_{0} \mathbb{D}\_{t}^{\beta\_{2}} \varepsilon(t) \big|\_{t = t\_{n+1}} \end{split}$$

which after applying (25) for the fractional derivatives of strains and stresses, leads to:

$$\sigma\_{n+1} = \left(1 + \mathbb{C}\_3^{PT} + \mathbb{C}\_4^{PT}\right)^{-1} \left[\mathbb{C}\_1 \left(\Delta\varepsilon\_{n+1} + \mathcal{H}^{\mathbb{S}\_1}\varepsilon\right) + \mathbb{C}\_2^{PT} \left(\Delta\varepsilon\_{n+1} + \mathcal{H}^{\mathbb{S}\_2}\varepsilon\right)\right. \\ \left. + \mathbb{C}\_3^{PT} \left(\sigma\_{n+1} + \mathcal{H}^{\mathbb{S}\_1 - \mathbb{S}\_3}\sigma\right)\right]. \tag{34}$$

$$+ \mathbb{C}\_3^{PT} \left(\sigma\_{n+1} + \mathcal{H}^{\mathbb{S}\_1 - \mathbb{S}\_3}\sigma\right) + \mathbb{C}\_4^{PT} \left(\sigma\_n - \mathcal{H}^{\mathbb{S}\_2 - \mathbb{S}\_3}\sigma\right)\right]. \tag{35}$$

**Fractional Quasi-Linear Viscoelastic Model:** The discretization for the FQLV model (18) has a slightly different development than the preceding models. It involves a slight modification of the fully implicit L1 difference approach by a trapezoidal rule, taken on the exponential factor. More specifically, we evaluate the FQLV operator as:

$$\sigma\_{n+1} = \frac{EAB}{\Gamma(1-\beta)} \sum\_{k=0}^{n} \int\_{t\_k}^{t\_{k+1}} (t\_{n+1} - s)^{-\beta} \exp(B\varepsilon\_{k+\frac{1}{2}}) \left(\frac{\varepsilon\_{k+1} - \varepsilon\_k}{\Delta t}\right) ds \tag{35}$$

with *εi*<sup>+</sup> <sup>1</sup> 2 = (*ε<sup>i</sup>* + *εi*+1)/2. Following similar steps as in [33], we obtain the following discretized stresses at *t* = *tn*+<sup>1</sup> for the FQLV model:

$$
\sigma\_{n+1} = \mathcal{C}\_1^{QLV} \left[ \exp(B\varepsilon\_{n+\frac{1}{2}}) (\varepsilon\_{n+1} - \varepsilon\_n) + \mathcal{H}^a \left( \varepsilon, \frac{\partial \sigma^c}{\partial \varepsilon} \right) \right] \tag{36}
$$

with constant *CQLV* <sup>1</sup> = *EAB*/(Δ*t <sup>α</sup>*Γ(<sup>2</sup> <sup>−</sup> *<sup>α</sup>*)). The discretized history load in this case is given by:

$$\mathcal{H}^{a}\left(\varepsilon,\frac{\partial\sigma^{c}}{\partial\varepsilon}\right) = \sum\_{k=1}^{n} \exp(B\varepsilon\_{n-k+\frac{1}{2}}) (\varepsilon\_{n-k+1} - \varepsilon\_{n-k}) b\_{k} \tag{37}$$

with weights *bk* = (*<sup>k</sup>* <sup>+</sup> <sup>1</sup>)1−*<sup>α</sup>* <sup>−</sup> *<sup>k</sup>*1−*α*. Since the trapezoid approximation of the strains in the exponential term are second-order accurate, the overall accuracy of the viscoelastic models is still bounded by the native L1-difference approach and, therefore, should be of O(Δ*t* <sup>2</sup>−*α*).

**Remark 1.** *We note that except for the FQLV model, any of the aforementioned discretizations for the linear models can recover the existing classical counterparts by properly setting β<sup>i</sup>* → 0 *or β<sup>i</sup>* → 1*. In these cases, to achieve a comparable performance to the integer-order models, history terms can be selectively disregarded, and the corresponding discretization constants can be adjusted to their integer-order counterparts.*

#### *5.2. Time-Fractional Integration of Visco-Plasticity*

We start with the discretization of internal variables. Following [14], we assume a strain-driven process with known total strains *εn*+<sup>1</sup> at time *tn*+1. The strain decomposition becomes:

$$
\varepsilon\_{n+1} = \varepsilon\_{n+1}^{\text{rc}} + \varepsilon\_{n+1}^{\text{vp}}.\tag{38}
$$

The flow rule (22) is discretized through a first-order backward-Euler approach, which yields:

$$
\varepsilon\_{n+1}^{vp} = \varepsilon\_n^{vp} + \text{sign}(\sigma\_{n+1})\Delta\gamma\_{n+1} \tag{39}
$$

with Δ*γn*+<sup>1</sup> = *γn*+<sup>1</sup> − *γ<sup>n</sup>* representing the plastic slip increment in the interval [*tn*, *tn*+1]. Similarly, the discretization of the hardening law (23) is given by

$$
\alpha\_{n+1} = \alpha\_n + \Delta \gamma\_{n+1}.\tag{40}
$$

Evaluating the yield function (20) at *tn*+<sup>1</sup> and employing discretization (25) for the hardening variable, we obtain:

$$\begin{split} f\_{n+1} &= \left| \sigma\_{n+1} \right| - \left[ \sigma^Y + \mathbb{K}\_0^C \mathcal{D}\_t^{\mathcal{G}\_K}(\boldsymbol{a}) \big|\_{t=t\_{n+1}} + H \boldsymbol{a}\_{n+1} \right] \\ &= \left| \sigma\_{n+1} \right| - \left[ \sigma^Y + \mathbb{K}^\* \left( \boldsymbol{a}\_{n+1} - \boldsymbol{a}\_n + \mathcal{H}^{\mathcal{G}\_K} \boldsymbol{a} \right) + H \boldsymbol{a}\_{n+1} \right] \end{split} \tag{41}$$

with <sup>K</sup><sup>∗</sup> <sup>=</sup> <sup>K</sup>/(Δ*tβ<sup>K</sup>* <sup>Γ</sup>(<sup>2</sup> <sup>−</sup> *<sup>β</sup>K*)).

The next step is to define the trial states for the stress and yield functions, which is the core idea to define the viscoelastic prediction phase, and the correction step after solving the internal visco-plastic variables. Therefore, we freeze the internal variables for the prediction step at *tn*+1. Accordingly, the trial visco-plastic strains and hardening are given by:

$$
\varepsilon\_{n+1}^{vp^{trial}} = \varepsilon\_n^{vp}, \quad \alpha\_{n+1}^{trial} = \alpha\_n. \tag{42}
$$

In this token, the trial yield function is given by setting the above relationship for the hardening variable into (41) to obtain:

$$f\_{n+1}^{trial} = |\sigma\_{n+1}^{trial}| - \left[\sigma^Y + \mathbb{K}^\* \left(\mathcal{H}^{\beta\_K} \mathfrak{a}\right) + H\alpha\_n\right].\tag{43}$$

In order to complete the return-mapping procedure, we need an explicit relationship between the stresses *σn*+<sup>1</sup> in terms of the known total strains *εn*<sup>+</sup>1. To achieve this, we solve for the plastic slip Δ*γ* using a discrete consistency condition *fn*+<sup>1</sup> = 0. We start with the trial stresses for each presented fractional viscoelastic model by substituting the visco-plastic trial strain (42) and (38) into (28)–(36), where we obtain the following for each discretized model:

*Scott-Blair:*

$$
\sigma\_{n+1}^{trial} = \mathbb{C}\_1^{SB} \left[ \varepsilon\_{n+1} - \varepsilon\_n + \mathcal{H}^{\mathcal{S}\_1} (\varepsilon - \varepsilon^{vp}) \right] \tag{44}
$$

*Fractional Kelvin–Voigt:*

$$\sigma\_{n+1}^{trial} = \mathbb{C}\_1^{KV} \left[ \varepsilon\_{n+1} - \varepsilon\_n + \mathcal{H}^{\mathbb{R}\_1} (\varepsilon - \varepsilon^{vp}) \right] + \mathbb{C}\_2^{KV} \left[ \varepsilon\_{n+1} - \varepsilon\_n + \mathcal{H}^{\mathbb{R}\_2} (\varepsilon - \varepsilon^{vp}) \right] \tag{45}$$

*Fractional Maxwell:*

$$\sigma\_{n+1}^{trial} = \frac{\mathbb{C}\_1^M \left[ \varepsilon\_{n+1} - \varepsilon\_n + \mathcal{H}^{\mathbb{R}\_2} (\varepsilon - \varepsilon^{vp}) \right] + \mathbb{C}\_2^M \left[ \sigma\_n - \mathcal{H}^{\mathbb{R}\_2 - \beta\_1} \sigma \right]}{1 + \mathbb{C}\_2^M} \tag{46}$$

*Fractional Kelvin–Zener:*

$$\begin{split} \sigma\_{n+1}^{trial} &= (1 + \mathsf{C}\_{4}^{KZ})^{-1} \left[ \mathsf{C}\_{1}^{KZ} \left( \Delta \varepsilon\_{n+1} + \mathcal{H}^{\beta\_{2}} (\varepsilon - \varepsilon^{\upsilon p}) \right) + \mathsf{C}\_{2}^{KZ} \left( \Delta \varepsilon\_{n+1} + \mathcal{H}^{\beta\_{3}} (\varepsilon - \varepsilon^{\upsilon p}) \right) \right. \\ &\left. + \mathsf{C}\_{3}^{KZ} \left( \Delta \varepsilon\_{n+1} + \mathcal{H}^{\beta\_{2} + \beta\_{3} - \beta\_{1}} (\varepsilon - \varepsilon^{\upsilon p}) \right) + \mathsf{C}\_{4}^{KZ} \left( \sigma\_{n} - \mathcal{H}^{\beta\_{2} - \beta\_{1}} \sigma \right) \right] \end{split} \tag{47}$$

*Fractional Poynting–Thomson:*

$$
\sigma\_{n+1}^{\text{trial}} = (1 + \mathbb{C}\_3^{\text{PT}} + \mathbb{C}\_4^{\text{PT}})^{-1} \left[ \mathbb{C}\_1^{\text{PT}} \left( \Delta \varepsilon\_{n+1} + \mathcal{H}^{\text{f}\_1} (\varepsilon - \varepsilon^{\text{vp}}) \right) + \mathbb{C}\_2^{\text{PT}} \left( \Delta \varepsilon\_{n+1} + \mathcal{H}^{\text{f}\_2} (\varepsilon - \varepsilon^{\text{vp}}) \right) \right]
$$

$$
+ \mathbb{C}\_3^{\text{PT}} \left( \sigma\_{n+1} + \mathcal{H}^{\text{f}\_1 - \beta\_3} \sigma \right) + \mathbb{C}\_4^{\text{PT}} \left( \sigma\_n - \mathcal{H}^{\text{f}\_2 - \beta\_3} \sigma \right) \right] \tag{48}
$$

*Fractional Quasi-Linear Viscoelastic Model:*

For this model, we follow a similar procedure of substituting the viscoelastic strains into (36), however, we evaluate the exponential term explicitly in time for all stages of the return-mapping algorithm. Therefore, the corresponding trial state becomes:

$$
\sigma\_{n+1}^{trial} = \mathbb{C}\_1^{QLV} \left[ \exp(B(\varepsilon\_n - \varepsilon\_n^{vp})) (\varepsilon\_{n+1} - \varepsilon\_n) + \mathcal{H}^a \left( \varepsilon - \varepsilon^{vp}, \frac{\partial \sigma^\varepsilon}{\partial \varepsilon} \right) \right]. \tag{49}
$$

#### *5.3. Generalized Fractional Return-Mapping Algorithm (Algorithm 1)*

From the aforementioned trial states, each discretized viscoelastic constitutive laws (28)–(36) and recalling (39), one can show the following stress correction onto the yield surface:

$$
\sigma\_{n+1} = \sigma\_{n+1}^{trial} - \text{sign}(\sigma\_{n+1}^{trial}) \mathbb{C}\_{RM}^{\text{rc}}(\mathbb{E}\_{\prime} \text{ } \Delta t \text{ } \varepsilon) \Delta \gamma\_{n+1 \prime} \tag{50}
$$

where all discretized aforementioned viscoelastic models change the return-mapping procedure by a scaling factor *Cve RM*(*C*, *εn*, *ε vp <sup>n</sup>* ) <sup>∈</sup> <sup>R</sup><sup>+</sup> acting on the Lagrange multiplier <sup>Δ</sup>*γ*, which is given by the following for each model:

$$\mathbf{C}\_{RM}^{\text{ref}} = \begin{cases} \mathbf{C}\_{1}^{SB} & (\text{Soft-Blair}) \\ \mathbf{C}\_{1}^{KV} + \mathbf{C}\_{2}^{KV} & (\text{Fractional Kelvin-Voigt}) \\ \mathbf{C}\_{1}^{M}/(1 + \mathbf{C}\_{2}^{M}) & (\text{Fractional Maxwell}) \\ \left(\mathbf{C}\_{1}^{KZ} + \mathbf{C}\_{2}^{KZ} + \mathbf{C}\_{3}^{KZ}\right)/(1 + \mathbf{C}\_{4}^{KZ}) & (\text{Fractional Kelvin-Zener}) \\ \left(\mathbf{C}\_{1}^{PT} + \mathbf{C}\_{2}^{PT}\right)/(1 + \mathbf{C}\_{3}^{PT} + \mathbf{C}\_{4}^{PT}) & (\text{Fractional Population-Thomas}) \\ \mathbf{C}\_{1}^{QLV} \exp(B(\varepsilon\_{n} - \varepsilon\_{n}^{vp})) & (\text{Fractional Quasis-Linear-Viscolastic}) \end{cases} \tag{51}$$

We show the derivation of (50) for the fractional Kelvin–Zener model in Appendix C, from which the Scott-Blair, fractional Maxwell, and fractional Kelvin–Voigt models can be directly recovered; note that the derivation for the fractional Poynting–Thomson and quasi-linear viscoelasticity follow similarly in a straightforward fashion. Substituting the updated stresses (50) into the discrete yield function (41) and recalling (43), we obtain:

$$f\_{n+1} = f\_{n+1}^{trial} - (\mathbb{C}\_{RM}^{\text{rc}} + \mathbb{K}^\* + H)\Delta\gamma. \tag{52}$$

Enforcing the discrete yield condition *fn*+<sup>1</sup> = 0, we obtain the solution for the discrete plastic slip:

$$
\Delta \gamma\_{n+1} = \frac{f\_{n+1}^{trial}}{\mathbb{C}\_{RM}^{\text{rc}} + \mathbb{K}^\* + H} \Bigg| \tag{53}
$$

### **Algorithm 1** Fractional return-mapping algorithm.

1: Database for *ε*, *εvp*, *σ*, *α*, and total strain *εn*+1. 2: *ε vptrial <sup>n</sup>*+<sup>1</sup> = *ε vp <sup>n</sup>* , *αtrial <sup>n</sup>*+<sup>1</sup> = *α<sup>n</sup>* 3: Compute *σtrial <sup>n</sup>*+<sup>1</sup> from (28)–(36) according to the selected fractional viscoelastic model. 4: *f trial <sup>n</sup>*+<sup>1</sup> <sup>=</sup> <sup>|</sup>*σtrial <sup>n</sup>*+<sup>1</sup> | − *σ<sup>Y</sup>* + K<sup>∗</sup> <sup>H</sup>*β<sup>K</sup> <sup>α</sup>* + *Hα<sup>n</sup>* 5: **if** *f trial <sup>n</sup>*+<sup>1</sup> ≤ 0 **then** 6: *ε vp <sup>n</sup>*+<sup>1</sup> = *ε vp <sup>n</sup>* , *αn*+<sup>1</sup> = *αn*, *σn*+<sup>1</sup> = *σtrial <sup>n</sup>*+<sup>1</sup> . 7: **else** 8: Return-Mapping: 9: Compute *Cve RM* from (51) according to the selected fractional viscoelastic model. 10: Δ*γn*+<sup>1</sup> = *f trial <sup>n</sup>*+<sup>1</sup> /(*Cve RM* + *K*<sup>∗</sup> + *H*) 11: *σn*+<sup>1</sup> = *σtrial <sup>n</sup>*+<sup>1</sup> <sup>−</sup> sign(*σtrial <sup>n</sup>*+<sup>1</sup> )*Cve RM*Δ*γ* 12: *ε vp <sup>n</sup>*+<sup>1</sup> = *ε vp <sup>n</sup>* + sign(*τn*+1)Δ*γ* 13: *αn*+<sup>1</sup> = *α<sup>n</sup>* + Δ*γ* 14: **end if**

Comparison of the Return-Mapping Algorithm to the Existing Approaches

In [14], trial states were defined prior to the discretization of fractional operators, and the corresponding trial variables were taken as continuous functions of time, therefore making the return-mapping procedure "semi-discrete." Let the quantities with bars, (¯·), be the corresponding solutions for the procedure developed in [14]. For the SB viscoelastic case, one has the following trial stresses at *t* = *tn*+1:

$$\bar{\sigma}\_{n+1}^{\text{trial}} = \mathbb{E}\, ^C\_0 \mathcal{D}\_t^{\mathcal{G}\_E} (\varepsilon - \mathbb{E}^{vp^{trial}}) \big|\_{t=t\_{n+1}} \tag{54}$$

in which, after employing the discretized plastic flow rule, the following relationship between the corrected and trial stresses is obtained:

$$
\bar{\sigma}\_{n+1} = \bar{\sigma}\_{n+1}^{\text{trial}} - \mathbb{E}\,\text{sign}(\bar{\sigma}\_{n+1}) \, ^C\_0 \mathcal{D}\_t^{\mathcal{G}\_E}(\Delta \gamma) \Big|\_{t=t\_{n+1}}.\tag{55}
$$

This equation can be explicitly inserted into the discrete yield function to solve for the plastic slip rate. While such a procedure is straightforward for SB and FKV viscoelastic elements, it is non-trivial for the serial combinations such as the FM, FKZ, and FPT models. For instance, if we follow the same procedure for the FM model, we obtain:

$$\left. \mathcal{D}\_{n+1}^{trial} + \frac{\mathbb{E}\_2}{\mathbb{E}\_1} \, ^{\mathbb{C}}\_0 \mathcal{D}\_t^{\mathcal{G}\_2 - \mathcal{G}\_1} (\mathcal{O}^{trial}) \right|\_{t=t\_{n+1}} = \mathbb{E}\_2 \, ^{\mathbb{C}}\_0 \mathcal{D}\_t^{\mathcal{G}\_2} (\varepsilon - \varepsilon^{wp^{trial}}) \vert\_{t=t\_{n+1} \prime} \tag{56}$$

which yields the following relationship between *σ*¯ and *σ*¯*trial*:

$$\mathbb{E}\left(\left.\bar{\sigma}\_{n+1} - \bar{\sigma}\_{n+1}^{trial}\right) + \frac{\mathbb{E}\_2}{\mathbb{E}\_1} \, ^C\_0 \mathcal{D}\_t^{\mathbb{S}\_2 - \mathbb{S}\_1} (\bar{\sigma} - \bar{\sigma}^{trial})\Big|\_{t=t\_{n+1}} = -\mathbb{E}\_2 \, ^C\_0 \mathcal{D}\_t^{\mathbb{S}\_2} (\varepsilon - \mathbb{E}^{\mathbb{s}p^{trial}})\Big|\_{t=t\_{n+1}}.\tag{57}$$

Except for the SB case, a fractional viscoelastic model involving a serial combination of SB elements cannot be incorporated into the yield function in a differential form unless a full discretization is performed at this stage. This happens since the discretized yield function (41) requires a closed description of *σn*<sup>+</sup>1, needing an equivalent Boltzmann representation for such models, which is impractical due to complex forms of relaxation kernels. Therefore, our approach in this work already carries the trial states with fully discretized fractional operators, which closely and completely resembles classical elasto-plastic approaches.

Regarding the obtained discretizations in this work, we note that the plastic slip (53) assumes a simple form similar to the rate-independent elasto-plasticity. As discussed above, in the return-mapping procedure developed in [14], the trial states and plastic slip were assumed to have memory in the discretization procedure; therefore, a fractional relaxation equation in the following form was obtained:

$$
\Delta\bar{\gamma}\_{n+1} = \frac{\mathbb{E}^\* \left(\Delta\bar{\gamma}\_n - \mathcal{H}^{\mathcal{S}\_E} \Delta\bar{\gamma}\right) + \mathbb{K}^\* \left(\Delta\bar{\gamma}\_n - \mathcal{H}^{\mathcal{S}\_K} \Delta\bar{\gamma}\right) + \bar{f}\_{n+1}^{\text{trial}}}{\mathbb{E}^\* + \mathbb{K}^\* + H}. \tag{58}
$$

Furthermore, we observe that the obtained plastic slip discretization in this work has two less history terms to be evaluated. Although this does not influence the computational complexity of the original scheme, we show in the numerical examples that this fact still leads to about 50% less CPU time. Regarding the difference in the stress solutions, let *t* = *tp* be the time step of onset of plasticity for the first time. Therefore, we have the following estimate:

$$|\sigma\_{p+1} - \sigma\_{p+1}| = \frac{\mathbb{E}^\*}{\mathbb{E}^\* + \mathbb{K}^\* + H} \left[ \mathbb{K}^\* \left( \mathcal{H}^{\mathbb{\beta}\_E} \Delta \tilde{\gamma} - \mathcal{H}^{\mathbb{\beta}\_K} \Delta \tilde{\gamma} \right) - H \left( \Delta \tilde{\gamma}\_p - \mathcal{H}^{\mathbb{\beta}\_E} \Delta \tilde{\gamma} \right) \right], \tag{59}$$

which shows that at such stage both discretizations coincide when *β<sup>E</sup>* = *β<sup>K</sup>* and *H* = 0. In the following Section, we verify such an estimate by obtaining an analytical solution with the aid of Proposition 1.

#### **6. Numerical Tests**

We present three convergence examples with different loading conditions to verify the employed fractional viscoelastic models, the validity of the new fractional visco-plastic return-mapping algorithm, and the full visco-elasto-plastic response of the models. For the convergence analyses, let **u**<sup>∗</sup> and **u***<sup>δ</sup>* be, respectively, the reference and approximate solutions in Ω = (0, *T*], for a specific time-step size Δ*t*. We define the following relative error measures:

$$\text{errr}\_N(\Delta t) = \frac{|\boldsymbol{u}\_N^\* - \boldsymbol{u}\_N^\delta|}{|\boldsymbol{u}\_N^\*|}, \text{err}(\Delta t) = \frac{||\mathbf{u}^\* - \mathbf{u}^\delta||\_{L^2(\Omega)}}{||\mathbf{u}^\*||\_{L^2(\Omega)}}, \text{Order} = \log\_2\left[\frac{\text{err}(\Delta t)}{\text{err}(\Delta t/2)}\right].\tag{60}$$

We consider homogeneous initial conditions for all model variables in all cases. The presented algorithms were implemented in MATLAB R2020b and were run in a system with Intel Core i7-8850H CPU with 2.60 GHz, 32 GB RAM and MacOS 11.5 operating system.

**Example 1** (Convergence of fractional viscoelastic algorithms)**.** *We perform a convergence study of the fractional viscoelastic component of our framework under the stress relaxation and monotone loading experiments. For this example, we set* (E1,E2,E3)=(1, 1, 1) *and* (*β*1, *β*2, *β*3) = (0.3, 0.7, 0.1) *for fractional linear viscoelastic models, ensuring all fractional derivatives are taken with an "equivalent order" β* ∈ (0, 1)*, i.e., the sum of fractional orders arising in the fractional derivatives of each linear viscoelastic model. For the FQLV model, we set E* = 1*, A* = 1*, and B* = 1*, and varying fractional order β.*

*For the stress relaxation test, we impose a step strain ε*(*t*) = *H*(*t*)*ε*<sup>0</sup> *with ε*<sup>0</sup> = 1 *for T* = 1000 [*s*]*, where H*(*t*) *denotes the Heaviside step function. We compare the obtained solutions at t* = *T for the SB, FKV, FM, and FKZ models to their corresponding relaxation functions* (7)*,* (9)*,* (11) *and* (13)*. The FPT and QLV models are not analyzed in this step since their timedependent stress relaxation functions are not readily available, and they are instead analyzed under the monotone strains. The obtained results are illustrated in Figure 3a, where an expected linear convergence behavior is obtained for all models when the error is evaluated at the end point, given the non-smooth nature of the applied step strain in the stress relaxation solution.*

**Figure 3.** Convergence analysis for the fractional viscoelastic models with known analytical solutions: (**a**) a stress relaxation test with non-smooth step-strains and material parameters (E1,E2,E3)=(1, 1, 1) and (*β*1, *β*2, *β*3)=(0.3, 0.7, 0.1) yielding first-order convergence; (**b**) convergence for the FQLV model with a fabricated solution of linearly increasing strains and material properties (*E*, *A*, *B*)=(1, 1, 1) and varying *β*. The slopes of the error curves are *q* ≈ 2 − *β*.

*For the monotone strain case, we set T* = 1 *and fabricate a solution for strains in the form ε*(*t*) = *εT*(*t*/*T*) *with the total applied strain fixed at ε<sup>f</sup>* = 1*. Since analytical solutions for all fractional viscoelastic models are difficult to obtain, we compute a reference solution for each model with* Δ*t* = 2−<sup>17</sup> [*s*]*. Particularly for the FQLV model, we utilize the fabricated strain function ε*(*t*) *to obtain the following analytical stress solution:*

$$\sigma^{QLV^\*} (t) = EAB^\beta \exp(Bt) \left[ 1 - \frac{\Gamma(1 - \beta\_\prime Bt)}{\Gamma(1 - \beta)} \right],\tag{61}$$

*where* Γ(·, ·) *denotes the upper incomplete gamma function. The convergence results for all fractional viscoelastic models with respect to the reference numerical solution are presented in Figure 4, while the results for the FQLV model with the analytical solution are illustrated in Figure 3b. We observe for both cases that the accuracy of the implemented and developed schemes is of order* O(Δ*t* <sup>2</sup>−*β*)*. We note that although we have employed a trapezoid rule for the exponential term in the FQLV model, we do not obtain a second-order convergence, since the accuracy is limited by the L1 approach. The difference in error slopes among models in Figure 4 is due to the highest fractional order assigned to each model. For the SB and FQLV models, the fractional order is set as β* = 0.3*, and therefore, the observed slope is q* ≈ 1.7*. For all remaining models and choice of fractional orders, the error slopes are determined by the fractional derivative of the highest order, which is β*<sup>2</sup> = 0.7 *in this example, yielding q* ≈ 1.3*.*

**Figure 4.** Convergence analysis for all fractional viscoelastic models with (E1,E2,E3)=(1, 1, 1) and (*β*1, *β*2, *β*3)=(0.3, 0.7, 0.1). A cubic strain function was employed with a reference solution using the time step size <sup>Δ</sup>*<sup>t</sup>* <sup>=</sup> <sup>2</sup>−17. A monotone loading test with the convergence rate of *<sup>q</sup>* <sup>≈</sup> 1.3 was applied for all models.

**Example 2** (Convergence of fractional visco-plastic algorithms)**.** *The purpose of this example is to demonstrate the conditions where the presented plastic slip discretization* (53)*, the form* (58) *from [14] and their associated return-mapping algorithms are equivalent, and also provide a numerical estimate for their difference when such conditions are not satisfied. For this purpose, we test a monotone load where an analytical solution is available and a case with a cyclic load under high strain rates. For both cases, we set an SB viscoelastic part with* E = 50 [*Pa*.*sβ<sup>E</sup>* ]*, and* K = 5 [*Pa*.*sβ<sup>K</sup>* ]*.*

*For the monotone strain case, we start with a fabricated solution for strains in the form ε*(*t*) = *At*<sup>3</sup> *with A* = *εT*/*T*<sup>3</sup> [*s*−3]*. Here, ε<sup>f</sup> denotes the total applied stress, and T represents the final simulation time. Utilizing the result of Proposition 1 and setting β<sup>E</sup>* = *β<sup>K</sup>* = 0 *and σ<sup>Y</sup>* = *H* = 0*, we obtain the following analytical solution for stresses:*

$$
\sigma^\*(t) = \frac{6A \to \mathbb{K}}{\mathbb{E} + \mathbb{K}} \frac{t^{3-\beta\_E}}{\Gamma(4-\beta\_E)}.\tag{62}
$$

*We should note that the proposed fabricated solution ensures no internal variable is a linear function and, therefore, not computed exactly by the L1 discretization. We set ε<sup>T</sup>* = 1 *and T* = 1 [*s*]*, and therefore, A* = 1 [*s*−3]*. Table 1 presents the obtained convergence results for the fabricated solution* (62) *for both return-mapping algorithms and under the same fractional-orders β<sup>E</sup> and βK. We observe that the errors coincide for this particular case, while the accuracy of order* O(Δ*t* <sup>2</sup>−*β*) *of the L1 approach is also achieved. The computational times are illustrated in Figure 5, where the developed fractional return-mapping approach, when using an SB viscoelastic element, is about* 50% *faster than the original return-mapping approach from [14] since about half of the history terms need to be computed.*

**Table 1.** Convergence behavior for the return-mapping Algorithm 1 obtained in this work and the original approach from [14] for an FVEP device with an SB element.


**Figure 5.** CPU times for the developed fractional return-mapping algorithm and the original one [14] for an SB viscoelastic part. The black line has slope *q* = 2.

*Similar results are obtained for the monotone loading condition; however, this is not the case under general loadings. To demonstrate the difference between the visco-elasto-plastic discretization σn*+<sup>1</sup> *developed in this work and σ*¯*n*+<sup>1</sup> *from [14], we take the latter as a reference solution with* Δ*t* = 2−<sup>19</sup> [*s*]*, and T* = 1 [*s*]*. We also consider σ<sup>Y</sup>* = 10 [*Pa*]*, β<sup>E</sup>* = 0.3*, and β<sup>K</sup>* = 0.7 *with the same pseudo-constants as in the previous test case. A constant rate loading/unloading cyclic strain test of the following form is employed:*

$$\varepsilon(t) = \frac{2\varepsilon\_A}{\pi} \arcsin(\sin(2\pi\omega t)),\tag{63}$$

*where we consider a strain amplitude ε<sup>A</sup>* = 0.25 *and two strain frequencies of ω* = 1 [*Hz*]*, and ω* = 60 [*Hz*]*. The difference between both approaches is illustrated in Figure 6. Here, higher frequencies result in higher strain rates and consequently a significant plastic strain history even after a number of hysteresis cycles. The obtained results confirm the estimates from* (59)*, which is already valid at the onset of plasticity. Furthermore, we observe that a tenfold increase in strain rates approximately leads to a tenfold increase in the difference between both algorithms.*

**Figure 6.** Comparison between the presented return-mapping algorithm and the reference approach from [14] under low- and high-frequency loading.

**Example 3** (Convergence of fractional visco-elasto-plasticity)**.** *Finally, we perform a verification on the entire fractional visco-elasto-plastic framework under cyclic strain. Since no fabricated solutions are available, we employ reference solutions with time step size* Δ*t* = 2−<sup>18</sup> [*s*]*. Let T* = 1 [*s*] *with the same applied strains* (63) *as in the previous example. The viscoelastic material parameters are set to* (E1,E2,E3)=(50, 50, 50)*, and* (*β*1, *β*2, *β*3)=(0.3, 0.7, 0.1)*. In addition, the visco-plastic parameters are taken as* K = 5*, β<sup>K</sup>* = 0.7*, H* = 0*, and σ<sup>Y</sup>* = 1*. Figures 7 and 8 illustrates the obtained convergence results, where all models except the FKV one showed a convergence rate of order q* ≈ 1.3*, which is compatible with the employed L1 discretization scheme and given β*<sup>2</sup> = *β<sup>K</sup>* = 0.7*. The FKV model achieved linear asymptotic convergence for the considered example, which is the expected worst case scenario from the backward-Euler discretization of internal variables. We believe the difference in convergence behavior between the FKV model and the others could be due to the sharper response of the FKV model because of the stiffer rheology combined with the nonlinear loading/unloading response. This combination of effects could result in a lower solution regularity and therefore, a lower convergence rate.*

**Figure 7.** Convergence analysis for the fractional visco-elasto-plastic models under cyclic loads. Due to the particular choice of fractional orders (with *β*<sup>2</sup> = *β<sup>K</sup>* = 0.7 being dominant), we observed the convergence rate of *q* ≈ 1.3 for all models except for the FKV. In the latter case, we observed a linear convergence to the reference solution.

**Figure 8.** Visco-elasto-plastic reference solutions for the employed models for the first 30 loading cycles. We noticed a similar behavior for most models under the choice of material parameters except for the FPT and FKV models. The FKV particularly yielded a very stiff response due to the combination of high fractional order values and high strain rates.

#### **7. Conclusions**

We proposed a general return-mapping procedure for multiple power-law, timefractional visco-elasto-plastic materials. The developed framework provided a flexible way to integrate multiple known fractional viscoelastic models that are representative of soft materials rheology to power-law visco-plastic hardening and permanent strains. Furthermore, a nonlinear viscoelasticity, suitable for bio-tissues, was considered through a fractional quasi-linear Fung's model, which allowed the possibility of plasticity onset after substantial amounts of viscoelastic strains. The main features of the proposed framework are:


Regarding the computational costs, the framework is computationally tractable since it does not involve history calculations for the plastic slip and is therefore about 50% faster than the existing fractional frameworks, regardless of the employed numerical discretization for fractional derivatives. Extensions on fast numerical schemes of order O(*N* log *N*) for the employed time-fractional derivatives would be straightforward to implement. We also note that the thermodynamics of all models in the developed framework can be analyzed through the approach developed in [16].

The modeling framework developed here could be applied to simulate the self-similar structures and memory-dependent behavior in human bio-tissues [3,35,36]. The viscoelasto-plastic characteristics can be observed in different bio-tissues in the human body, specifically due to the process of aging. Aging results in the oxidation or loss of elastin, which leads to the loss of tissue elasticity such as in the vocal fold tissues [37,38]. Furthermore, in terms of multi-scale modeling, the lumped plastic behavior introduced here could potentially be coupled with existing discrete dislocation dynamics (DDD) models [39,40]. The models developed in this work uniquely qualify for simulating such characteristics, which will be undertaken in our future work.

**Author Contributions:** Conceptualization, J.L.S. and M.Z.; Methodology, J.L.S., M.N. and M.Z.; Investigation, M.N. and M.Z.; Resources, M.N. and M.Z.; Data curation, J.L.S.; Writing—original draft, J.L.S., M.N. and M.Z.; Visualization, J.L.S.; Supervision, M.N. and M.Z.; Project administration, M.N. and M.Z.; Funding acquisition, J.L.S. and M.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the ARO YIP Award (W911NF-19-1-0444), the NSF Award (DMS-1923201), the MURI/ARO Award (W911NF-15-1-0562), the AFOSR YIP Award (FA9550-17-1- 0150), and NIH NIDCD K01DC017751 and R21DC020003.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Appendix A. Proof of Proposition 1**

**Proof.** Similar to the derivation of the tangent elasto-plastic modulus in classical plasticity [32], we start by taking the time derivative of the yield function to enforce the persistency condition:

$$\begin{split} \dot{f}(\sigma, a) &= \frac{d}{dt} \left\{ |\sigma(t)| - \left[ \sigma^{\bar{Y}} + \mathbb{K}\_0^{\mathbb{C}} \mathcal{D}\_t^{\mathcal{G}\_K} a(t) + H a(t) \right] \right\} \\ &= \text{sign}(\sigma) \vartheta(t) - \left[ \mathbb{K} \frac{d}{dt} \mathbb{C} \mathcal{D}\_t^{\mathcal{G}\_K} a(t) + H a(t) \right]. \end{split} \tag{A1}$$

Using the SB stressd–strain relationship (6), we obtain:

$$\dot{f}(\sigma, a) = \text{sign}(\sigma) \mathbb{E} \left[ \frac{d}{dt} \prescript{\text{C}}{}{\mathcal{D}}\_{t}^{\mathcal{S}\_{E}} \varepsilon(t) - \frac{d}{dt} \prescript{\text{C}}{}{\mathcal{D}}\_{t}^{\mathcal{S}\_{E}} \varepsilon^{\text{up}}(t) \right] - \left[ \mathbb{K} \frac{d}{dt} \prescript{\text{C}}{}{\mathcal{D}}\_{t}^{\mathcal{S}\_{E}} a(t) + H\dot{a}(t) \right]. \tag{A2}$$

Employing definition (3) for the Caputo derivative, performing integration by parts and employing the Leibniz integral rule, we obtain:

$$\begin{split} \frac{d}{dt} \, ^C\_0 \mathcal{D}\_t^{\beta} u(t) &= \frac{1}{\Gamma(1-\beta)} \frac{d}{dt} \int\_0^t \frac{\dot{u}(s)}{(t-s)^{\beta}} ds \quad \text{(from(3))}\\ &= \frac{1}{\Gamma(1-\beta)} \frac{d}{dt} \left[ \dot{u}(s) \frac{(t-s)^{1-\beta}}{1-\beta} \Big|\_t^0 + \int\_0^t \frac{(t-s)^{1-\beta} \ddot{u}(s)}{1-\beta} ds \right] \\ &= \frac{\dot{u}(0)t^{-\beta}}{\Gamma(1-\beta)} + \frac{1}{\Gamma(1-\beta)} \int\_0^t \frac{\ddot{u}(s)}{(t-s)^{\beta}} ds \\ &= \frac{\dot{u}(0)t^{-\beta}}{\Gamma(1-\beta)} + \,\_0^C \mathcal{D}\_t^{\beta} \dot{u}(t). \end{split} \tag{A3}$$

Substituting (A3) into (A2), setting *γ*˙(0) = 0, and therefore *α*˙(0) = 0 from (23) and *ε*˙ *vp*(0) = 0 from (22), we obtain:

$$\dot{f}(\sigma, a) = \text{sign}(\sigma) \mathbb{E} \left[ \frac{\dot{\varepsilon}(0)t^{-\beta\_E}}{\Gamma(1-\beta\_E)} + {}\_0^C \mathcal{D}\_t^{\beta\_E} \dot{\varepsilon}(t) - {}\_0^C \mathcal{D}\_t^{\beta\_E} \dot{\varepsilon}^{vp}(t) \right] - \mathbb{K}\_0^C \mathcal{D}\_t^{\beta\_K} \dot{u}(t) - H\dot{u}(t). \tag{A4}$$

Finally, substituting (23) and (22) into (A4), and enforcing the persistency condition ˙ *f*(*σ*, *α*) = 0, we obtain:

$$\mathbb{E}^{\mathbb{C}}\_{0}\mathcal{D}^{\delta\_{E}}\_{t}\dot{\gamma}(t) + \mathbb{K}^{\mathbb{C}}\_{0}\mathcal{D}^{\delta\_{K}}\_{t}\dot{\gamma}(t) + H\dot{\gamma}(t) = \text{sign}(\sigma)\mathbb{E}\left[\frac{\boldsymbol{\varepsilon}(0)t^{-\beta\_{E}}}{\Gamma(1-\beta\_{E})} + ^{\mathbb{C}}\_{0}\mathcal{D}^{\delta\_{E}}\_{t}\dot{\varepsilon}(t)\right].\tag{A5}$$

**Appendix B. Discretization Constants and Terms for Fractional Viscoelastic Models** *Scott-Blair:*

$$\mathbf{C}\_{1}^{SB} = \frac{\mathbf{E}}{\Delta t^{\beta\_1} \Gamma(2-\beta\_1)}\tag{A6}$$

*Fractional Kelvin–Voigt:*

$$\mathbb{C}\_{1}^{KV} = \frac{\mathbb{E}\_{1}}{\Delta t^{\beta\_{1}} \Gamma(2-\beta\_{1})}, \quad \mathbb{C}\_{2}^{KV} = \frac{\mathbb{E}\_{2}}{\Delta t^{\beta\_{2}} \Gamma(2-\beta\_{2})} \tag{A7}$$

*Fractional Maxwell:*

$$\mathbb{C}\_1^M = \frac{\mathbb{E}\_2}{\Delta t^{\beta\_2} \Gamma(2-\beta\_2)}, \quad \mathbb{C}\_2^M = \frac{\mathbb{E}\_2/\mathbb{E}\_1}{\Delta t^{\beta\_2-\beta\_1} \Gamma(2-\beta\_2+\beta\_1)}\tag{A8}$$

*Fractional Kelvin–Zener:*

$$\mathbb{C}\_{1}^{KZ} = \frac{\mathbb{E}\_{2}}{\Delta t^{\beta\_{2}} \Gamma(2-\beta\_{2})}, \quad \mathbb{C}\_{2}^{KZ} = \frac{\mathbb{E}\_{3}}{\Delta t^{\beta\_{3}} \Gamma(2-\beta\_{3})} \tag{A9}$$

$$\mathbb{C}\_{3}^{KZ} = \frac{\mathbb{E}\_{2}\mathbb{E}\_{3}/\mathbb{E}\_{1}}{\Delta t^{\beta\_{2} + \beta\_{3} - \beta\_{1}}\Gamma(2 - \beta\_{1} - \beta\_{3} + \beta\_{2})}, \quad \mathbb{C}\_{4}^{KZ} = \frac{\mathbb{E}\_{2}/\mathbb{E}\_{1}}{\Delta t^{\beta\_{2} - \beta\_{1}}\Gamma(2 - \beta\_{2} + \beta\_{1})} \tag{A10}$$

*Fractional Poynting–Thomson:*

$$\mathbb{C}\_{1}^{PT} = \frac{\mathbb{E}\_{1}}{\Delta t^{\beta\_{1}} \Gamma(2-\beta\_{1})}, \quad \mathbb{C}\_{2}^{PT} = \frac{\mathbb{E}\_{2}}{\Delta t^{\beta\_{2}} \Gamma(2-\beta\_{2})} \tag{A11}$$

$$\mathbb{C}\_{3}^{PT} = \frac{\mathbb{E}\_{1}/\mathbb{E}\_{3}}{\Delta t^{\beta\_{1} - \beta\_{3}} \Gamma(2 - \beta\_{1} + \beta\_{3})}, \quad \mathbb{C}\_{4}^{PT} = \frac{\mathbb{E}\_{2}/\mathbb{E}\_{3}}{\Delta t^{\beta\_{2} - \beta\_{3}} \Gamma(2 - \beta\_{2} + \beta\_{3})} \tag{A12}$$

*Fractional Quasi-Linear viscoelastic:*

$$C\_1^{QLV} = \frac{EAB}{\Delta t^{\beta} \Gamma(2-\beta)}\tag{A13}$$

#### **Appendix C. Return-Mapping Derivation for the Fractional Kelvin–Zener Model**

Recalling the discretized FKV model (33) employed as the viscoelastic part of the visco-elasto-plastic model:

$$\begin{split} \sigma\_{n+1} &= (1 + \mathsf{C}\_{4}^{\mathrm{KZ}})^{-1} \Big[ \mathsf{C}\_{1}^{\mathrm{KZ}} \Big( \Delta \varepsilon\_{n+1}^{\mathrm{xe}} + \mathcal{H}^{\beta\_{2}} \varepsilon \Big) + \mathsf{C}\_{2}^{\mathrm{KZ}} \Big( \Delta \varepsilon\_{n+1}^{\mathrm{xe}} + \mathcal{H}^{\beta\_{3}} \varepsilon^{\mathrm{xe}} \Big) \\ &+ \mathsf{C}\_{3}^{\mathrm{KZ}} \Big( \Delta \varepsilon\_{n+1}^{\mathrm{xe}} + \mathcal{H}^{\beta\_{2} + \beta\_{3} - \beta\_{1}} \varepsilon^{\mathrm{xe}} \Big) + \mathsf{C}\_{4}^{\mathrm{KZ}} \Big( \sigma\_{n} - \mathcal{H}^{\beta\_{2} - \beta\_{1}} \sigma \Big) \Big], \end{split} \tag{A14}$$

where with the kinematic relationship (38) and the viscoplastic strain evolution (39), we note that Δ*εve <sup>n</sup>*+<sup>1</sup> = Δ*εn*+<sup>1</sup> − Δ*γn*+<sup>1</sup> sign(*σn*+1). Therefore, (A14) becomes:

$$
\begin{split}
\sigma\_{n+1} &= (1 + \mathsf{C}\_{4}^{KZ})^{-1} \Big[ \mathsf{C}\_{1}^{KZ} \Big( \Delta \varepsilon\_{n+1} + \mathsf{H}^{\beta\_{2}} \varepsilon \Big) + \mathsf{C}\_{2}^{KZ} \Big( \Delta \varepsilon\_{n+1} + \mathsf{H}^{\beta\_{3}} \varepsilon^{\mathrm{vc}} \Big) \\ &+ \mathsf{C}\_{3}^{KZ} \Big( \Delta \varepsilon\_{n+1} + \mathsf{H}^{\beta\_{2} + \beta\_{3} - \beta\_{1}} \varepsilon^{\mathrm{vc}} \Big) + \mathsf{C}\_{4}^{KZ} \Big( \sigma\_{n} - \mathsf{H}^{\beta\_{2} - \beta\_{1}} \sigma \Big) \\ &- \Big( \mathsf{C}\_{1}^{KZ} + \mathsf{C}\_{2}^{KZ} + \mathsf{C}\_{3}^{KZ} \Big) \Delta \gamma\_{n+1} \operatorname{sign}(\sigma\_{n+1}) \Big]. \end{split} \tag{A15}
$$

Recalling the trial state for the FKZ model:

$$\begin{split} \boldsymbol{\sigma}\_{n+1}^{\text{trial}} &= (1 + \mathsf{C}\_{4}^{\text{KZ}})^{-1} \Big[ \mathsf{C}\_{1}^{\text{KZ}} \Big( \Delta \boldsymbol{\varepsilon}\_{n+1} + \mathsf{H}^{\text{f}\_{2}} (\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^{\text{vp}}) \Big) + \mathsf{C}\_{2}^{\text{KZ}} \Big( \Delta \boldsymbol{\varepsilon}\_{n+1} + \mathsf{H}^{\text{f}\_{3}} (\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^{\text{vp}}) \Big) \\ &+ \mathsf{C}\_{3}^{\text{KZ}} \Big( \Delta \boldsymbol{\varepsilon}\_{n+1} + \mathsf{H}^{\text{f}\_{2} + \beta\_{3} - \beta\_{1}} (\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^{\text{vp}}) \Big) + \mathsf{C}\_{4}^{\text{KZ}} \Big( \boldsymbol{\sigma}\_{n} - \mathsf{H}^{\text{f}\_{2} - \beta\_{1}} \boldsymbol{\sigma} \Big) \Big] \end{split} \tag{A16}$$

which, by combining the above two equations, we find:

$$
\sigma\_{n+1} = \sigma\_{n+1}^{trial} - \text{sign}(\sigma\_{n+1}) \left( \frac{\mathbb{C}\_1^{KZ} + \mathbb{C}\_2^{KZ} + \mathbb{C}\_3^{KZ}}{1 + \mathbb{C}\_4^{KZ}} \right) \Delta \gamma\_{n+1}. \tag{A17}
$$

Finally, we obtain the loading/unloading sign consistency by following standard plasticity procedures:

$$\log(\sigma\_{n+1})|\sigma\_{n+1}| = \text{sign}(\sigma\_{n+1}^{trial})|\sigma\_{n+1}^{trial}| - \text{sign}(\sigma\_{n+1})\left(\frac{\mathbb{C}\_1^{KZ} + \mathbb{C}\_2^{KZ} + \mathbb{C}\_3^{KZ}}{1 + \mathbb{C}\_4^{KZ}}\right) \Delta \gamma\_{n+1}. \tag{A18}$$

therefore,

$$\text{sign}(\sigma\_{n+1})\left[|\sigma\_{n+1}|+\left(\frac{\mathbf{C}\_1^{KZ}+\mathbf{C}\_2^{KZ}+\mathbf{C}\_3^{KZ}}{1+\mathbf{C}\_4^{KZ}}\right)\Delta\gamma\_{n+1}\right] = \text{sign}(\sigma\_{n+1}^{trial})|\sigma\_{n+1}^{trial}|,\tag{A19}$$

since both terms multiplying the sign functions on the left and right sides of the above equation are positive, we therefore conclude that sign(*σn*+1) = sign(*σtrial <sup>n</sup>*+<sup>1</sup> ), and hence (A17) becomes:

$$
\sigma\_{n+1} = \sigma\_{n+1}^{trial} - \text{sign}(\sigma\_{n+1}^{trial}) \left( \frac{\mathbb{C}\_1^{KZ} + \mathbb{C}\_2^{KZ} + \mathbb{C}\_3^{KZ}}{1 + \mathbb{C}\_4^{KZ}} \right) \Delta \gamma\_{n+1},
\tag{A20}
$$

which completes the derivation.

#### **References**

