**1. Introduction**

In search of a generalization of existing models describing experimentally observed phenomena, the concept of fractional calculus [1,2] emerged as a tool that in the recent years became widely applied. Among the areas in which this theory has found application, it is worth mentioning mechanics where one can distinguish: (i) time-fractional models; (ii) space-fractional models; and (iii) stress-fractional models. For example, in [3], the time-fractional model was used to describe the time-dependent mechanical property evolution in ductile metals. The fractional oscillators were analyzed in [4], whereas the heat and mass transfer analysis in the framework of fractional calculus was presented in [5,6]. Furthermore, the analysis and modeling of turbulent flow in a porous medium [7], fluid transport induced by the osmotic pressure of glucose and albumin [8], wave propagation in the viscoelastic material [9], non-local boundary value problem [10], and evolution for the damage variable for hyperelastic materials [11] with an application of time-fractional derivative suggests grea<sup>t</sup> versatility of this approach. On the other hand, the space-fractional models are successfully used in mechanics to describe the deformation of a harmonic oscillator [12], deformation of an infinite bar subjected to a self-equilibrated load distribution [13], modeling plane strain and plane stress elasticity [14], Euler–Bernoulli beam [15], Darcy's flow in porous media [16] and fractional strain formulation [17]. Finally, the stress-fractional models [18,19] and their finite element implementations were used to study the granular soils under drained cyclic loading [20], and monotonic triaxial compression [21]. Concluding, one should emphasize that regardless of the specific formulation, the fractional operators have one common feature, namely, the 'change' of a selected variable is based on integration over a closed interval, thus extending the definition of integer order derivative (defined in a point) and simultaneously introducing a non-locality in a given space.

It is commonly accepted that the Theory of Thermo-Viscoplasticity (TTV), which plays a central role in the following considerations, began with the publication of Perzyna [22], which, until the present day, serves a basis for many efforts in linking experimental and numerical results for different types

of materials. The main results of this theory were discussed in a grea<sup>t</sup> number of papers that focused on phenomena such as propagation of mechanical and thermal waves [23,24], viscosity controlled by material parameter [25,26], dispersion [26], or implicit non-locality in the time variable [27]. Nonetheless, the classical TTV formulation does not include directional viscosity, and to include the non-normality extension needs, as all classical plasticity theories, postulation of an additional potential, which is not straightforward and causes the increase of material parameters. Furthermore, the same concern is relevant to the plastic anisotropy effects in terms of the original Perzyna model; to include this effect additional variables and evolution equations for them are needed to be postulated. This limitations were resolved by the generalization of the Perzyna formulation by definition of the fractional flow rule, first proposed in [18] and later developed in [19,28–30].

The implementation of the fractional plastic (rate independent) rule, for the Huber-Mises-Hencky (HMH) yield criterion, in the framework of implicit and explicit procedures and with examples on material point level, was presented in [28]. This was further developed in the subsequent article [29] to any smooth and convex yield criterion but still focusing on rate independent plastic flow. Concluding, in both these articles the non-locality in the stress state was present, however the implicit time non-locality common for the viscoplastic flow was not included in them.

This paper extends the concept of FV, which was first reported in [18], for the Initial Boundary Value Problem (IBVP), and provides a detailed discussion of the model material parameters. The parametric study includes the influence of the overstress power and the relaxation time (which is understood as implicit length scale parameter, as mentioned in [27]) on the dynamic properties of the FV model. Moreover, additional fractional material parameters, which induce the directional viscosity, the non-associative, and the anisotropic plastic flow, were also discussed.

### **2. Fractional Viscoplasticity**

### *2.1. Remarks on Fractional Calculus*

Fractional calculus (FC) introduces a new, universal method for calculating the intensity of changes of various quantities in mathematical models describing experimentally observed phenomena. FC implies a generalization of integer order derivatives, by fractional derivatives (FD). The selection of the FD definition (from an infinite number) can use a type of material as a criterion to obtain the best fitting of the constitutive model to a given experimental evidence. All definitions of the FD have a common property, namely they include summation over an interval abandoning the integer order derivative definition given at a single point; therefore they are called non-local. The classical derivative can be regarded a special case of the FD when its order becomes integer.

In order to explain the FD concept, let us consider a generalized fractional differential operator *BαP* as a composition of fractional integral *<sup>K</sup><sup>α</sup>P* with classical integer (*n*-th) differential operator [31]

$$B\_P^a = K\_P^{\mu - \alpha} \circ \frac{d^n}{dt^n} \tag{1}$$

where *α* is the order of the derivative, *n* = *α* + 1, · denotes the floor function, *P* is a parameter set (described below) and ◦ denotes the composition operator. *BαP* is referred to as the fractional differential operator *B* (*B*-op) of order *α* and *p*-set *P*, and analogously *<sup>K</sup><sup>α</sup>P* identifies the *K* (*K*-op) fractional integer operator of order *α* and *p*-set *P*.

The definition of K for the parameter set *P* = *a*, *t*, *b*, *p*, *q* can be given as

$$\left(\left(\mathbf{K}\_{\mathrm{P}}^{a}f\right)(\mathbf{t}) = p \int\_{a}^{t} k\_{a}(\mathbf{t}, \boldsymbol{\tau}) f(\boldsymbol{\tau}) d\boldsymbol{\tau} + q \int\_{t}^{b} k\_{a}(\boldsymbol{\tau}, \mathbf{t}) f(\boldsymbol{\tau}) d\boldsymbol{\tau},\tag{2}$$

where *t* ∈ *a*, *b* and *a* < *t* < *b*, *p*, *q* are real numbers, and *kα*(*<sup>t</sup>*, *τ*) is a kernel that depends on the order of the derivative *α*. It can be shown that if *kα* is a difference kernel, i.e., *kα*(*<sup>t</sup>*, *τ*) = *kα*(*<sup>t</sup>* − *τ*) and *kα* ∈ *L*1 ([0, *b* − *a*]) then *L*1 ([*b*, *a*]) → *L*1 ([*b*, *a*]) is well defined, bounded and linear. For explicit definition, the special form of the kernel function can be assumed

$$k\_a(t - \tau) = \frac{1}{\Gamma(a)} \left( t - \tau \right)^{a-1},\tag{3}$$

then for *P* = *a*, *t*, *b*, 1, 0 

$$(\left(K\_P^a f\right)(t) = \frac{1}{\Gamma(a)} \int\_a^t (t - \tau)^{a-1} f(\tau) d\tau = (\left. \left. I\_t^a f \right)(t) \right| \tag{4}$$

is obtained or, if *P* = *a*, *t*, *b*, 0, 1 then

$$(\left(K\_{\rm P}^{a}f\right)(t) = \frac{1}{\Gamma(a)} \int\_{t}^{b} (\tau - t)^{a-1} f(\tau) d\tau = (\iota I\_{b}^{a}f)(t),\tag{5}$$

where Γ is the Euler gamma function. Equations (4) and (5) describe the left and right Riemann-Liouville fractional integrals of the order *α*, respectively. The application of these operators in Equation (1) leads to the following fractional derivative definitions:

$$\left(\left(B\_P^a\right)f\left(t\right) = \,\_a^C D\_t^a f\left(t\right) = \frac{1}{\Gamma(n-a)} \int\_a^t \frac{f^{(n)}\left(\tau\right)}{(t-\tau)^{a-n+1}} d\tau,\tag{6}$$

for *t* > *a*, and

$$-\left(B\_{\mathcal{P}}^{a}\right)f\left(t\right) = \,^{\mathcal{C}}\_{t}D\_{b}^{a}f\left(t\right) = \frac{(-1)^{n}}{\Gamma(n-a)}\int\_{t}^{b}\frac{f^{(n)}\left(\tau\right)}{(\tau-t)^{a-n+1}}d\tau,\tag{7}$$

for *t* < *b*. The FD operators *Ca <sup>D</sup><sup>α</sup>t f* (*t*) and *Ct <sup>D</sup><sup>α</sup>b f* (*t*) are known as the left- and right-sided Caupto fractional integrals.

Finally, for the purpose of further definition of the FV, the Riesz-Caputo (RC) derivative can be expressed as a linear combination of previously given left and right Caputo derivatives

$$\, \, \_a^{RC}D\_b^a f\left(t\right) = \frac{1}{2} \left( \, \_a^C D\_t^a f\left(t\right) + (-1)^n \, \_t^C D\_b^a f\left(t\right) \right). \tag{8}$$

It can be shown that for the RC derivative the fundamental property of integer order derivatives is preserved, that is, the derivative of a constant is zero.

### *2.2. Basic Concepts*

In the following section Voigt notation is applied, thus the second rank tensors are ordered as (6 × 1 column matrix)

$$\mathbf{t} = \begin{pmatrix} t\_{11} \ t\_{22} \ t\_{33} \ t\_{23} \ t\_{13} \ t\_{12} \end{pmatrix}^{\mathrm{T}} = \begin{pmatrix} t\_1 \ t\_2 \ t\_3 \ t\_4 \ t\_5 \ t\_6 \end{pmatrix}^{\mathrm{T}},\tag{9}$$

whereas the fourth order tensors are represented by 6 × 6 matrices ordered in accordance with the rule used in Equation (9).

Deformation assumes the additive decomposition of total strain, therefore

$$
\mathfrak{e} = \mathfrak{e}^{\mathfrak{e}} + \mathfrak{e}^{vp},
\tag{10}
$$

or in a rate form

$$
\dot{\mathfrak{e}} = \dot{\mathfrak{e}}^{\mathfrak{e}} + \dot{\mathfrak{e}}^{vp},
\tag{11}
$$

where *ε* is the total strain, *εe* is the elastic strain and *εvp* is the viscoplastic strain. Next, due to thermodynamic restrictions, the elastic strain is related to elastic stress through Hooke's law

$$
\sigma^{\varepsilon} = \mathcal{L}^{\varepsilon} \mathfrak{e}^{\varepsilon},
\tag{12}
$$

where *σ<sup>e</sup>* denotes the Cauchy stress tensor and L*e* denotes the elastic constitutive tensor. The rate of viscoplastic strain is analogous to the classical viscoplastic definition, namely

$$
\dot{\varepsilon}^{vp} = \Lambda \mathbf{p}\_{\prime} \tag{13}
$$

where Λ is a scalar multiplier and **p** is the second order unit tensor which governs the direction of viscoplastic flow. As the **p** tensor is normalized, the magnitude of *ε*˙ *vp* depends solely on the Λ parameter.

Following the concept introduced by Perzyna [22], this parameter is expressed as

$$
\Lambda = \gamma \left< \Phi(F) \right> , \tag{14}
$$

where *γ* = 1*Tm* is the viscosity parameter, Φ is the overstress function that depends on the rate-independent yield surface *F*, and · is Macaulay brackets. It is well-known that *γ* introduces *implicit time non-locality* in the viscoplastic model [27]. Furthermore, the function Φ has the following form

$$\Phi(F) = F^{m\_{\text{op}}} = \left(\frac{\sqrt{f\_2}}{\kappa} - 1\right)^{m\_{\text{rep}}},\tag{15}$$

where √*J*2 denotes the second invariant of stress deviator and *κ* is the static yield stress in simple shear.

Finally, the remaining object needed to be defined is the tensor **p**. In this place, the difference between the classical theory of viscoplasticity and the new approach is most evident. Let us recall, that in the classical formulation the direction of yield is normal to yield surface and **p** can be written as

$$\mathbf{p} = \frac{\partial F}{\partial \sigma} \left( \left| \left| \frac{\partial F}{\partial \sigma} \right| \right| \right)^{-1}. \tag{16}$$

It is also well known, that Equations (16) and (15) indicate that the viscoplastic strain is coaxial with the deviatoric stress tensor (associated flow). As a result, the volume change can occur in the range of elastic deformations only. For modern materials such as metal-matrix composites, this assumption is no longer valid. Therefore, the constitutive model should be modified to capture this phenomenon.

The fractional approach assumes the application of the RC operator to **p** definition [18]. In such a case, Equation (16) is generalized to the form

$$\mathbf{p} = D^{\mathbf{a}} F \left| \left| D^{\mathbf{a}} F \right| \right|^{-1} \,, \tag{17}$$

where *Dα* stands for the RC operator (see Equation (8)). It is worth noting that the proposed formulation of **p** introduces the anisotropy of viscoplastic flow and furthermore (due to non-associativity) develops a tool to control the volume change in the plastic range of material behaviour [18]. Another essential remark is that Equation (17) introduces *explicit stress-fractional non-locality* in the overall model. It is important that the thermodynamic restrictions are formulated in a standard manner, and because of complicated structure of Equation (17) they are checked incrementally in the numerical procedure (see [29] for a detailed discussion).
