**1. Introduction**

In recent decades, extensive research efforts have been made to develop a multilevel approach to constructing constitutive models for metals and alloys via the introduction of internal variables [1–4]. The crystal plasticity models obtained within the framework of this approach provide an opportunity to explicitly describe changes in the structure of the material (its anisotropic physical and mechanical macro-properties) and deformation mechanisms [5–9]. This validates the application of these models for studying and improvement of materials processing technologies, as well as for solving the problems associated with the development of new functional materials.

Crystal plasticity models can be divided into three classes. In statistical models aimed at describing the behavior of representative macrovolumes, consideration is given to a sample of uniformly deformed grains (reviews of the models of this class are provided in [10–12]. The application of modified statistical models makes it possible to take into account the interaction of neighboring grains, e.g., during the explicit description of the motion of dislocations across the boundary [13]. Self-consistent models are used to investigate the behavior of uniformly deformed grains in a matrix representing the

**Citation:** Trusov, P.; Shveykin, A.; Kondratev, N. Some Issues on Crystal Plasticity Models Formulation: Motion Decomposition and Constitutive Law Variants. *Crystals* **2021**, *11*, 1392. https://doi.org/ 10.3390/cryst11111392

Academic Editor: Wojciech Polkowski

Received: 13 October 2021 Accepted: 11 November 2021 Published: 15 November 2021

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

**Copyright:** © 2021 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/).

effective characteristics of polycrystals. These models provide the fulfillment of equilibrium conditions; the effective characteristics are determined using statistical averaging methods. Reviews of the self-consistent models are given in [7,10,14]. Direct models based on the finite element method are recognized to be the most effective models for analyzing the nonuniform stress-strain state of grains (reviews of these models can be found in [7–9,15,16].

The crystal plasticity models mentioned above differ in the way of aggregation (combining the elements of a lower scale level into one element of a higher scale level), but all these models are based on the mesolevel relations used to describe crystallites, which led to the name "crystal plasticity".

In the last decade, numerous attempts have been made to develop crystal plasticity models based on generalized, mainly gradient, continua [17–20]. However, most authors assumed that the material under study is a simple (first-order) material [21]. In this case, the basic structure of the mesolevel constitutive model aimed at taking into account intragranular dislocation slip and lattice rotation is composed of the following equations: constitutive law for the relationship between stresses and elastic component of deformation, kinematic relations (including equations for determining the lattice spin), relation for the determination of the inelastic component of deformation and equations for the description of hardening.

The description of hardening (usually through the evolutionary relations for the introduced internal variables "critical shear stresses") is the important element of the model; there are many works devoted to the development of crystal plasticity in this direction (we mention here only some of the well-known works [22–24]). To date, other ratios of the basic structure of the model can be considered as the typical ratios, and therefore they are often dropped from the consideration in papers. However, researchers use several different types of basic formulations that differ in kinematic and constitutive relations. In our opinion, it is important to undertake a thorough analysis of these relations, which will enable their comparison in the context of geometric nonlinearity description and fulfillment of thermodynamic requirements. This is especially important given the fact that an intensively developing direction in crystal plasticity is related to the application of constitutive models for describing other, side by side with intragranular dislocation sliding, significant mechanisms and processes, such as twinning [8,25–27], phase transitions [8,28,29], grain boundary sliding [13,30–32], recrystallization [33–35]. In order to create advanced models, it is necessary to make a reasonable choice of a platform (the basic structure of the model).

In Section 2, popular formulations of mesolevel models are presented and analyzed. In the finite form, the models are formulated in the unloaded configuration, and in the rate form, in the actual configuration. The incorporation of a corotational derivative into the models formulated in the actual configuration means that the motion is decomposed into a rigid motion of the moving coordinate system and a deformation motion of the medium relative to it; the stress rate relative to the moving coordinate system is associated with the rate of this deformation. Section 3 considers the structure of kinematic and constitutive relations in the case when the component responsible for the rigid rotation of a moving coordinate system, relative to which the elastic distortions of the lattice are determined, is explicitly separated in the deformation gradient. An analytical comparison of different formulations and the results of numerical calculations are given in Section 4.

#### **2. Some Well-Known Variants of Base Mesolevel Relations**

In a number of works, linear constitutive relations with isotropic elastic law are used to describe small deformations [36]. However, since crystal plasticity models are usually focused on describing the behavior of materials at large displacement gradients, including changes in anisotropic properties at the macrolevel (texture formation), such approximation is rarely used. In addition, the application of anisotropic elastic law causes significant errors in the calculation of residual stress [37].

Let us consider the well-known variants of geometrically nonlinear anisotropic relations of the basic structure of crystal plasticity models, including the description of

intragranular dislocation slip and crystallite lattice rotation. Single and double contractions are denoted by "·" and ":", respectively, and the outer tensor product is denoted by "⊗".

#### *2.1. Relations in Terms of the Unloaded Configuration (U-Model)*

In many studies, the classical decomposition of the Kröner–Lee deformation gradient *f* is explicitly utilized [38–40]:

$$f = f' \cdot f'',\tag{1}$$

where *f e* , *f <sup>p</sup>* denote the elastic and plastic components of the deformation gradient.

The crystal lattice rotation is associated with the orthogonal vector *r<sup>e</sup>* from the polar decomposition *f <sup>e</sup>* = *r<sup>e</sup>* · *u<sup>e</sup>* = *v<sup>e</sup>* · *r<sup>e</sup>* , so the lattice spin is described by the relation:

$$
\overline{\boldsymbol{\omega}} = \dot{\boldsymbol{r}}^{\varepsilon} \cdot \boldsymbol{r}^{\varepsilon T}, \tag{2}
$$

where the dot indicates the time derivative.

Constitutive relations are formulated in the "classical" unloaded configuration found from the actual configuration by applying an affine transformation with (*f e* ) −1 :

$$k = \pi\_0 : \mathfrak{c}',\tag{3}$$

where *k* = *J*(*f e* ) <sup>−</sup><sup>1</sup> · *σ* · (*f e* ) <sup>−</sup>*<sup>T</sup>* is the second Piola–Kirchhoff tensor defined in the unloaded configuration (*σ* is the Cauchy stress tensor, *J* = det*f e* ), *c<sup>e</sup>* = 1/2 - (*f e* ) *<sup>T</sup>* · *f <sup>e</sup>* − *I* is the elastic Green-Lagrange strain defined in the unloaded configuration, *π*<sup>0</sup> = *πijmnk*0*<sup>i</sup>* ⊗ *k*0*<sup>j</sup>* ⊗ *k*0*<sup>m</sup>* ⊗ *k*0*<sup>n</sup>* is the elastic tensor, in which *k*0*<sup>i</sup>* denotes the crystallographic basis in the reference configuration.

The plastic component of the deformation gradient *f <sup>p</sup>* is obtained from the relation:

$$\dot{\boldsymbol{f}}^{\boldsymbol{p}} \cdot (\boldsymbol{f}^{\boldsymbol{p}})^{-1} = \sum\_{k=1}^{K} \dot{\boldsymbol{\gamma}}^{(k)} \overset{\boldsymbol{\varrho}^{(k)}}{\boldsymbol{b}} \overset{\boldsymbol{\varrho}^{(k)}}{\otimes \overset{\boldsymbol{\varrho}^{(k)}}{\boldsymbol{n}}},\tag{4}$$

where *<sup>o</sup> b* (*k*) , *o n* (*k*) are the unit vectors of the slip direction and the normal to the slip plane in the reference configuration, and *K* is the number of the slip systems in a crystallite. Here, positive and negative dislocations (doubling of slip systems) are considered separately, which makes the specification of hardening laws more variable, and therefore we make use of this technique in our study.

In most models, the shear rates . *<sup>γ</sup>*(*k*) on the slip systems are determined by viscoplastic relations [41–43]:

$$\dot{\gamma}^{(k)} = \dot{\gamma}\_0 \left( \frac{\tau^{(k)}}{\tau\_c^{(k)}} \right)^m H(\tau^{(k)} - \tau\_c^{(k)}),\tag{5}$$

where *<sup>τ</sup>*(*k*), *<sup>τ</sup>*(*k*) *<sup>c</sup>* are the resolved and critical shear stresses on the *<sup>k</sup>*-th slip system, . *γ*<sup>0</sup> is the shear rate in the slip system in the case when the tangential stress approaches the critical shear stress, and *m* is the strain rate sensitivity exponent of the material, and *H*(·) is the Heaviside function. There are some papers, the authors of which propose different ways to determine the resolved shear stress; this issue deserves a separate discussion because it lies, strictly speaking, beyond the scope of this paper. For definiteness, we suppose that *τ*(*k*) = *k* : *n*(*k*) ⊗ *b*(*k*) , where *k* is the weighted Kirchhoff stress tensor (*k* = *Jσ*), *b*(*k*) , *n*(*k*) are the unit vectors of the slip direction and the normal to the slip plane in the actual configuration.

In many works, Relation (5) are often supplemented with a term exp-− *<sup>U</sup> κθ* , where *κ* is the Boltzman constant, and *U* is the model parameter (activation energy). This allows one to consider the dislocation motion as a function of temperature. There is also an approach that involves determination of the temperature-dependent parameters . *γ*0,m[44] or formulations of corresponding relations for the critical stresses *<sup>τ</sup>*(*k*) *<sup>c</sup>* [24]. In any case, the initial critical stress values *<sup>τ</sup>*(*k*) *<sup>c</sup>* (0) should appear to be temperature-dependent.

Based on the literature review, we can conclude that the strain rate dependence of stresses in crystal plasticity can be taken into account in different ways. In most of the works . *γ*<sup>0</sup> = *const*, and therefore the strain rate sensitivity is defined by the parameter m. The value of . *γ*<sup>0</sup> should match the strain rate in order to avoid the sawtooth stress-strain response (curve) of the material during the numerical implementation of the model. It should also be noted that in the general case, there is an interaction between the parameters . *γ*<sup>0</sup> and m. At high m, the elastoviscoplastic model is almost similar to the elastoplastic model. Relation (5) provides a point in the stress space in a small neighborhood of the yield surface defined by Schmid's criterion: when the point is trying to move at great distance away from the neighborhood, significant relaxation of stresses is realized, and the point gets back to it. One can also use the relation . *γ*<sup>0</sup> = *ηdu*, *du* = <sup>2</sup> <sup>3</sup>*<sup>d</sup>* : *<sup>d</sup>* is the strain rate intensity, *d* is the deviator of the stretching tensor *d*, *η* is the model parameter (usually *η* = 1) [24,45]. For the active slip systems, the shear stresses tend towards critical stresses [45]. In the framework of this approach, the effect of the strain rate on the material response is embedded into the evolutionary relations for critical stresses *<sup>τ</sup>*(*k*) *<sup>c</sup>* .

The aim of this study is to analyze the kinematic relations and constitutive laws. To this end, we use, for definiteness, Relation (5) with constant . *γ*<sup>0</sup> in all cases under study and describe a hardening behavior according to a simple and well-known law [46–48] as:

$$\begin{aligned} \dot{\boldsymbol{\tau}}\_{\mathbf{c}}^{(k)} &= \sum\_{l=1}^{K} h^{(kl)} \dot{\boldsymbol{\gamma}}^{(l)}, \\ \boldsymbol{h}^{(kl)} &= \left[ q\_{l\,\,\rm at} + (1 - q\_{l\,\,\rm at}) \delta^{(kl)} \right] \boldsymbol{h}^{(l)} \boldsymbol{\lambda}^{(l)} = \boldsymbol{h}\_0 \Big| \boldsymbol{1} - \boldsymbol{\tau}\_{\mathbf{c}}^{(l)} / \boldsymbol{\tau}\_{\mathbf{s}\mathbf{a}t} \Big|^{a} \end{aligned} \tag{6}$$

where the latent hardening parameter *ql at* takes the values 1 and 1.4 for coplanar and noncoplanar slip systems (with numbers *k* and *l*), respectively, and *δ*(*kl*) is the Kronecker delta.

Thus, for definiteness, all the base models considered in our study involve Relations (5) and (6) so that the inelastic strain rate can be determined.

We will call the model, which uses the kinematic relations and constitutive laws (1)–(4), the "U-model" (the unloading configuration type of the model). The given base structure of the relations is employed in numerous works, for instance, in [8,47–53]. In [54], this model is termed the Kalidindi–Bronkhorst–Anand formulation [55].
