**3. Implementation**

Introduction of the three-dimensional fractional viscoplastic model requires a numerical procedure that governs the solution. Since our considerations are focused on extreme dynamic processes, the explicit time integration was chosen for finite element method. Therefore, the ABAQUS/Explicit code was utilized together with the user subroutine VUMAT. The critical steps of the implementation are presented below.

In the first step, Hooke's law is written as

$$\begin{Bmatrix} \sigma\_{11} \\ \sigma\_{22} \\ \sigma\_{33} \\ \sigma\_{23} \\ \sigma\_{13} \\ \sigma\_{12} \\ \sigma\_{12} \end{Bmatrix} = \begin{pmatrix} 2G+\lambda & \lambda & \lambda & 0 & 0 & 0 \\ \lambda & 2G+\lambda & \lambda & 0 & 0 & 0 \\ & \lambda & 2G+\lambda & 0 & 0 & 0 \\ & 0 & 0 & 0 & G & 0 & 0 \\ & 0 & 0 & 0 & 0 & G & 0 \\ & 0 & 0 & 0 & 0 & 0 & G \end{pmatrix} \begin{Bmatrix} \varepsilon\_{11} \\ \varepsilon\_{22} \\ \varepsilon\_{33} \\ \varepsilon\_{34} \\ 2\varepsilon\_{23} \\ 2\varepsilon\_{13} \\ 2\varepsilon\_{12} \end{Bmatrix},\tag{18}$$

where *G* = *E*/2(1 + *ν*) and *λ* = *<sup>E</sup>ν*/(<sup>1</sup> + *ν*)(<sup>1</sup> − <sup>2</sup>*ν*) are elastic constants and *E* and *ν* denote Young's modulus and Poisson's ratio, respectively. Next, because of application of the HMH yield criterion, the yield function in a matrix form is as follows

$$F(\sigma) = \begin{Bmatrix} \sigma\_{11} \\ \sigma\_{22} \\ \sigma\_{33} \\ \sigma\_{23} \\ \sigma\_{13} \\ \sigma\_{12} \\ \sigma\_{12} \end{Bmatrix}^{T} \begin{pmatrix} 1 & -\frac{1}{2} & -\frac{1}{2} & 0 & 0 & 0 \\ -\frac{1}{2} & 1 & -\frac{1}{2} & 0 & 0 & 0 \\ -\frac{1}{2} & -\frac{1}{2} & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 3 & 0 & 0 \\ 0 & 0 & 0 & 0 & 3 & 0 \\ 0 & 0 & 0 & 0 & 0 & 3 \end{pmatrix} \begin{Bmatrix} \sigma\_{11} \\ \sigma\_{22} \\ \sigma\_{33} \\ \sigma\_{33} \\ \sigma\_{13} \\ \sigma\_{12} \\ \sigma\_{13} \end{Bmatrix} - 3\kappa^{2} = 0. \tag{19}$$

Afterwards, the increment of viscoplastic strain (see Equations (13) and (17)) can be written in the form 

$$
\Delta \varepsilon^{vp} = \Delta t \Lambda \mathbf{p} = \Delta t \Lambda \frac{\left\{ \begin{array}{ll} D^a F & D^a F & D^a F & D^a F & D^a F & D^a F \end{array} \right\}^T}{||D^a F||} \tag{20}
$$

where (cf. [28])

$$\begin{pmatrix} p\_{11} \\ p\_{22} \\ p\_{33} \\ p\_{32} \\ p\_{13} \\ p\_{13} \\ p\_{12} \\ p\_{12} \end{pmatrix} = \begin{pmatrix} k\_{(11)}^{M} & -\frac{1}{2} k\_{(11)}^{M} & -\frac{1}{2} k\_{(11)}^{M} & 0 & 0 & 0 \\ -\frac{1}{2} k\_{(22)}^{M} & k\_{(22)}^{M} & -\frac{1}{2} k\_{(22)}^{M} & 0 & 0 & 0 \\ -\frac{1}{2} k\_{(33)}^{M} & -\frac{1}{2} k\_{(33)}^{M} & k\_{(33)}^{M} & 0 & 0 & 0 \\ & 0 & 0 & 0 & 3 k\_{(23)}^{M} & 0 & 0 \\ & 0 & 0 & 0 & 0 & 3 k\_{(13)}^{M} & 0 \\ & 0 & 0 & 0 & 0 & 0 & 3 k\_{(12)}^{M} \end{pmatrix} \begin{cases} \sigma\_{11} \\ \sigma\_{22} \\ \sigma\_{33} \\ \sigma\_{33} \\ \sigma\_{32} \\ \sigma\_{13} \\ \sigma\_{12} \end{cases} + \begin{cases} k\_{(11)}^{Q} \\ k\_{(21)}^{Q} \\ k\_{(21)}^{Q} \\ 3 k\_{(33)}^{Q} \\ 3 k\_{(23)}^{Q} \\ 3 k\_{(13)}^{Q} \\ 3 k\_{(12)}^{Q} \end{cases} \tag{21}$$

Symbols in Equation (21) denotes

$$k\_{\left(\dot{\boldsymbol{\eta}}\right)}^{\mathcal{M}} = \frac{\Gamma(2)}{\Gamma(2-\alpha)} \left[ \left(\Delta\_{\left(\dot{\boldsymbol{\eta}}\right)}^{L}\right)^{1-\alpha} + \left(\Delta\_{\left(\dot{\boldsymbol{\eta}}\right)}^{R}\right)^{1-\alpha} \right],\tag{22}$$

$$k\_{\left(\dot{\boldsymbol{\eta}}\right)}^{Q} = \left(\frac{\Gamma(2)}{\Gamma(2-a)} - \frac{1}{2} \frac{\Gamma(3)}{\Gamma(3-a)}\right) \left[ \left(\Delta\_{\left(\dot{\boldsymbol{\eta}}\right)}^{R}\right)^{2-a} - \left(\Delta\_{\left(\dot{\boldsymbol{\eta}}\right)}^{L}\right)^{2-a} \right],\tag{23}$$

$$
\Delta^{\mathbf{L}} = \left( \begin{array}{cccc} \Delta\_{(11)}^{\mathbf{L}} & \Delta\_{(22)}^{\mathbf{L}} & \Delta\_{(33)}^{\mathbf{L}} & \Delta\_{(23)}^{\mathbf{L}} & \Delta\_{(13)}^{\mathbf{L}} & \Delta\_{(12)}^{\mathbf{L}} \end{array} \right)^{T}, \tag{24}
$$

and

$$
\boldsymbol{\Delta}^{\mathbf{R}} = \left( \begin{array}{cccc} \boldsymbol{\Delta}^{\boldsymbol{R}}\_{(11)} & \boldsymbol{\Delta}^{\boldsymbol{R}}\_{(22)} & \boldsymbol{\Delta}^{\boldsymbol{R}}\_{(33)} & \boldsymbol{\Delta}^{\boldsymbol{R}}\_{(23)} & \boldsymbol{\Delta}^{\boldsymbol{R}}\_{(13)} & \boldsymbol{\Delta}^{\boldsymbol{R}}\_{(12)} \end{array} \right)^{\boldsymbol{T}},\tag{25}
$$

where

$$a\_{(ij)} = \sigma\_{\bar{i}\bar{j}} - \Delta\_{(\bar{i}\bar{j})'}^{L} \quad b\_{(\bar{i}\bar{j})} = \sigma\_{\bar{i}\bar{j}} + \Delta\_{(\bar{i}\bar{j})}^{\mathbb{K}}.\tag{26}$$

Terminals *<sup>a</sup>*(*ij*), *<sup>b</sup>*(*ij*) are needed to define the partial fractional derivatives in Equation (17) that enforce the directional nature of the fractional viscoplastic flow–the subscripts *L* and *R* corresponds to the left and the right Caputo derivatives, respectively. In addition, by introducing sections that extend the calculation beyond the material point, a virtual neighbourhood is obtained that results in a non-locality in a stress state.

The interpretation of the virtual surrounding in a stress state depends on the specific material (see [28]), but in general it could be understood as a (homogenized) phenomenological measure of some instability, e.g., for metallic materials it is connected with dislocation nucleation [32–34], nucleation of voids [35] or breakup of grains [36–38] (see review paper [39]). By way of illustration, Figure 1 shows the cross-section of this virtual neighbourhood in the *σ*2 − *σ*3 plane.

**Figure 1.** Virtual surrounding of a material point.

The analysis of Equation (21) shows differences in relation to the classical viscoplasticity where the change in volume may only occur in the elastic range—in the classical case, the trace of the **p** tensor equals 0. This condition is abandoned in the fractional formulation (when *α* ∈ (0, 1)), thus explicitly providing a tool to control the evolution of volume in the plastic range through *α* and parametric vectors **Δ<sup>L</sup>** and **ΔR**. Moreover, the versatility of the fractional approach is proven for *α* = 1, for which the associated plastic flow as a special case is obtained.

Finally, the flowchart was formulated that presents the general calculation scheme for the elasto-viscoplastic material in the framework of the fractional viscoplastic flow rule for explicit time integration in VUAMT subroutine (see Figure 2). The VUAMT subroutine aims at determination of the values of Cauchy stresses and updating strains and internal variables at time *tn*+<sup>1</sup> based on the knowledge of these parameters at the previous moment *tn*. The procedure starts with the calculation of the elastic trial stress, which is later used to establish the value of the yield criterion. If this criterion is fulfilled, the plastic multiplier Λ and the direction of plastic flow **p** are calculated according to the flow rule; otherwise the elastic step is conducted.

### **Given parameters:** *E*, *ν*, *κ*, *Tm*, *m*

**Figure 2.** VUMAT subroutine flowchart for the fractional viscoplastic rule.

### **4. Parametric Study: Uniaxial Tension**

### *4.1. Description of the Numerical Experiment*

The conducted parametric study is focused on the material point level represented by a unit cube with dimensions of 1 × 1 × 1 mm discretized by a single finite element C38DR (linear, eight-node brick with reduce integration). The boundary conditions required to achieve uniaxial constraints are shown in Figure 3. Basic mechanical properties were assumed as for the carbon steel, therefore the elastic range was characterized by Young's modulus *E* = 205 GPa and Poisson's ratio *ν* = 0.27. The fractional flow rule presented in Section 2 was applied in the plastic range. The static yield stress in simple shear for the selected material was *κ* = 605 MPa. Other material parameters, depending on the studied case, were chosen as described below.

**Figure 3.** Unit cubic model restricted to uniaxial tension.

The analyzed cases of fractional flow were divided into two groups to show how various combinations of model parameters influence plastic deformation. The first group is focused on the value of the stress-fractional non-locality spread (Δ*L*,*<sup>R</sup>*) and the order (*α*) of the fractional flow. In the second group, the influence of the material parameters *Tm* and *m* under various speeds of the imposed displacement is closely studied. Anticipating the anisotropic behaviour of the fractional material model, these two groups were further subdivided according to the direction where the dominant viscoplastic flow was expected. Hence, two cases were formed for tension direction (Δ*L*,*<sup>R</sup>* 22 = 0.005*κ* ≈ 3.0 MPa) and direction perpendicular to tension (Δ*L*,*<sup>R</sup>* 11 = 0.005*κ* ≈ 3.0 MPa). In each of those cases, other values of the Δ were set to 0.0017*κ* ≈ 1.0 MPa.

Two kinds of plots were used to present the results of the parametric study. The first kind exemplify the relation between three normal strains *ε*11,*ε*22,*ε*33. The second type shows the stress–strain relation in the tension (2) direction (it should be pointed out that for this kind of plots, a 'softening' is observed, especially for highest tension velocities, however, this effect is due to the lateral stresses induced by the inertia effects and is not due to constitutive model *κ* = *const*. See Figure 4, where this effect is negligible due to relatively small tension velocity *v* = 1 *ms* ). The research on influence of the fractional derivative order was performed for a set *α* ∈ {0.1, 0.25, 0.5, 0.75, 0.99, 1.0}—as mentioned earlier, fractional generalization of the viscoplasticity reduces to the classical solution for *α* = 1. The study of *Tm* and *m* was conducted for three different velocities of tension, i.e., *v* = 1, 25 and 50 *ms* .

**Figure 4.** Influence of the order *α* and the value of material parameter Δ22 on the stress–strain relation, for: *v* = 1 *ms*, *Tm* = 2.5*e*-6 *s*, *m* = 1.

### *4.2. Influence of the Order of FV and Non-Locality in a Stress State on Plastic Flow*

4.2.1. Study of Intensified Plastic Flow in Tension Direction for Different Orders of Flow

Figure 4 presents the material response to the applied tension velocity of *v* = 1 *ms* for different flow intensities in tension direction and flow orders. Increasing the flow intensity parameter (Δ*L*,*<sup>R</sup>*) causes higher evolution of the plastic flow in the chosen direction but in this case the velocity of the load is not sufficient to reveal different behaviour in the stress–strain relation for different values of *α* (see Figure 4). Next, for the same configuration of material parameters higher velocity is applied, namely *v* = 25 *ms* At this speed (Figure 5), a slight waveform begins to be visible for Δ*L*,*<sup>R</sup>* 22 = 1.0 MPa. The amplitude of the stress signal increases with the increases Δ*L*,*<sup>R</sup>* 22 . Additionally, the influence of *α* is shown because the decrease in its value translates into greater amplitude of oscillations. So, we conclude that both fractional parameters, control the dynamic properties of a fractional model.

.

**Figure 5.** Influence of the order *α* and the value of material parameter Δ22 on the stress–strain relation, for: *v* = 25 *ms*, *Tm* = 2.5*e*-6 *s*, *m* = 1.

4.2.2. Study of Intensified Plastic Flow Perpendicular to the Tension Direction for Different Orders of Flow

Results presented in this section were obtained for a parameter set similar to this in Section 4.2.1 with the difference that flow intensity is increased in the direction perpendicular to tension, namely Δ*L*,*<sup>R</sup>* 11 = 3.0 MPa. Others components of vectors in Equations (24) and (25) equal 1. As in the discussion in the previous section, the tension velocity of *v* = 1 *m s* is insufficient to reveal the influence of *α* on the stress–strain relation (see Figure 6). However, Figure 7 shows that the material prefers to deform in (1) direction when the magnitude of Δ*L*,*<sup>R</sup>* 11 growths, hence the *<sup>ε</sup>*11/*<sup>ε</sup>*33 ratio is greater then 1. Moreover, the intensity of the flow in the preferred direction increases as the value of *α* diminishes to 0. Next, as before, the velocity is increased to *v* = 25 *m s* . Figure 8 shows that when the value of Δ*L*,*<sup>R</sup>* rises, greater amplitude of the oscillation and hardening of the material can be observed. This last effect is inversely proportional to the order of the fractional flow. The relation between *ε*11,*ε*22 and *ε*33 is presented in Figure 9 and is very similar to Figure 7 with the only distinction that slight oscillation occurs as a result of higher velocity. So, we conclude that both fractional parameters control the anisotropic properties of a fractional model in the plastic range.

**Figure 6.** Influence of the order *α* and the value of material parameter Δ11 on the stress–strain relation, for: *v* = 1 *m s*, *Tm* = 2.5*e*-6 *s*, *m* = 1.

**Figure 7.** Influence of the order *α* and the value of material parameter Δ11 on the relation between three normal stresses, for: *v* = 1 *m s*, *Tm* = 2.5*e*-6 *s*, *m* = 1.

**Figure 8.** Influence of the order *α* and the value of material parameter Δ11 on the stress–strain relation, for: *v* = 25 *ms*, *Tm* = 2.5*e*-6 *s*, *m* = 1.

**Figure 9.** Influence of the order *α* and the value of material parameter Δ11 on the relation between three normal stresses, for: *v* = 25 *ms*, *Tm* = 2.5*e*-6 *s*, *m* = 1.

### *4.3. Influence of the Relaxation Time and the Overstress Power*

4.3.1. Study of the Fractional Flow Under Different Dynamic Loading Rates for Intensified Plastic Flow in Tension Direction

Here we assume that the intensified plastic flow, determined by Δ*L*,*<sup>R</sup>* 22 = 3.0 MPa, is in tension direction. Figure 10 presents the effect of different relaxation times for various velocities of tension. It should be pointed out that in order to increase clarity of interpretation both the classical (*α* = 1) and the fractional (*α* = 0.75) solutions are compared on each graph. As can be seen, when the relaxation time grows, the hardening of the material as well as the stress waves oscillations increase. The latter is especially pronounced for the relaxation time *Tm* = 2.5*e*-5 *s*.

**Figure 10.** Influence of the relaxation parameter *Tm* and the value of applied velocity field *v* on the stress–strain relation, for: *α* = 0.75, *m* = 1, Δ22 = 3.0.

Figure 11 presents the result of increasing the value of the overstress parameter *m*. For the fractional (*α* = 0.75) viscoplastic material the stress level is smaller in relation to the classical (*α* = 1.0) viscoplastic solution (Figures 10 and 11).

**Figure 11.** Influence of the material parameter *m* and the value of applied velocity field *v* on the stress–strain relation, for: *α* = 0.75, *Tm* = 2.5*e*-6 *s*, Δ22 = 3.0.

Results discussed above indicate that the relaxation time and overstress power, together with fractional parameters, control the level of strain rate hardening and stress waves oscillation amplitude.

4.3.2. Study of the Fractional Flow Under Different Dynamic Loading for the Intensified Plastic Flow Perpendicular to the Tension Direction

In this section, it is assumed that fractional flow is intensified in the direction perpendicular to tension load (Δ*L*,*<sup>R</sup>* 11 = 3.0 *MPa*). Figure 12 shows that raising relaxation time causes similar effects to those discussed in the previous section. The biggest change occurs for the relaxation time *Tm* = 2.5*e*-5 *s* both in the stress level and the stress wave oscillations frequency. In Figure 13 a clear anisotropy of the plastic deformation for *α* = 0.75 can be noticed. As before, it is observed that increasing the value of *m* causes a growth in the material strain hardening without any apparent influence on the frequency of the stress wave (see Figure 14). By analogy to what is presented for *Tm*, in Figure 15 the dominant nature of *ε*11 can be observed, in addition to more pronounced oscillations for greater values of *m*. The analysis of the *σ* − *ε* relation for *α* = 1 and *α* = 0.75 also revealed that the stress levels of the fractional model are generally greater than for the classical solution (see Figures 12 and 14). This last observation is different from what was discovered in the previous section.

**Figure 12.** Influence of the relaxation parameter *Tm* and the value of applied velocity field *v* on the stress–strain relation, for: *α* = 0.75, *m* = 1, Δ11 = 3.0.

**Figure 13.** Influence of the relaxation parameter *Tm* and the value of applied velocity field *v* on the relation between three normal stresses, for: *α* = 0.75, *m* = 1, Δ11= 3.0.

In conclusion, as before, a prominent influence of the relaxation time and overstress power, together with fractional parameters, on the level of strain rate hardening and stress waves oscillation amplitude is observed. Additionally, the impact on dynamic properties of deformation anisotropy for fractional material was confirmed.

**Figure 14.** Influence of the material parameter *m* and the value of applied velocity field *v* on the stress–strain relation, for: *α* = 0.75, *Tm* = 2.5*e*-6 *s*, Δ11 = 3.0.

**Figure 15.** Influence of the material parameter *m* and the value of applied velocity field *v* on the relation between three normal stresses, for: *α* = 0.75, *Tm* = 2.5*e*-6 *s*, Δ11 = 3.0.

### *4.4. Study of the Disperse Character of the Fractional Viscoplastic Stress Waves*

The last study performed was the analysis of the disperse character of the fractional viscoplastic stress waves. As shown in the previous sections, the uniaxial dynamic deformation induces a stress waves. The regularity of the oscillations was determined by averaging intervals between peaks and then calculating the frequency. Table 1 lists the stress wave frequencies for Δ*L*,*<sup>R</sup>* 22 = 3.0 MPa and Δ*L*,*<sup>R</sup>* 11 = 3.0 MPa, which were depicted in the middle graphs of Figures 10 and 12. The material parameter Δ*L*,*<sup>R</sup>* does not appear in the classical viscoplasticity, so there is no change in the stress wave frequency for various values of this parameter when *α* = 1, thus tension velocity is only important. However, for fractional material, i.e., when *α* = 0.75, both flow intensity parameters and tension velocities modulate the frequency of stress waves. For *α* = 0.75 the frequencies are higher when the distinguished direction is co-linear with tension (Δ22 = 3.0 MPa).


.

**Table 1.** The stress wave frequencies for *Tm* = 2.5*e*-6 *s*, *m* = 1

The results presented in Table 2 correspond to the investigation of the role of the relaxation time and overstress power discussed in Section 4.3. Regardless of the value of *α*, the largest change in frequency can be observed between columns 2 and 3, that is for *Tm* = 2.5*e*-6 *s* and *Tm* = 2.5*e*-5 *s*. The comparison of the values in columns 1 (*Tm* = 2.5*e*-7 *s*) and 2 (*Tm* = 2.5*e*-6 *s*) shows that the stress wave frequencies are the same or very similar. The impact of *m* shows that the significant influence of this parameter was only recorded for *Tm* = 2.5*e*-5 *s*. For both *α* = 0.75 and *α* = 1 the increase of the overstress power results in the increase of the frequency of the stress wave.

**Table 2.** The stress wave frequencies for Δ22 = 3.0 *MPa*, *v* = 50 *ms*


One concludes that the relaxation time and overstress power, together with fractional parameters, control the dispersive character of stress waves, and even more, makes this attribute directional.

.
