Modified elastic-plastic model: implementation algorithm and comparison of computational efficiency with the elastic-viscoplastic model

Despite the intense integration of polymer and composite materials into modern manufacturing, metal and alloy products remain in high demand in different industries. The constantly evolving range of alloys and alloy-based products, the severe requirements for their operational characteristics, and the need to rapidly develop manufacturing technologies (mostly, by applying plastic deformation methods) are a strong motive for creating mathematical models (MM), which can adequately describe thermomechanical processing (TMP) of metallic materials. The "heart" of such MMs, determining their adequacy and accuracy, are the constitutive models (CM) (constitutive relations (CR)) used for their formulation.
Up to now, the processes of plastic working of metals and alloys are normally described by the mathematical models that are based on classical macrophenomenological theories of plasticity (MPT) [1 and others]. Unfortunately, the classical theories of plasticity (elastoviscoplasticity, creep, etc.) do not allow an explicit description of the evolving structure of materials, which determines the mechanical properties of the material.
This circumstance acted as an impetus to the creation of alternative approach, which was developed in the 20th century based on the introduction of internal variables [2 and others]. The latter are currently taken to mean parameters, characterizing the structure of the material at various structural-scale levels. This approach provides the framework for the development of multilevel CM based on the crystal plasticity (elastoplasticity, elastoviscoplasticity) [3].
The first two-level CM was proposed by J.I. Taylor [4], a rigorous mathematical justification of which was developed by J. Bishop and R. Hill [5]. In the literature, different versions of the models, which are based on the basic statements of the pioneering works, are generally called the models of Taylor-Bishop-Hill type (TBH models). The original TBH model is based on the assumption of rigid-plastic behavior of the material, which leads to the appearance of constraint-the material incompressibility. As is well known, for materials with constraints, it is impossible to determine the response (stress) solely in terms of its motion (strain) (K. Truesdell [6]). Under the imposed condition of material incompressibility only the deviator stress tensor can be determined in terms of the strain, whereas solving of boundary value problems (equations of motion) requires knowledge of all six components of the (symmetric) stress tensors.
In TBH-type models, the implementation of which can be treated as the process of solving a linear programming problem, the stress deviator is determined in a five-dimensional space of deviators at the vertices of the constraint polyhedron, which is the yield polyhedron established by Schmid's law. In this case, to determine all components of the stress deviator, it is necessary to use the coordinates of the vertex of the yield polyhedron, which are formed by the intersection of five linearly independent hyperplanes in a 5 -dimensional space. Note however that more than 5 hyperplanes can intersect at the vertices of the yield polyhedron; their number is called the vertex order; for example, in FCC crystals the vertices in the reference natural configuration can be of the 5th, 6th, and 8th orders. At the same time, the deviatoric space has only five linearly independent basis tensors, that is, any set of five orientation tensors can serve as a basis. However, in the most commonly encountered crystallites (with FCC, BCC and HCP lattices), there can be several such sets for the vertices of the yield polyhedron, which leads to ambiguity in the choice of active systems and the
determination of slips along them. Taylor suggested using the set that minimizes the dissipation intensity (the rate of work done by the shear stress at slip velocities). This eliminates the uncertainty problem, but has no clear physical justification. Any slip system, for which the Schmid condition is satisfied, "enjoys equal rights" with all others for which the latter is true at the moment of deformation. As shown in some works, depending on the selection of sets of active slip systems (SS), the solution of problems on the deformation of single- and polycrystalline samples can yield different results [7].
The later work by T.G. Lin [8] took into account elastic deformation, eliminating thereby the disadvantages associated with the presence of a constraint (incompressibility) and the possibility of implementing elastic-plastic deformation by activating less than five slip systems. However, the most important drawback - the uncertainty in the selection of a set of active slip systems still exists. It should be emphasized that the limitation of the number of slip systems to five in the case, when the representative point in the stress space falls on the vertex of higher than the fifth order, is determined only by the procedure of determining the rates (or increments) of shears and stresses. There is no physical substantiation of such limitation.
Two-level elastic-viscoplastic (i.e., strain rate-sensitive) models that appeared in the 1970s are devoid of this limitation [910 and others]. It has been shown that when the strain-rate sensitivity parameter tends to zero, the required solution converges to solving the problem based on the elastic-plastic model [11]. However, in this case the system of equations becomes stiff, which necessitates the application of implicit integration schemes and significantly reduces the computational efficiency [11]. In view of the above circumstance extensive studies have been attempted to remove the major drawback of TBH- type models.
The simplest methods for solving the problem of uncertainty are as follows: random selection of slip systems [12], usage of the average over all possible sets of slips over active slip systems [13], determination of active slip systems by introducing random disturbances into critical stresses [14], selection of 5 active slip systems based on the amount of their potential energy [15], iterative procedure for the selection of 5 SSs , providing the maximum contribution when expanding the velocity deformation tensor into SS orientation tensors [16], the use of generalized inverse (or pseudoinverse) matrix for solving the system of equations from which shear rates can be determined [11, 17].
Another approach to solving this problem involves modification of the hardening law. In [18], the authors consider the condition for uniqueness of the solution to the problem of shear strain determination in relation to the form of the hardening law used. It is shown that in the case when the matrix of hardening moduli $H_{ij}$$H_{ij}$H_(ij)H_{i j}$H_{ij}$ is positive definite, there exists a unique solution to the problem of determining shear rates. However, the cited work does not concern the selection of active SSs and determination of the shear rates on them. It is implicitly assumed that no more than 5 SSs can be concurrently active. Moreover, the laws of hardening depend on the properties of specific materials and cannot be "adapted" to mathematical procedures.
In some works, energy considerations are used to eliminate the above disadvantage of TBH-type models. In [19], the determination of the best set of shear increments for potentially active SSs at a loading step involves the fulfillment of an additional condition of "second-order minimality" of work done on plastic deformation (in authors' notation), i.e. $\mathrm{\Delta }(\mathrm{\Delta }\mathrm{A})=\mathrm{\Delta }(\sum _{k=1}^K\mathit{\sigma }:{\mathbf{M}}^{(k)}\mathrm{\Delta }{\gamma }^{(k)})\to min$$\mathrm{\Delta }(\mathrm{\Delta }\mathrm{A})=\mathrm{\Delta }\left(\sum _{k=1}^K \mathit{\sigma }:{\mathbf{M}}^{(k)}\mathrm{\Delta }{\gamma }^{(k)}\right)\to min$Delta(DeltaA)=Delta(sum_(k=1)^(K)sigma:M^((k))Deltagamma^((k)))rarr min\Delta(\Delta \mathrm{A})=\Delta\left(\sum_{k=1}^{K} \boldsymbol{\sigma}: \mathbf{M}^{(k)} \Delta \gamma^{(k)}\right) \rightarrow \min$\mathrm{\Delta }(\mathrm{\Delta }\mathrm{A})=\mathrm{\Delta }(\sum _{k=1}^K\mathit{\sigma }:{\mathbf{M}}^{(k)}\mathrm{\Delta }{\gamma }^{(k)})\to min$, where $\mathrm{\Delta }\mathrm{A}$$\mathrm{\Delta }\mathrm{A}$DeltaA\Delta \mathrm{A}$\mathrm{\Delta }\mathrm{A}$ is the increment of work at the loading step, $\mathrm{\Delta }{\gamma }^{(k)}$$\mathrm{\Delta }{\gamma }^{(k)}$Deltagamma^((k))\Delta \gamma^{(k)}$\mathrm{\Delta }{\gamma }^{(k)}$ is the increment of slip on the $k-$$k-$k-k-$k-$ th slip system, $\mathit{\sigma }$$\mathit{\sigma }$sigma\boldsymbol{\sigma}$\mathit{\sigma }$ is the Cauchy stress tensor, ${\mathbf{M}}^{(k)}$${\mathbf{M}}^{(k)}$M^((k))\mathbf{M}^{(k)}${\mathbf{M}}^{(k)}$ is the orientation tensor of the $k-$$k-$k-k-$k-$ th SS . The study fails to provide rigorous proof of this hypothesis, although it demonstrates good agreement with experimental data.
In [20], to eliminate the ambiguity in determining the set of SSs in the models of the TBH type, a numerical procedure for maximizing the Lagrangian (stress power - the rate of change of the strain energy density on plastic deformation) is proposed. Two additional terms are introduced into the Lagrangian; the first term, which is introduced via Lagrange multipliers, is responsible for the fulfillment of plasticity condition. The second term is introduced as a penalty function (the effective stress exceeds the value of the flow stress) with a penalty parameter represented by fictitious viscosity; this term is responsible for the fulfillment of condition of consistency of plastic deformation. At the same time, it should be noted that an unambiguous determination of shear rates is impossible due to the linear dependence of the relations, establishing limitations. They are calculated using a special algorithm, which the author calls "pseudo-inverse".
In [21], two algorithms are developed based on the incremental variational principles, the minimization of which yields both the equation for conservation of momentum and the constitutive relations. To provide continuous active loading, the procedure of "return to yield polyhedron" (return mapping algorithm) or the analogous stress design method is used. The problems of constraint degeneracy and solution ambiguity are eliminated by regularization and penalty functions, which implies that at each load increment it is necessary to solve a nonlinear optimization problem. However, the final solution strongly depends on the selected parameters (regularization strength, penalty multiplier), for which there are no physically justified criteria of selection.

In [22] it is suggested to use the principle of maximum work on shear rates to determine sets of active SSs at each loading step in the framework of elastoplastic model. The results of deformation calculations for a FCC polycrystal are compared with the data obtained within the framework of elastoviscoplastic model with a power-law viscoplasticity and a high exponent. It is noted that the results are close only for strain hardening, but in the case of significant latent hardening there is some difference.
Another technique borrowed from the macrophenomenological theories of plasticity is the replacement of the singular yield surface (polyhedron) with a smooth surface [23]. In this technique, the elastoplastic model is an analogue of the elastoviscoplastic CM.
Based on the above brief review of works dealing with non-uniqueness of the choice of a set of active SSs in such models as TBH (physical elastoplastic theories), we can ascertain that currently there is no generally accepted approach to resolution of this issue. Most of the studies attempt to justify the selection of the set consisting of five active SSs. In the articles, which allow for simultaneous activation of a larger number of SS, the determination of shear rates (or increments) is based on the complex mathematical procedures for which there is no proper physical justification
It seems that the main problem is the inconsistency between the dimension of space, in which velocities or shear increments are sought (per se, the five-dimensional space of Cauchy stress deviators) and the dimension of space of deviators of the plastic component of the velocity gradient (or displacement increment), the dimension of which is eight. It should be noted that the definition of Cauchy stress does not imply its symmetry. Due to the presence of an elastic law written in terms of asymmetric measures of stress and deformation, there is simply no space for the problem of nonuniqueness under consideration, yet thus far there are no such formulations of the analogue of Hooke's law.
In this work, we consider a possible variant of overcoming the indicated disadvantage of elastoplastic models of the TBH type by making use of the well-known provisions of crystal plasticity [24]. A two-level statistical constitutive elastoplastic model is used to describe the behavior of a representative macrovolume element (mono- or polycrystalline body) [3]. The constitutive model is developed based on Schmid's law, any of the possible hardening laws for slip systems, the hypoelastic law, and the decomposition of motion into quasi-solid and deformation. Due to the nonlinearity of CM, a stepwise procedure (time step or nondecreasing parameter increment) is applied. At the current time of loading, the activity of the SS is determined by the Schmid condition, and regardless of the number of activated systems, they all are considered "equal." To determine shear rates (increments) at each loading step, in addition to the hypoelastic relation and the law of hardening, it is proposed to use an iterative procedure that involves the decomposition of the plastic component of the strain rate measure into linearly independent dyads of slip system ("oblique basis" of the space of asymmetric deviators). The maximum number of basic dyads is eight. It should be noted that due to the latter, an excess of the number of simultaneously active slip systems (equal to the degree of vertices of the yield polyhedra, i.e. 8) is improbable regardless of the type of lattice. In the known works on physical theories of elastoplasticity, there are no examples of yield polyhedra with vertex degrees higher than eighth.
It should be noted that up to now, the object of study of most works on multi-level physical-oriented constitutive models are limited to a representative macrovolume. At the same time, for solving real problems of creation of new technologies for processing metals and alloys by plastic deformation (especially for producing functional materials and products), the issue of computational efficiency of mathematical models of technological processes, and thus of constitutive models, is very important. The papers cited in the review lay emphasis on the fact that the use of elastoviscoplastic models involves significantly higher computational costs compared to elastoplastic ones. These circumstances determine the authors' interest in the development of computationally efficient modification of the elastoplastic constitutive model devoid of the above-mentioned disadvantages of TBH-type models.

Processing of metals and alloys by methods of plastic deformation is implemented, as a rule, under conditions of large displacement gradients, so that for formulating and solving the corresponding boundary value problems, it is necessary to take into account the contact-type boundary conditions. In this regard, boundary value problems arising in the study of these processes belong to the class of geometrically and physically nonlinear, which requires the formulation of specific CM. One of the main difficulties in constructing geometrically nonlinear CMs is the decomposition of motion into a quasi-solid and deformational motion; for a continuum, this problem still does not have a single-valued solution. A possible variant of such decomposition for crystalline solids was proposed earlier in the work of [3]. Due to a priori unknown actual configuration of the examined region at each moment of the deformation process,
solutions to geometrically nonlinear contact problems should be developed in the velocity formulation using the stepwise (in time or non-decreasing parameter) procedures.
The constitutive model under consideration is the rate analogue of Lin's model [8]. It includes submodels of two structural-scale levels; at the macro level, the response of a representative macrovolume element is considered (a single crystal or a polycrystalline aggregate containing a sufficient number of crystallites for statistical averaging), and at the meso-level - the response of crystallites (grains, subgrains, fragments). To simplify the description of the model, the case of isothermal loading is considered. The behavior of the material at the upper level is considered in terms of mechanical variables - stress tensors, strains and their rates, satisfying the requirement of independence from the choice of reference system [6]. At the meso-level, the description of the deformation process is carried out in terms of discrete-continuous variables - shear rates, shear stresses on slip systems (SS), measures of distortion and rotation tensors of the crystallite lattice. The main mechanism of inelastic deformation is the movement of edge dislocations along the slip systems, the position of which is known for each type of crystal lattice. "Relational" variables are designated by the same letters, at the macro level - in capital letters, at the meso-level - in lower case letters. The link of submodels of the meso- and macrolevels is accomplished through the generalized Voigt hypothesis and the operations of averaging the stress tensors, elastic properties, spin, and the plastic part of the velocity deformation tensor.
The system of equations used to describe the behavior of the material of a representative volume at the macro level is written as

where $\mathbf{\Sigma }$$\mathbf{\Sigma }$Sigma\boldsymbol{\Sigma}$\mathbf{\Sigma }$ is the stress tensor of the macrolevel, $r$$r$rr$r$ is the designation of the co-rotational derivative, $\dot{X}$$\dot{X}$X^(˙)\dot{X}$\dot{X}$ is the designation of the material time derivative of the variable $\mathbf{X},\mathbf{\Pi }$$\mathbf{X},\mathbf{\Pi }$X,Pi\mathbf{X}, \boldsymbol{\Pi}$\mathbf{X},\mathbf{\Pi }$ is the tensor of elastic properties of the macrolevel, $\mathbf{Z},{\mathbf{Z}}^{\mathrm{e}},{\mathbf{Z}}_{\mathrm{p}}$$\mathbf{Z},{\mathbf{Z}}^{\mathrm{e}},{\mathbf{Z}}_{\mathrm{p}}$Z,Z^(e),Z_(p)\mathbf{Z}, \mathbf{Z}^{\mathrm{e}}, \mathbf{Z}_{\mathrm{p}}$\mathbf{Z},{\mathbf{Z}}^{\mathrm{e}},{\mathbf{Z}}_{\mathrm{p}}$ is the asymmetrical measure of the strain rate of the macrolevel, and its elastic and plastic components, respectively, $\mathbf{F}$$\mathbf{F}$F\mathbf{F}$\mathbf{F}$ is the deformation affinor (transposed deformation gradient), $\mathbf{L}$$\mathbf{L}$L\mathbf{L}$\mathbf{L}$ is the transposed gradient of the velocity of macrolevel movements, $\mathbf{\Omega }$$\mathbf{\Omega }$Omega\boldsymbol{\Omega}$\mathbf{\Omega }$ is the spin of the macrolevel, $<\cdot >$$<\cdot >$< * ><\cdot>$<\cdot >$ is the designation of the averaging over the representative volume of the macrolevel.
The system of equations describing the behavior of meso-level elements has the following form (designation of crystallite numbers is omitted):

where ${\dot{\gamma }}^{(k)}$${\dot{\gamma }}^{(k)}$gamma^(˙)^((k))\dot{\gamma}^{(k)}${\dot{\gamma }}^{(k)}$ is the shear rate along the $k$$k$kk$k$-th slip system; ${\tau }_c^{(k)},{\dot{\tau }}_c^{(k)}$${\tau }_c^{(k)},{\dot{\tau }}_c^{(k)}$tau_(c)^((k)),tau^(˙)_(c)^((k))\tau_{c}^{(k)}, \dot{\tau}_{c}^{(k)}${\tau }_c^{(k)},{\dot{\tau }}_c^{(k)}$ are the critical stress on the $k$$k$kk$k$-th slip system and the rate of its change determined by the adopted hardening model (see below); ${\mathbf{m}}^{(k)}={\mathbf{b}}^{(k)}{\mathbf{n}}^{(k)}$${\mathbf{m}}^{(k)}={\mathbf{b}}^{(k)}{\mathbf{n}}^{(k)}$m^((k))=b^((k))n^((k))\mathbf{m}^{(k)}=\mathbf{b}^{(k)} \mathbf{n}^{(k)}${\mathbf{m}}^{(k)}={\mathbf{b}}^{(k)}{\mathbf{n}}^{(k)}$ is the orientation tensor of the $k_{-}$$k_{-}$k_(-)k_{-}$k_{-}$ th slip system; ${\mathbf{b}}^{(k)},{\mathbf{n}}^{(k)}$${\mathbf{b}}^{(k)},{\mathbf{n}}^{(k)}$b^((k)),n^((k))\mathbf{b}^{(k)}, \mathbf{n}^{(k)}${\mathbf{b}}^{(k)},{\mathbf{n}}^{(k)}$ are the unit vectors of the sliding direction and the normals to the sliding plane; $K$$K$KK$K$ is the SS number (duplicated number of systems is used, so that in active systems shear rates and shear stresses are always positive). Note that Eqn. (2)6 is valid only for active SS (i.e. SS, for which the Schmid condition is satisfied). To decompose the motion into quasi-solid and deformation ones, a rigid Cartesian moving coordinate system (MCS) with a basis ${\mathbf{k}}_i$${\mathbf{k}}_i$k_(i)\mathbf{k}_{i}${\mathbf{k}}_i$ associated with the lattice is introduced in each crystallite [3]; in the numerical implementation of the algorithm at the mesolevel, the components of the stress tensors are determined in the orthonormal basis of the MCS. To determine the spin of the MSC, the lattice rotation model is used [3]. The initial conditions are set based on the hypothesis of a natural configuration.
Within the framework of this work, the quasilinear anisotropic law is adopted as the hardening law:

where $q_{\text{lat}},\alpha ,{\mathrm{h}}_0,{\tau }_s$$q_{\text{lat }},\alpha ,{\mathrm{h}}_0,{\tau }_s$q_("lat "),alpha,h_(0),tau_(s)q_{\text {lat }}, \alpha, \mathrm{h}_{0}, \tau_{s}$q_{\text{lat }},\alpha ,{\mathrm{h}}_0,{\tau }_s$ are the parameters of the hardening law, $\delta (k\mapsto )$$\delta (k\mapsto )$delta(k|->)\delta(k \mapsto)$\delta (k\mapsto )$ is the Kronecker delta, ${\tau }_s$${\tau }_s$tau_(s)\tau_{s}${\tau }_s$ is the saturation stress (for the FCC lattices are assumed to be the same for all SSs).
The reference configuration is assumed to be natural (unstressed), the displacements are assumed to be zero, and the critical stresses are determined by the properties of the lattice. The (transposed) gradient of the movement velocity $\mathbf{L}=\hat{\mathrm{\nabla }}{\mathbf{V}}^{\mathrm{T}}$$\mathbf{L}=\widehat{\mathrm{\nabla }}{\mathbf{V}}^{\mathrm{T}}$L= hat(grad)V^(T)\mathbf{L}=\hat{\nabla} \mathbf{V}^{\mathrm{T}}$\mathbf{L}=\hat{\mathrm{\nabla }}{\mathbf{V}}^{\mathrm{T}}$ ( $\hat{\mathrm{\nabla }}$$\widehat{\mathrm{\nabla }}$hat(grad)\hat{\nabla}$\hat{\mathrm{\nabla }}$ is the nabla operator defined in the current configuration, $\mathbf{V}$$\mathbf{V}$V\mathbf{V}$\mathbf{V}$ is the movement velocity vector), which is uniform for the examined representative macrovolume, is taken as the imposed action and specified as a continuous tensor-valued function (of the second rank) of time. Determination of the response of a representative macrovolume is carried out using a stepwise procedure. At each step, the solution is implemented in three stages: 1) solving the problem in terms of velocity, 2) integration, determination of the values of the sought variables at the end of the time step (all operations with tensors are performed in terms of the components defined in the MCS basis), 3) redefinition of the MCS basis orientations in crystallites, determination of tensor variables for all crystallites taking into account the rotation of the MCS axes, averaging of variables over the representative macrovolume. The algorithm for the implementation of the basic elastoplastic model is described in more detail in the article [7].

As noted above, the problem of ambiguity (non-uniqueness) in the selection of sets of active slip systems is due to the numerical procedure for determining shear rates in elastic-plastic models (such as the Lin model). In the given formulation of the mesoscale submodel, this circumstance is associated with the solution of Eqn. (2)6. Despite the fact that this equation uses asymmetric orientation dyads (deviators) ${\mathbf{m}}^{(k)}$${\mathbf{m}}^{(k)}$m^((k))\mathbf{m}^{(k)}${\mathbf{m}}^{(k)}$, their symmetrization is carried out implicitly (due to the symmetry of the tensor of elastic characteristics $\mathbf{\Pi }$$\mathbf{\Pi }$Pi\boldsymbol{\Pi}$\mathbf{\Pi }$ with respect to the indices of the first and second pairs ${\mathrm{\Pi }}_{ijkl}={\mathrm{\Pi }}_{jik}={\mathrm{\Pi }}_{ijlk}$${\mathrm{\Pi }}_{ijkl}={\mathrm{\Pi }}_{jik}={\mathrm{\Pi }}_{ijlk}$Pi_(ijkl)=Pi_(jik)=Pi_(ijlk)\Pi_{i j k l}=\Pi_{j i k}=\Pi_{i j l k}${\mathrm{\Pi }}_{ijkl}={\mathrm{\Pi }}_{jik}={\mathrm{\Pi }}_{ijlk}$ ). Therefore, the number of linearly independent equations in system (2) ${}_6$${}_6$_(6){ }_{6}${}_6$ does not exceed five (the dimensions of the space of symmetric deviators of the 2nd rank).
At the same time, from physical considerations, all SSs, for which the Schmid criterion is satisfied at the current moment of deformation, should be considered to have equal rights of being recognized as active, which is used in constructing an iterative procedure for solving the problem. The latter is integrated into the general algorithm for the implementation of a two-level elastic-plastic model destined for determining the shear rates in the mesoscale submodel at the first stage (solution in term of rates).
It should be noted that in case of using rigid-plastic models (just as the Taylor model), the Voigt hypotheses ( $\mathbf{l}=\mathbf{L}$$\mathbf{l}=\mathbf{L}$l=L\mathbf{l}=\mathbf{L}$\mathbf{l}=\mathbf{L}$ ) and the spin tensor $\omega$$\omega$omega\omega$\omega$, which was established by the Taylor rotation model, the shear rates at the current moment of deformation with a known number of active SSs can be determined from the system of equations:

where ${\mathrm{K}}_{\mathrm{A}}$${\mathrm{K}}_{\mathrm{A}}$K_(A)\mathrm{K}_{\mathrm{A}}${\mathrm{K}}_{\mathrm{A}}$ is the number of active slip systems. In this case, additional questions arise about the sufficiency of the number of equations for determining the stress deviator for $K_A\ne 5$$K_A\ne 5$K_(A)!=5K_{A} \neq 5$K_A\ne 5$, and the fulfillment of the conditions for consistency of the stress state with the yield conditions. However, this requires separate consideration, which goes beyond the scope of the present article.
Let us consider the first stage of the algorithm for implementing the mesoscale submodel at an arbitrary $n$$n$nn$n$-th time step. From the solution at the previous step, prior to the start of the $n$$n$nn$n$-th step (for time $t_n$$t_n$t_(n)t_{n}$t_n$ ), all mesoscale variables (accumulated shears, critical stresses, operating stresses, and orientations of the MCS of each crystallite) are known. Recall that all operations with tensors (tensors of stresses, elastic characteristics) are carried out in terms of the components of the MCS bases of the corresponding crystallites; shears, shear rates, and critical stresses are "linked" to slip systems, whose orientation dyads are also completely defined in the MCS. This allows us to move from co-rotational differentiation and integration to the corresponding usual operations.
In the examined crystallite, the belonging of the SS to the active ones and the linear independence of their orientation tensors are verified based on the stresses, orientation tensors and critical stresses determined at the time ${\mathrm{t}}_n$${\mathrm{t}}_n$t_(n)\mathrm{t}_{n}${\mathrm{t}}_n$ using the Schmid condition. The set of such SSs is designated as A, and their number - as $K_{\mathrm{A}}$$K_{\mathrm{A}}$K_(A)K_{\mathrm{A}}$K_{\mathrm{A}}$. If $K_{\mathrm{A}}\le 5$$K_{\mathrm{A}}\le 5$K_(A) <= 5K_{\mathrm{A}} \leq 5$K_{\mathrm{A}}\le 5$, then using the prescribed kinematic loading (the specified tensor $\mathbf{1}=\mathbf{L}$$\mathbf{1}=\mathbf{L}$1=L\mathbf{1}=\mathbf{L}$\mathbf{1}=\mathbf{L}$ ), the usual algorithm of the Lin model is implemented, i.e. the system of Eqns. (2) is solved and the transition to the next crystallite is carried out.

If the number of active SS in the crystallite under consideration is $K_{\mathrm{A}}>5$$K_{\mathrm{A}}>5$K_(A) > 5K_{\mathrm{A}}>5$K_{\mathrm{A}}>5$, the calculations for them are performed using the iterative procedure, which is described below. It should be noted that the dimension of the space of asymmetric deviators is 8 . This implies that the maximum possible number of linearly independent orientation tensors SSs, the totality of which can be considered as a basis for decomposing the plastic component of the displacement velocity gradient, is also does not exceed 8 . Accordingly, the vertices of the yield polyhedron which are formed by the intersection of hyperplanes corresponding to the independent orientation tensors of the SS should also be no more than the eighth order. The formation of vertices of higher order is improbable; this will require the fulfillment of special conditions of compatibility between the (linearly dependent) orientation tensors and the values of critical stresses. As a result, it is further assumed that the stress state in the case of activation of more than 5 SSs corresponds to the vertex of the yield polyhedron of the order of $K_{\mathrm{A}}$$K_{\mathrm{A}}$K_(A)K_{\mathrm{A}}$K_{\mathrm{A}}$ determined by crossing linearly independent hyperplanes by $K_{\mathrm{A}}$$K_{\mathrm{A}}$K_(A)K_{\mathrm{A}}$K_{\mathrm{A}}$. Thus, all these active SSs are linearly independent, and by virtue of this fact, any of their subsets are also linearly independent. At each iteration, the set A is divided into two linearly independent nonintersecting subsets D (with dimension $K_{\mathrm{D}}=5$$K_{\mathrm{D}}=5$K_(D)=5K_{\mathrm{D}}=5$K_{\mathrm{D}}=5$ ) and $\mathrm{U}(K_{\mathrm{U}}=K_{\mathrm{A}}-5),\mathrm{A}=\mathrm{D}\cup \mathrm{U},\mathrm{D}\cap \mathrm{U}=\varnothing$$\mathrm{U}\left(K_{\mathrm{U}}=K_{\mathrm{A}}-5\right),\mathrm{A}=\mathrm{D}\cup \mathrm{U},\mathrm{D}\cap \mathrm{U}=\varnothing$U(K_(U)=K_(A)-5),A=DuuU,DnnU=O/\mathrm{U}\left(K_{\mathrm{U}}=K_{\mathrm{A}}-5\right), \mathrm{A}= \mathrm{D} \cup \mathrm{U}, \mathrm{D} \cap \mathrm{U}=\varnothing$\mathrm{U}(K_{\mathrm{U}}=K_{\mathrm{A}}-5),\mathrm{A}=\mathrm{D}\cup \mathrm{U},\mathrm{D}\cap \mathrm{U}=\varnothing$. The rules for assigning SS to subsets D and U at each $q$$q$qq$q$-th iteration are heuristic in nature and are set by a researcher, for example, using a cyclic interchange or random searching through the systems. It is also possible to assign to the number ${\mathrm{D}}_{\mathrm{q}}$${\mathrm{D}}_{\mathrm{q}}$D_(q)\mathrm{D}_{\mathrm{q}}${\mathrm{D}}_{\mathrm{q}}$ the slip systems with maximum deviations of the shear stresses determined at the previous iteration from the critical ones at the end of the loading step, i.e.

To describe the iterative procedure, we introduce some additional notation:

The need for an iterative procedure for each crystallite is mainly associated with the elastic components of the velocity gradient ${\mathbf{z}}^{\mathrm{e}}$${\mathbf{z}}^{\mathrm{e}}$z^(e)\mathbf{z}^{\mathrm{e}}${\mathbf{z}}^{\mathrm{e}}$, which are unknown at the loading step. If ${\mathbf{z}}^{\mathrm{e}}$${\mathbf{z}}^{\mathrm{e}}$z^(e)\mathbf{z}^{\mathrm{e}}${\mathbf{z}}^{\mathrm{e}}$ could be determined by any of the known methods, then it would be possible to use a system of equations of the form (4) with $\mathbf{z}$$\mathbf{z}$z\mathbf{z}$\mathbf{z}$ replaced by ${\mathbf{z}}^{\mathrm{p}}=\mathbf{z}-{\mathbf{z}}^{\mathrm{e}}$${\mathbf{z}}^{\mathrm{p}}=\mathbf{z}-{\mathbf{z}}^{\mathrm{e}}$z^(p)=z-z^(e)\mathbf{z}^{\mathrm{p}}=\mathbf{z}-\mathbf{z}^{\mathrm{e}}${\mathbf{z}}^{\mathrm{p}}=\mathbf{z}-{\mathbf{z}}^{\mathrm{e}}$, which would allow us to directly proceed to system (2), eliminating Eqn. (2)॰ from it. However, ${\mathbf{z}}^{\mathrm{e}}$${\mathbf{z}}^{\mathrm{e}}$z^(e)\mathbf{z}^{\mathrm{e}}${\mathbf{z}}^{\mathrm{e}}$ is unknown at the beginning of the step, so that in the first iteration of each loading step, ${\mathbf{z}}^{\mathrm{e}}$${\mathbf{z}}^{\mathrm{e}}$z^(e)\mathbf{z}^{\mathrm{e}}${\mathbf{z}}^{\mathrm{e}}$ is assumed equal to the zero tensor.
Let us consider the iterative procedure for an arbitrary iteration q , assuming that all variables entering into the system of Eqns. (2) from the previous iteration are known. At each q-th iteration, the following operations are performed:

The values of tangential ${\tilde{\tau }}_{(\mathrm{q})}^{(k)}$${\tilde{\tau }}_{(\mathrm{q})}^{(k)}$tilde(tau)_((q))^((k))\tilde{\tau}_{(\mathrm{q})}^{(k)}${\tilde{\tau }}_{(\mathrm{q})}^{(k)}$ and critical ${\tilde{\tau }}_{c(\mathrm{q})}^{(k)}$${\tilde{\tau }}_{c(\mathrm{q})}^{(k)}$tilde(tau)_(c(q))^((k))\tilde{\tau}_{c(\mathrm{q})}^{(k)}${\tilde{\tau }}_{c(\mathrm{q})}^{(k)}$ stresses on the SS at the end of the step are determined based on the trial values of shear rates.
2. The set of active slip systems is divided into two subsets ${\mathrm{D}}_{\mathrm{q}}$${\mathrm{D}}_{\mathrm{q}}$D_(q)\mathrm{D}_{\mathrm{q}}${\mathrm{D}}_{\mathrm{q}}$ and ${\mathrm{U}}_{\mathrm{q}}$${\mathrm{U}}_{\mathrm{q}}$U_(q)\mathrm{U}_{\mathrm{q}}${\mathrm{U}}_{\mathrm{q}}$. In the present paper, random division of the set $A$$A$AA$A$ into the subsets $D_q$$D_q$D_(q)D_{q}$D_q$ and $U_q$$U_q$U_(q)U_{q}$U_q$ was used during the procedure implementation.
3. For the SS from the set ${\mathrm{U}}_{\mathrm{q}}$${\mathrm{U}}_{\mathrm{q}}$U_(q)\mathrm{U}_{\mathrm{q}}${\mathrm{U}}_{\mathrm{q}}$, the part of the inelastic component of the asymmetric measure of the strain rate is determined based on the trial shear rates of the current iteration:

where $\beta$$\beta$beta\beta$\beta$ is the relaxation coefficient used to improve convergence.
6. The parts of inelastic component of the asymmetric measure of deformation are determined for the SSs from the set ${\mathrm{D}}_{\mathrm{q}},{\mathrm{U}}_{\mathrm{q}}$${\mathrm{D}}_{\mathrm{q}},{\mathrm{U}}_{\mathrm{q}}$D_(q),U_(q)\mathrm{D}_{\mathrm{q}}, \mathrm{U}_{\mathrm{q}}${\mathrm{D}}_{\mathrm{q}},{\mathrm{U}}_{\mathrm{q}}$ :

To assess the computational efficiency of the modified elastoplastic Lin model, we performed a series of numerical experiments for complex and simple loading of a representative macrovolume (analogue to a macrosample of polycrystalline material) at a constant strain rate of $0.002{\mathrm{s}}^{-1}$$0.002{\mathrm{s}}^{-1}$0.002s^(-1)0.002 \mathrm{~s}^{-1}$0.002{\mathrm{s}}^{-1}$. The examined macrovolume consists of 343 equal volumes of grains with a uniform distribution of orientations and FCC lattice (aluminum). Numerical experiments allowed us to compare the computational efficiency of the elastoplastic (EP) and elastovicoplastic (EVP) models and the compatibility between the results obtained with the use of these models under loading regimes and for virtual samples with identical initial configurations. The main difference between the EVP model and the EP model was the replacement of relation (2). with Hutchinson's relation [9]:

where $m,{\dot{\gamma }}_0$$m,{\dot{\gamma }}_0$m,gamma^(˙)_(0)m, \dot{\gamma}_{0}$m,{\dot{\gamma }}_0$ are the model parameters, in the framework of this work ${\dot{\gamma }}_0$${\dot{\gamma }}_0$gamma^(˙)_(0)\dot{\gamma}_{0}${\dot{\gamma }}_0$ is taken equal to $0.00118{\mathrm{c}}^{-1}$$0.00118{\mathrm{c}}^{-1}$0.00118c^(-1)0.00118 \mathrm{c}^{-1}$0.00118{\mathrm{c}}^{-1}$, and the exponent varies in the range [30, 300]. The plan of numerical experiments is given in Tab. 1.
The parameters were taken from the article [25]:

To implement the above mentioned loading pattern, a kinematic method of specifying deformation in terms of velocity gradients is adopted,

The results of calculation are compared using the $L_2$$L_2$L_(2)L_{2}$L_2$ metric for grain displacements