Some Issues on Crystal Plasticity Models Formulation: Motion Decomposition and Constitutive Law Variants

: In this paper, kinematic relations and constitutive laws in crystal plasticity are analyzed in the context of geometric nonlinearity description and fulﬁllment of thermodynamic requirements in the case of elastic deformation. We consider the most popular relations: in ﬁnite form, written in terms of the unloaded conﬁguration, and in rate form, written in terms of the current conﬁguration. The presence of a corotational derivative in the relations formulated in terms of the current conﬁguration testiﬁes to the fact that the model is based on the decomposition of motion into the deformation motion and the rigid motion of a moving coordinate system, and precisely the stress rate with respect to this coordinate system is associated with the strain rate. We also examine the relations of the mesolevel model with an explicit separation of a moving coordinate system and the elastic distortion of crystallites relative to it in the deformation gradient. These relations are compared with the above formulations, which makes it possible to determine how close they are. The results of the performed analytical calculations show the equivalence or similarity (in the sense of the response determined under the same inﬂuences) of the formulation and are supported by the results of numerical calculation. It is shown that the formulation based on the decomposition of motion with an explicit separation of the moving coordinate system motion provides a theoretical framework for the transition to a similar formulation in rate form written in terms of the current conﬁguration. The formulation of this kind is preferable for the numerical solution of boundary value problems (in a case when the current conﬁguration and, consequently, contact boundaries, are not known a priori) used to model the technological treatment processes.


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][2][3][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][6][7][8][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][11][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 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][8][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][18][19][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][23][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][26][27], phase transitions [8,28,29], grain boundary sliding [13,[30][31][32], recrystallization [33][34][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.

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 "⊗".

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][39][40]: where f e , f p denote the elastic and plastic components of the deformation gradient. The crystal lattice rotation is associated with the orthogonal vector r e from the polar decomposition f e = r e · u e = v e · r e , so the lattice spin is described by the relation: 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 : where k = J( f e ) −1 · σ · ( f e ) −T is the second Piola-Kirchhoff tensor defined in the unloaded configuration (σ is the Cauchy stress tensor, J = det f e ), c e = 1/2 ( f e ) T · f e − I is the elastic Green-Lagrange strain defined in the unloaded configuration, π 0 = π ijmn k 0i ⊗ k 0j ⊗ k 0m ⊗ k 0n is the elastic tensor, in which k 0i denotes the crystallographic basis in the reference configuration.
The plastic component of the deformation gradient f p is obtained from the relation: 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 . γ (k) on the slip systems are determined by viscoplastic relations [41][42][43]: . where c are the resolved and critical shear stresses on the k-th slip system, . γ 0 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 − U κθ , where κ is the Boltzman constant, and U is the model parameter (activation energy). This allows Crystals 2021, 11, 1392 4 of 19 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 τ (k) c [24]. In any case, the initial critical stress values τ (k) c (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 . γ 0 = const, and therefore the strain rate sensitivity is defined by the parameter m. The value of . γ 0 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 . γ 0 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 . γ 0 = ηd u , d u = 2 3 d : d 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 τ (k) c . 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 . γ 0 in all cases under study and describe a hardening behavior according to a simple and well-known law [46][47][48] as: where the latent hardening parameter q l 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.

Relations in Terms of the Actual Configuration with a Taylor Spin (T-Model)
The models are often formulated using the Jaumann type rate of the Cauchy stress. The law of elasticity is written in this case as [56,57]: where π is the elastic tensor in the actual configuration (defined for the current crystallographic orientation), d is the stretching tensor (presented on the basis of the additive decomposition d = d p + d e ), the plastic part of which is determined as: The spin ω T is calculated in the framework of the Taylor's rotation model [58]: where l is the velocity gradient. Some recent studies (e.g., [59,60]) have used the law of elasticity written as: Apparently, these studies suggest that any strain rate can be defined as .
ε, yet . ε = d remains the most common option. The . ε p is determined by making use of (8).
It is easy to see that, for . ε = d, Relation (10) can be written in the form: where π is the elastic tensor for establishing the relation between the stress rate and elastic strain rate shown above, π = Jπ. Relation (11) and equivalent Relation (10) are supposed to be used instead of (7) since the internal energy density (per unit mass) is defined as: According to (12), at small deformations (at ω T = 0 and for in case of change in volume) when (11) is fulfilled, there will be no energy dissipation during the purely elastic deformation, which does not provide (7). A detailed analysis of the fulfillment of the thermodynamic constraints is given in Section 2.3. Note that, due to the smallness of volume changes in metals and alloys, (11) and (7) give similar results.
The model that involves the kinematic relations and constitutive laws (8), (9), and (11) will be further called the "T-model" (Taylor spin model). The relations of such conceptual structure have been used in many papers, in particular, in [59][60][61]. Some variants of selfconsistent models [7] are based on the mesolevel equations of the same structure but with correction of spin (9) via subtraction of the antisymmetric part of the Eshelby tensor for considered grain.

Analysis of the Formulations with an Emphasis on Describing Geometric Nonlinearity and Fulfillment of Thermodynamic Constraints
In many works cited above, it was shown that, due to the structure of the model relations, the dissipation of energy under inelastic deformation is positive within the framework of both U-and T-models. The U-model elastic law (3) involves the energetic conjugate stress-strain measures, and this guarantees the absence of energy dissipation during purely elastic deformation. Besides, the finite form of the relations ensures that, after the purely elastic deformation of the material (being in the unloaded configuration) along any closed deformation trajectories, the stresses will be zero. Thus, the conservation conditions for purely elastic deformation [62] are fulfilled within the framework of the U-model. At the same time, the U-model implies the necessity to perform numerical integration of (4) so that f p can be determined. This can violate the condition of isohoricity of plastic deformation and require corrective techniques. Therefore, ref. [54] proposes an alternative formulation for crystal plasticity that is based on the Kröner-Lee decomposition (also) and logarithmic strain measure, whose rate was decomposed additionally by using elastic corrector rates (defined on the basis of the plastic flow characteristics-shear rates on the slip systems). A major advantage is the lack of a necessity to calculate the exponent when determining f p at the end of each step. Meanwhile, we need to point out here that the interpretation of the physical meaning of the elastic law when considering it as a relation between the introduced complex measures of the stress-strain state (not only in the context of thermodynamics) is rather complicated. In addition, the fulfillment of the additivity condition for the components of the strain rate measure is guaranteed only if one slip system is active; in the general case, this condition is fulfilled approximately [54].
On the other hand, for the study of technological processes of thermomechanical processing, it is more convenient to formulate the boundary value problem in terms of the actual configuration in the rate form. This enables one to use numerical methods: a step-by-step solution with a redefinition of the configuration of the computational domain (including the contacting surfaces) can be realized. Another significant advantage of this formulation is an additive decomposition of the rate of inelastic deformation into contributions from various mechanisms. So, advanced crystal plasticity models taking into account grain boundary sliding were constructed on the basis of this structure [13,[30][31][32]. The T-model, referring to this type of formulation, uses a transparent measure of the strain rate-stretching tensor d.
It should be noted that the use of the corotational derivative in Relation (11) suggests that the linear (verified) constitutive relation is assumed to be valid for the observer related to a moving coordinate system rotating with the Taylor spin [63]. This suggests that the motion is decomposed into the rigid motion of the moving coordinate system and the deformation motion with respect to it [64]. Note that, to be fully consistent with the issue, we should write the measure of the rate of elastic deformation on the right-hand side of (11) [65]. However, due to the tensor symmetry, the results will be similar to those obtained when using (d − d p ). In some studies, the authors explicitly say that they use linear constitutive relations written in the local coordinate systems of crystallites [66,67]. It is evident that, for the relations in terms of the actual configuration in the rate form, the key question is how to determine the spin of the moving coordinate system.
The determination of a hypoelastic law of the form (11), for which the above-mentioned conservation conditions are fulfilled under purely elastic deformation, is a well-known problem in the mechanics of solids [62]. The popular Zaremba-Jaumann [68,69] and Green-Naghdi [70] corotational derivatives do not guarantee the fulfillment of the conservation conditions. Therefore, a logarithmic corotational derivative without this drawback was developed [62,71,72]. However, in [63,73], one can find examples that illustrate the difficulties which can be encountered in modeling the multistage deformation processes (including unloading): if the purely elastic cyclic loading is realized upon completion of elastoplastic deformation and unloading, then the stress hysteresis will be observed. In [73,74], the logarithmic spin, called the kinetic logarithmic spin (k-logspin), was improved. The modernized spin makes it possible to meet the conservation conditions in the simulation of multi-stage deformation processes, including unloading, but the consideration is limited to anisotropic material. Thus, the problem of formulating a hypoelastic constitutive relation of the Form (11) for an anisotropic material that ensures the fulfillment of the conservative conditions has not yet been completely solved. Obviously, this is also the case for the crystal plasticity T-model. So, both the U-model and the T-model have "a certain extent roughness" formulations, which still arouses the interest of researchers in finding a more elegant solution, with the leveling of these "roughness" formulations (see, e.g., [54]).
Note that, when formulating the T-model, the moving coordinate system (MCS) must be related to the crystallographic coordinate system since the components of the elastic tensor are assumed to be constant in the MCS. In other words, the MCS must be associated with the symmetry elements of the crystal.
According to the Taylor constrained rotation model, the displacement gradient at the mesoscale (grain level) is represented in the form [58]: whence, by taking the antisymmetric part, we obtain the relation for the lattice spin (9). Traditionally, the Taylor model is used for rigid-plastic models (this follows from the symmetrized Relation (13)-all deformation is associated with the shears on the slip systems). The closeness of the spin (9) to the "material rotation" spin (2) was demonstrated [9,75]. In both cases, one cannot speak of a direct relation between the MCS and the symmetry elements of the lattice. We also note that during deformation, the crystallographic coordinate system (CSC) can be distorted, while the MCS remains rigid, and therefore, strictly speaking, their identification is an approximation of the mathematical model.
In [63,65,76], an approach to the formulation of geometrically nonlinear relations is proposed. The approach is based on a modified multiplicative decomposition of the deformation gradient with an explicit selection of rigid rotation of the MCS, which is associated with the symmetry elements of crystallites. In the next section, the equations of the corresponding mesoscale formulations are given in a brief form (in the finite form in the lattice unloaded configuration and in the rate form in the current configuration). In Section 4, these ratios are compared with the U-model and T-model and are utilized to compare these formulations to each other.

Relations in Terms of Lattice Unloaded Configuration (LU-Model)
In [63,76], it was proposed to modify the decomposition of motion at a mesolevel, i.e., to realize a multiplicative representation of the deformation gradient f with an explicit separation of the rigid moving coordinate system (MCS) Op 1 p 2 p 3 : Here, r =k i ⊗ k i t=0 ≡k i ⊗ k i is the proper orthogonal tensor which converts the reference basis MCS) k i into the current basisk i ; in other words, r is the tensor of rotation of the MCS from reference to actual configuration; f p is the component of the deformation gradient that transforms the initial configuration into the plastically deformed configuration, f e is the component of the deformation gradient component that transforms the plastically deformed configuration into the actual configuration (this tensor also characterizes crystalline lattice distortion), and f e , f p are, in the general case, the non-symmetric second rank tensors. The motion of the MCS is a rigid motion, and the motion relative to the MCS is a deformation motion: geometrically generalized nonlinear constitutive relations are formulated by an observer related to the MCS [65]. Relation (14) is a modification of the classical Kröner-Lee decomposition (1) via the explicit separation of the quasi-rigid motion r. Note that, in the general case, the tensor does not coincide with the orthogonal tensor r e in the polar expansion f e = r e · u e = v e · r e in the case when the spin of the corresponding MCS is defined (2).
For analytical calculations, we also introduce the crystallographic coordinate system (CCS) Oy 1 y 2 y 3 with basis q i . Without loss of generality, we assume that, in the reference configuration, the CCS basis q i (0) is an orthonormal basis that coincides with k i . Note that the lengths of the CCS basis vectors and the angles between them vary during deformation.
The constitutive laws are formulated in the context of the observer related to the MCS. This MCS is assumed to be linked to the symmetry axes of the material. Probably, the idea of the necessity of determining corotational derivatives with reference to the symmetric elements of the material (directors) was initially put forward by J. Mandel [77], yet, unfortunately, his study contains no specific relations for describing a crystallite lattice spin. Based on the introduced multiplicative decomposition (14) and assuming that the elastic properties of the crystallite in MCS are constant, the constitutive elastoviscoplastic relation was formulated in terms of the unloaded configuration (unloading is carried out , and the MCS remains fixed) [63,76]: which includes the second Piola-Kirchhoff tensor k, defined in the unloaded configura- and the four-rank elastic tensor π = π ijmnk i ⊗ k j ⊗k m ⊗k n , c e is the elastic constituent of the right Cauchy-Green deformation tensor , and I is the unit tensor. The accepted definition of the elastic tensor (its components are constant in the MCS) conforms to the principle of independence of constitutive relations from reference choice [21]. The plastic component of the deformation gradient f p is determined using the viscoplastic or elastoplastic relations applied to describe the intragranular dislocation slip. In the model, as in other compared models, we use Relations (4)- (6).
An important component of the constitutive model is the kinetic relation for the motion of the MCS. As noted above, it seems reasonable to relate the mesoscale MCSs and the symmetry elements of the crystal lattice [65,76,77]. Plastic deformations occur due to the movement of dislocations and they do not cause symmetry distortion, and therefore r can be completely determined by the tensor f e .
The elastic gradient f e = f e · r contains both the quasi-rigid motion (rotation) of the MCS and the distortion of the lattice relative to it. Generally speaking, there are different variants of this relation. In order to overcome this challenge, it is necessary to accept a hypothesis for the definition of r. We consider here one of the options. Let a relation between CCS and MCS be as follows: (1) the axes Oy 1 and Op 1 coincide at every instant of deformation (the vectork 1 is directed along the vector q 1 ); (2) the vectork 2 is located in the plane Oy 1 y 2 at each instant of deformation. In [65], the following relation for the MCS spin is presented: k i ⊗k i = − k 2 · l e ·k 1 k 1 ⊗k 2 + k 2 · l e ·k 1 k 2 ⊗k 1 − − k 3 · l e ·k 1 k 1 ⊗k 3 + k 3 · l e ·k 1 k 3 ⊗k 1 − − k 3 · l e ·k 2 k 2 ⊗k 3 + k 3 · l e ·k 2 k 3 ⊗k 2 .
where l e = .
f e · f e−1 is the elastic component of the velocity gradient. In [78], the authors described the way to define r in finite form taking into account the above relation between the CCS and the MCS and showed the compliance of the obtained result with that obtained by integrating (16).
Besides, the MCS can be related to other symmetry elements of the crystal lattice. As an alternative, we consider the case when the axes Oy 3 and Op 3 coincide at each instant of time (the vectork 3 is directed along the vector q 3 ), the vector k 1 is positioned at each instant of deformation in the plane Oy 1 y 3 . Under these conditions, the spin can be calculated: We call the modified version of the LU-model with spin (18) the "LU2-model".
Consider the issue regarding the difference between LU1-model and LU2-model. We assume that, at equal kinematic effects f , the current positions of the MCS related to the symmetry axes of the lattice in different ways can be described by different tensors r and r * (due to the smallness of elastic distortions, these tensors are similar). So, the deformation gradient can be written as f = f e · r · f p = f e * ·r * · f p * .
Based on the constitutive Relation (15), in the case when the orientation of the MCS is defined via the application of the tensor r, the Cauchy stress are established as: where π = π ijmnk ikjkmkn = π ijmn (r · k i )(r · k j )(r · k m )(r · k n ). Relation (19) is easy to rearrange [80] into the form: Analogously, when the position of the MCS position is defined via the application of the tensor r * , the Cauchy stress can be obtained from the relation: σ * = 1/J * f e * · π ijmn k i ⊗ k j ⊗ k m ⊗ k n : 1/2 f e T * · f e * −I · f e T * (21) Comparative analysis of (20) and (21) shows that the difference in the definitions of stresses obtained using the LU1-and LU2-models can be attributed to the fact that f e and f e * are determined in different ways in case of inelastic deformation. As noted above, the definition for shear stress can be expressed as τ (k) = k : n (k) ⊗ b (k) . The unit vectors of the slip direction and the normal to the plane are found in terms of the actual configuration . It is evident that, at the elastic deformation stage and in a transition to plasticity (as f e = f e * ), the crystallographic characteristics will be the same, and therefore the reasons for the differences between f e and f e * that may occur are absent.
Thus, the mesostresses obtained within the framework of both LU1-and LU2-models are equal. The final orientations of the lattices determined using the models are different, though due to the smallness of elastic distortions, they differ only slightly [80]. Since the LU1-and LU2-models give almost the same definitions of stresses, we will call them an "LU-model" (if the case in point concerns the definition of stresses).
It is easy to show that the constitutive relation of the LU-model is independent of the imposed rigid motion [21]. Actually, let's consider deformation gradient is the orthogonal tensor describing the rigid rotation added to the initial motion. For pure elastic deformation, we get f e = O T · f e , and the analysis of (20) . Thus, measures of the stress-strain state with an imposed rigid movement are related by exactly the same function. It means that the principle of independence of the constitutive relation from the imposed rigid motion is fulfilled. Hence, it follows that at the onset of the plastic flow f p = f p is correct, and, as a result, the principle of independence of the constitutive law from the imposed rigid motion (i.e., σ = σ( f e ) will also be fulfilled in the case of inelastic deformation as well.
By virtue of the fact that the used stress and strain measures are energy conjugate, the requirement for the absence of stress hysteresis and energy dissipation in the arbitrary closed cycles of elastic deformations is fulfilled automatically in the LU-model [63], as in the U-model.
In the general case, other physically justified rotation models (in particular, those capable of taking into account the interaction of defects in neighboring crystallites [81] or the contribution of grain boundary sliding [13]) can be incorporated into the proposed approach to determine the MCS spin.

Rate Form Relations in Terms of the Actual Configuration (LR-Model)
From the formulation of the LU-model with the constitutive law in finite Form (15) through the corotational differentiation [64,76], we can pass to the rate analogue for the observer in the MCS: where c e − ω 1 · c e + c e · ω 1 (for implementation of this transition, we use the property of constancy of π for the observer in the MCS).
It was shown in [63,82] that, at f e ≈ I which is characteristic of metals, the following formulation in the rate form can be derived in terms of the current configuration.
Since, at f e ≈ I, the unloaded and actual lattice configurations are close, the following estimate is valid for the velocity gradient: where b (k) , n (k) are the vectors of the slip direction and the normal to the slip plane in the unloaded lattice configuration, i.e., they are in the actual position of the MCS. If the unloaded and actual lattice configurations are in close proximity, then we can assume b (k) = b (k) , n (k) = n (k) for the latter relation and use as an approximation for the inelastic strain rate. According to (23), the clear motion decomposition is established in the rate form: the velocity gradient l is represented as a set of the MCS spin and the rate of elastic distortions and inelastic deformations determined by the observer related to this MCS.
Taking into account the described approximations, the constitutive law (22) is close to where k is the weighted Kirchhoff stress tensor, and k cr = . k + k · ω 1 − ω 1 · k is its corotational derivative. Let us call the model involving the kinematic Relations (8), (17), and (23) and a constitutive law (24) the "LR1-model" (lattice rate type of model with spin 1). We also consider an analogous LR2-model in which Relation (18) is used for the spin. It was shown in [80] that the Cauchy stress determined using the formulation of this kind depends only slightly on the way they are related to specific material symmetry elements.
From the above discussion, we can draw a conclusion that the formulations of the Uand LU-models can provide the fulfillment of the conservative conditions for elastic deformation. On the other hand, for the study of technological processes of thermomechanical treatment with an a priori unknown configuration of the computational domain (including contact with the tool), the formulations in terms of the actual configuration in the rate form offer the most promise. Another advantage of such formulations is the possibility of an additive decomposition of the inelastic deformation rate into the contributions of different deformation mechanisms. This explains why the T-model is also popular.
The formulation based on the decomposition of motion with an explicit separation of the motion of the moving coordinate system (LU-model) makes possible a theoretically substantiated transition to a similar formulation in the rate form in terms of the actual configuration (LR-model). As will be illustrated below, this helps to easily show that all the models considered here give a similar response under equal influences.

Results and Discussion
In this section, the issue regarding a comparison of the formulations for crystal plasticity models described above is considered.

Analytical Comparison
Based on the Equation of motion (14), constitutive law (3) for the U-model can be written as: where k i is the reference basis of the MCS related to the lattice symmetry. The rearranged form of Relation (25) is: Relation (26) can be represented as: Relation (27) corresponds to the constitutive law of the LU-model (15). Thus, the stresses in the U-and LU-models are equally determined. As shown in Section 3.1, the determined stresses are independent of the portion of the MCS motion in the deformation gradient. When we use models of this kind, it is important to keep in mind this distinction and verify separately texture description.
If, by analogy with the transition from the LU-model to the rate LR-model described in Section 3.2, we take for the U-model an MCS that moves with spin (2), then we get the rate formulation of the LR-model type, but with a corotational derivative of the Green-Naghdi type (the GN-model used, for instance, in [75,83,84]). In line with this assumption, the proposed model will give stresses close to those determined in the U-model.
In [75], the following relation between spins (2) and (9) is given: where B is the four rank elastic compliance tensor, d p is defined using Formula (8). According to (28), by virtue of the smallness of elastic deformations (B : σ), the use of spin Relations (2) and (9) in the mesolevel model will cause slightly different stresses to occur (in other words, the proximity of the T-model and the GN-model is their characteristic feature). Hence, we can suggest that at small elastic distortions, all the models examined should give a close response (stresses). Table 1 presents information obtained when performing the analytical comparison of determination of stresses within the framework of different models.
In order to confirm the analytical results summarized in Table 1, we performed numerical calculations. Table 1. Comparison of determination of stresses within the framework of different models (notation used here: «=»precise coincidence of stresse obtained in the indicated models, «≈»-stresses are close; by "described" is meant the description of the essence of calculations with references to publications).

Illustrative Numerical Examples
We now turn to the consideration of the results obtained via the analysis of the application of kinematic relations and constitutive laws in the statistical constitutive model for simulating some kinematic loads imposed on the fcc polycrystal.
The representative volume included a sample of 343 crystallites, the initial orientations of which were distributed randomly according to a uniform law. The nominal properties of the polycrystal corresponded to those of copper. The mesolevel elastic property tensor contained the following independent components (constant for the observer in a rigid moving coordinate system linked with the lattice): π 1111 = 168.4 GPa, π 1122 = 121.4 GPa, π 1212 = 75.4 GPa [85]; the viscoplastic relations included . γ 0 = 0.001 s -1 , 1/m = 0.012; the hardening law parameters were h 0 = 180 MPa, a = 2.25, τ sat = 148 MPa, and the initial values of critical stresses for all slip systems were τ (k) [46,47]. Due to the significant nonlinearity of the systems of equations for all considered models, in particular, because of the presence of the Heaviside function in the viscoplastic law (5), which led to the necessity of discretization with a small time step to trace the activity of slip systems in crystallites, these equations were integrated using an explicit Euler method (in the calculations, the time step was 0.002 s). The results obtained with the aid of the LR1 model correspond satisfactorily to the experimental data for uniaxial compression and shear; the comparison of the results is given in [63].
We consider the affine deformations of the sample (corresponding to the representative volume of a polycrystal) subjected to uniformly distributed kinematic loads. Due to the uniform deformation, the radius vector of the material point at an arbitrary instant of time f (t) = 1 1+4h sin(ωt) p 1 p 1 + (1 + 4h sin(ωt))p 2 p 2 + p 3 p 3 + +3h(1 − cos(0.7ωt))p 1 p 2 , where ω = 0.001π, h = 0.1 is the constant parameter (deformation is considered within the time interval t = [0, T], T = 1000 s). For illustration, the trajectory of the point with Lagrangian coordinates (0, L, 0) m is given in Figure 1.
tive volume of a polycrystal) subjected to uniformly distributed kinematic loads. Due to the uniform deformation, the radius vector of the material point at an arbitrary instant of time t is determined according to ( ) = ( ) i is the deformation gradient, i q stands for the Lagrangian coordinates of the considered point, and i p is the basis of the fixed laboratory coordinate system. The motion is determined by the deformation gradient (the choice of motion can be arbitrary):  (0.7 )) , where ω = 0.001π, h  is the constant parameter (deformation is considered within the time interval ], s t      ). For illustration, the trajectory of the point with Lagrangian coordinates (0, L ,0) m is given in Figure 1.  We note the proximity of the results obtained (the curves plotted for the varying macrostress components are visually almost indistinguishable for all models examined in this study). For numerical assessment of the deviation in the results, we take the U-model as a base model and introduce the norm for comparing the results obtained using other models:   We note the proximity of the results obtained (the curves plotted for the varying macrostress components are visually almost indistinguishable for all models examined in this study). For numerical assessment of the deviation in the results, we take the U-model as a base model and introduce the norm for comparing the results obtained using other models: We note the proximity of the results obtained (the curves plotted for the varying macrostress components are visually almost indistinguishable for all models examined in this study). For numerical assessment of the deviation in the results, we take the U-model as a base model and introduce the norm for comparing the results obtained using other models: in the space C n L 2 continuous at t ∈ [0, T] on the vector-function dimension n; in the calculations n = 9 (the number of the macrostress tensor component Σ that is equal to averaged mesostresses in crystallites). Table 2 summarizes the results for deviations in the macrostresses obtained by the U-model from those determined by other models. Taking into account that the norms are integral over a sufficiently long period of time (T = 1000 s), the results summarized in Table 2 demonstrate the proximity of stresses obtained via the use of different formulations of mesolevel models, supporting thus the efficacy of theoretical statements given above (Table 1). Similar textures are obtained, which is indirectly confirmed by the proximity of macrostresses for the considered complex loading path. In [80], the issue of comparing the results of using some different spins (for other complex loading paths) was considered with an explicit analysis of the simulated lattice rotations. Figure 3 shows the time dependence of the mesolevel stress components for a single crystallite (chosen arbitrarily). The initial orientations of the crystallographic and moving coordinate systems coincide and can be determined through the sequential rotation of the crystallographic system initially coincident with the coordinate system for the fcccrystal around the Ox 1 axis with an angle φ 1 = 4.078, around the Ox 2 axis with an angle φ 2 = 0.0739, around the Ox 3 axis with an angle φ 3 = 2.149 for the U-model (Figure 3a) and the LR1-model (Figure 3b).  We note the proximity of the results obtained (the curves plotted for the varying macrostress components are visually almost indistinguishable for all models examined in this study). For numerical assessment of the deviation in the results, we take the U-model as a base model and introduce the norm for comparing the results obtained using other models:  Figure 3 demonstrates that the mesostresses for the crystallite under study, which were obtained within the framework of the models mentioned above, are similar. Table 3 gives the norms for deviations in the mesostresses

Analysis of the results shown in
for the considered models, which are similar to those for deviations in macrostresses (30). The results shown in Figure 3 and in Table 3 indicate the proximity of mesostresses. This is typical for the overwhelming number of crystallites (i.e., for almost all random orientations). In special cases, deviations of the mesostress components (difference modulus) may, at a certain instant of time, exceed 50 MPa. This can be attributed to some specific features in the motion of a representative point in the stress space near the yield surface [86] in transitions between the neighborhoods of the vertices. After transitions, in a great number of crystallites, the stresses become similar (there is a separate publication devoted to the study of these special cases). Due to the small fraction of such crystallites (37 crystallites out of 343) and averaging, the macroscale stresses obtained on the basis of the statistical constitutive model data using mesolevel models differ insignificantly ( Table 2).
The differences in stresses obtained at the macrolevel are significantly lower than the experimental data scatter for samples from one production batch. Indeed, for each material, there is no complete identity of samples even from the same batch due to the presence of many random factors during manufacturing, which leads to differences in the microstructure (distributions of crystallite orientations, dislocation densities, etc.) and, as a consequence, to different experimental dependences of stresses on deformations. For example, such scatters for copper alloy samples are given in [87,88].
It should be noted here that the U-model and LU-model offer the most promise for modeling purely elastic behavior (e.g., elastic cyclic loading during exploitation) since the constitutive law is presented in the finite form, and there is no stress hysteresis in closed elastic cycles. Moreover, due to the use of energetically conjugate stress and strain measures, no energy is dissipated during the elastic deformation.
On the other hand, crystal plasticity models are primarily used to describe large inelastic deformations and material structure evolution. In this regard, formulations in terms of the actual configuration hold more promise because the physical meaning of the stress tensor is clearer in this case, and, consequently, the simulation of hardening using evolutionary equations is facilitated. These arguments lend additional support to using this formulation in case of the boundary value problem and in some advanced constitutive models at different deformation mechanisms. The theoretical predictions and the results of the numerical experiments indicate that the response produced by these models is similar to that obtained using the U-model. So, it turns out that the formulation in rate form would be most suitable for constructing advanced constitutive models (e.g., with complex hardening laws). The point is that the hypotheses adopted for constructing the evolutionary equations of the model via the use of the internal variables determined in the unloaded configuration can give more significant errors than the small deviations from the strict fulfillment of conservative conditions for the formulations written in terms of the actual configuration.
Geometrically nonlinear crystal plasticity models are applied to anisotropic crystallites, and therefore the decomposition of motion, used in rate formulations in terms of the current configuration, must be constructed by taking into account the symmetry of the material. In other words, the moving coordinate system determining the rigid motion in motion decomposition must be associated with the elements of material symmetry so that its symmetry properties can be taken into account correctly. Above, some options for constructing constitutive models with the introduction of such moving coordinate systems have been considered.

Conclusions
We have analyzed the most popular relations of crystal plasticity models in the context of geometric nonlinearity description and fulfillment of thermodynamic constraints in case of elastic deformation. These models are focused on the description of significant inelastic deformations and material structure evolution caused by these deformations. That is why the formulations in the rate form written in terms of the current configuration seem preferable to constructing advanced constitutive models able to take into account different deformation mechanisms and their complex interactions. The main advantages of these formulations over the formulation in the finite form are as follows: • clear physical meaning of the stress tensor, which simplifies simulations with evolutionary hardening equations; • use of the clear measure of the strain rate (the stretching tensor) and the possibility of an additive decomposition of the strain rate into contributions from various mechanisms; • ease of their use in the rate formulation of the boundary value problem under varying, a priori unknown, contact boundary conditions.
The use of a corotational derivative in such a formulation means that the model is based on the decomposition of motion into the deformation motion and the rigid motion of a moving coordinate system, the rate of stress change relative to which is associated with the strain rate.
We have obtained the relations of the mesolevel model with an explicit separation of the motion of the moving coordinate system and elastic distortion of crystallites with respect to this system in the deformation gradient. The MCS is related to the elements of the material symmetry of the anisotropic crystallite. The proposed formalism with an explicit consideration of the MCS allows us to reasonably pass from the formulation in terms of the unloaded lattice configuration in the finite form (where the thermodynamic constraints for purely elastic deformation are strictly fulfilled) to the close formulation in rate form written in terms of the current configuration.
We have compared these relations with a popular crystal plasticity formulation, which makes it possible to establish their proximity to one another. The results of the performed analytical calculations show the equivalence or similarity (in the sense of the response determined under the same loads) of the formulations under consideration. These conclusions were supported by the results of the numerical calculations. It should be noted that the proposed approach aimed to determine the spin of the MCS also could include other physically justified rotation models, e.g., those taking into account the interaction of defects of neighboring crystallites or the contribution from the mechanism of grain boundary sliding. Such models cannot be introduced into other known approaches to constructing constitutive relations via the use of kinematic variables without explicit consideration of the material symmetric properties.

Funding:
The study was carried out with financial support from the Ministry of Education and Science of the Russian Federation as part of the implementation of the national project "Science and Universities" (the state task fulfillment in the Laboratory of Multilevel Structural and Functional Materials Modeling).

Conflicts of Interest:
The authors declare no conflict of interest.