Statistical Crystal Plasticity Model Advanced for Grain Boundary Sliding Description

: Grain boundary sliding is an important deformation mechanism, and therefore its description is essential for modeling di ﬀ erent technological processes of thermomechanical treatment, in particular the superplasticity forming of metallic materials. For this purpose, we have developed a three-level statistical crystal plasticity constitutive model of polycrystalline metals and alloys, which takes into account intragranular dislocation sliding, crystallite lattice rotation and grain boundary sliding. A key advantage of our model over the classical Taylor-type models is that it also includes a consideration of grain boundaries and possible changes in their mutual arrangement. The constitutive relations are deﬁned in rate form and in current conﬁguration, which makes it possible to use additive contributions of intragranular sliding and grain boundary sliding to the strain rate at the macrolevel. In describing grain boundary sliding, displacements along the grain boundaries are considered explicitly, and changes in the neighboring grains are taken into account. In addition, the transition from displacements to deformation (shear) characteristics is done for the macrolevel representative volume via averaging, and the grain boundary sliding submodel is attributed to a separate structural level. We have also analyzed the interaction between grain boundary sliding and intragranular inelastic deformation. The inﬂux of intragranular dislocations into the boundary increases the number of defects in it and the boundary energy, and promotes grain boundary sliding. The constitutive equation for grain boundary sliding describes boundary smoothing caused by di ﬀ usion e ﬀ ects. The results of the numerical experiments are in good agreement with the known experimental data. The numerical simulation demonstrates that analysis of grain boundary sliding has a signiﬁcant impact on the results, and the multilevel constitutive model proposed in this study can be used to describe di ﬀ erent inelastic deformation regimes, including superplasticity and transitions between conventional plasticity and superplasticity.


Introduction
Superplastic (SP) deformation is promising for manufacturing complex-shaped parts with the required properties from metals and alloys. Using this approach, one can produce lightweight large-size structures without welded seams or with a small number of welded seams and joints. Its application makes it possible to reduce the number of technological operations, to carry out the forming process with low stresses and to decrease the consumption of materials. The structures produced by this method have smooth surfaces and exhibit only slight surface deviations from the desired geometry, which is of particular importance for thickness variation minimization and for high-precision filling of dies [1][2][3][4][5][6][7][8][9][10].
In technological processes, the SP deformation regime develops at relatively moderate (about 0.5-0.6) homologous temperatures so that the size and equiaxed shape of individual grains remain unchanged [11][12][13]. At a higher temperature, the SP deformation can be realized under continuous dynamic recrystallization (DRX) conditions. However, in this case, expensive heat-resistant equipment is required, and it is almost impossible to guarantee the necessary continuity (absence of excessive porosity) of the complex-shaped parts during their manufacture. Besides, the residual stresses arising after high-temperature processing may be high, and the shape of an object may distort. Our study addresses the 1st type SP regime, called "structural SP with steady grain's state". This term has a more precise meaning compared to the pioneer works on superplasticity (e.g., [14]), in which only the average grain structure parameters were kept unchanged, and the material was redistributed among grains during the continuous DRX (the term "steady state structural SP" is used in the literature [15] in this broad sense).
In most works, the properties of the materials in the SP regime are determined on the basis of the results of uniaxial tensile tests. For many alloys, including those widely used in industry, within some temperature and strain rate ranges, there occurs a stress-strain dependence of staged nature, which can be attributed to the effects of different deformation mechanisms, their interaction with each other and a change in their roles during deformation. Similar dependences (bell-shaped curves) are given in the literature for aluminum alloys close to single-phase alloys [16][17][18][19][20], magnesium alloys [21,22] and titanium alloys [23,24]. These analyses reported the absence of a pronounced deformation localization zone (a neck), which allowed researchers to make an assumption about the uniformity of the stress-strain state.
In [25,26], basing on the analysis of information from most of the reviewed sources, the authors set forth a scenario that has been realized in deformation testing alloys close to single-phase alloys with a transition to the SP regime under prescribed temperature and strain rate conditions (the test temperature does not exceed 0.7 of the homologous temperature).
At the initial deformation stage (an ascending branch of the tension curve), the effects of intragranular dislocation sliding (IDS) and grain rotation mechanisms become stronger, which causes a crystallographic texture to occur, whereas the role of grain boundary sliding (GBS) remains insignificant.
At the transient stage (the tension curve gradually bends), the role of GBS increases. The grain boundaries are more prepared for this sliding because of the influx of lattice dislocations that increase the boundary internal energy and improve the smoothing of the grain surfaces by grain boundary diffusion. Depending on the initial structure of the material, test temperature and strain rate conditions, recrystallized grains begin to grow at the expense of the non-recrystallized phase and an integral decrease in the density of intragranular dislocations. During the deformation process, a grain rotation occurs. The GBS causes a shear of adjacent grains against one another.
The specific features of the structural SP regime (part of the curve that corresponds to a steady-state flow stress or its gradual decrease) are as follows. The GBS, accompanied by the accommodative processes of IDS and grain boundary diffusion, grain rotations and shear of adjacent grains against one another, is prevailed. The stability of the structure (all grains remain unchanged, and the structure becomes an equiaxial ultrafine-grained structure) is achieved via the simultaneous realization of accommodative IDS and grain rotation. GBS (as well as dislocation sliding at an transient stage) leads to weakening or disappearance of the texture formed after the first stage of testing.
Therefore, the available complex deformation scenario for the implementation of even the simplest (uniaxial at the macrolevel) SP tests includes several interacting mechanisms with changing roles and significant structure evolution of the material. A similar situation is typical for the technological processes involving SP. Temperature and strain rate conditions may change significantly in different parts of the structure, and thus these processes will proceed in different ways in these objects. This causes a need for the development of mathematical models for describing the peculiar features of existing metal-forming technologies. Much research has been published regarding the multilevel models developed within the framework of the crystal plasticity for describing IDS, twinning and rotation of crystallite lattices [27][28][29][30][31][32][33][34][35][36]. An important direction towards obtaining this class of constitutive models and an accurate description of traditional mechanisms (dislocation motion and interaction) in the context of solid-state physics includes a consideration of other significant deformation mechanisms. In order to describe the structural SP with steady grain's state and preparatory stage, it is necessary to constitute an improved single-crystal plasticity model capable of accounting for GBS, grain growth and some other deformation mechanisms mentioned above. GBS is the leading mechanism in the SP deformation and in the transition to it, and therefore its description is an issue of primary importance. However, to the best of our knowledge, this problem has not received due attention in the literature.
In [37], the activity of GBS was estimated in terms of the direct crystal plasticity model by analyzing tangential stresses on grain boundaries. In [38,39], where the deformation of grains was described with allowance for crystal plasticity, the grain boundary was modelled as a separate deformable medium (phenomenological relations were applied), but the grain sliding was not taken into account explicitly. In [40], almost the same approach is proposed, the mixture model with a smooth transition between phases is considered (the grain is described by the crystal plasticity model, the boundary is modeled by the J 2 -plasticity theory).
As noted above, the grain growth can play a significant role at the preparatory stage before the onset of SP regime. Recently, direct (based on solving the boundary volume problem at mesolevel via FEM) [41][42][43] and self-consistent crystal plasticity models [44,45] have been developed for describing the discontinuous DRX. However, it should be noted that these models include rather complicated procedures needed in grained restructuring and to redefine the stress field so that they satisfy a balance equation.
A self-consistent crystal plasticity model that is able to account for GBS and some other SP mechanisms has been introduced in a recently published article [46], the authors of which set forth an idea to describe "superplastic mechanisms (GBS, diffusion etc.)" in a complex manner by making use of an additional inelastic component for each grain. Although this approach shows some promise, it can be assumed that it hardly be used to formulate the kinetic equations needed to determine a given component via physical analysis (taking into account the interaction of GBS with another deformation mechanisms), also, it will lead to a significant complication of numerical procedures; besides, it is difficult to interpret the grain interaction effect using this approach. Direct crystal plasticity models are difficult to develop for GBS describing because of the need to describe the loss of continuity. For this reason, there occurs a necessity to apply special methods for solving boundary value problems. In recent years, some attempts have been made towards developing such mathematical apparatus (for example, in [47]), yet its successful practical application requires serious additions, including those promoting an increase in its computational efficiency. The issue of resource intensity is of great importance in solving boundary value problems applied to modeling SP technologies. Due to the totality of the arguments given, we have chosen a statistical model.
The structure of the model, in which part of inelastic strain rate due to GBS is separated at the macrolevel, is proposed in [48]. The models of similar structure are also used in [49,50], where the authors analyzed GBS in nanomaterials without considering physical processes that govern it (for instance, the equations for critical stresses were dropped from the consideration).
In [26], the first version of a multilevel statistical constitutive model for describing SP regime and transitions to it was developed. The model was formulated in the rate form in the current configuration, which made it possible to apply the hypothesis of additive contributions of IDS and GBS to strain rate at the macroscale level. In a series of works [51][52][53][54], the choice of this approach for formulating geometrically nonlinear equations was justified. It turned out that, in the absence of GBS, the model offers the results, which agree well with the data obtained using other mesoscale models, including those with the most popular formulation of constitutive equations in terms of the unloaded configuration [28,55]. In this paper, we describe in detail the modifications necessary to adapt the statistical crystal plasticity model to the consideration of GBS. In contrast to the model presented in [26], the structure of the model proposed here, as well as the kinetic equations used to describe SP and GBS effects, are changed. This modified model is a basis for developing a constitutive model advanced by an explicit account of discontinuous DRX (grain growth) and its effect on other deformation mechanisms, which are needed to describe both SP and transitions to it (an appropriate publication is scheduled for the near future). Since many deformation mechanisms operate and interact during the SP regime, we use here the simplest phenomenological relations (at the corresponding scales). Probable modifications of these relations are planned to perform in our further study on this topic. On the other hand, the presented constitutive model and numerical computations have revealed very important and interesting aspects regarding the crystal plasticity development for the simulation of deformation with active GBS. Section 2 includes a description of a three-level constitutive model and presents mesolevel equations for IDS and crystallite lattice rotations, as well as relations for GBS. Section 3 discusses the simulation results obtained in this study and gives their analysis.

Three-Level Model with Description of GBS, IDS and Grain Lattice Rotations
In contrast to the classical Taylor-type statistical models, the advanced model proposed in this study is able to account for the interfaces between crystallites. Intergranular boundaries are considered in the model as separate objects, connected at each moment of time with boundaries (grain facets). For these objects, internal variables are introduced, and kinetic equations for their evolution are formulated. A set of intergranular boundaries forms a separate structural level (level for description of GBS), and the model is modified into a three-level model: macrolevel (representative macrovolume), structural level for description of GBS that can be considered as a macrolevel submodel, and mesolevel (grains with boundaries). The schematic diagram of the model structure is given in Figure 1. At the initial moment, the aggregates of grains and intergranular boundaries are constructed, and the relations between them are assigned (for the intergranular boundary, the numbers of the grains adjacent to it are set). Due to GBS, the mutual arrangement of grains in the statistical model will change and this will gradually change the numbers of grains that form intergranular boundaries. The initial structure (samples of grains and intergranular boundaries) can be formed both coincident with spatial packing (e.g., the grain structure can be approximated by using Voronoi tessellation) and in a random way (due to the statistical equality of the model elements). The calculation results confirm the closeness of the integral parameters obtained at the macrolevel, including the change in the GBS proportion in the total strain rate and the change in the distribution function of grain lattice orientations, for different generations of the sample relations. This approach can be called a modification of the "statistical model which takes into account intergranular boundaries, mutual grain arrangement and its change". As regards the possibility of redefining the relations between the elements, it is similar to the movable cellular automata method [56].
For simplicity of the presentation, the formulations are given for isothermal deformations because a maintained constant temperature is observed in most tests and technological modes of SP forming. The changes necessary to consider the non-isothermal case are easy to introduce into the constitutive model.
We suppose that the macroscale stresses integrally characterize the elastic bonds in the material, and therefore the macroscopic stresses are equal to average mesoscale stresses where K, k is the weighted Kirchhoff stress tensor at macro-and mesolevels (the grain index is omitted), k = o ρ /ρ σ, σ is the Cauchy stress tensor, o ρ is the crystallite density in the reference configuration, andρ is the crystallite density in the current configuration. The continuity conditions are assumed to be fulfilled along the grain contact boundaries (no material discontinuities are present along the normal to the boundary).
The strain rate due to GBS Z in gb at the macrolevel is determined by averaging the GBS shear rates along the grain contact boundaries. Instead of standard (for the statistical Taylor type model) transmission of the velocity gradient (l = L) to the lower level, the following relation is used (2) where l, L is the transposed velocity gradient at meso and macrolevels, respectively. It can be said that some of the applied kinematic impacts are realized due to GBS at the macrolevel without contributing to the deformation of grains. Theoretically, there may be a situation where, after reaching a certain state of boundaries, the entire instantaneous deformation is realized due to GBS (among the grain layers), i.e., L = Z in gb , l = 0. This situation may occur in the bicrystal subjected to simple shear loading in the boundary plane. Relations describing GBS and an example that illustrates (2) are given below. First, let us consider the relations at the mesolevel.
To describe the deformation of each crystallite (the crystallite index is omitted), the following system of elastoviscoplastic equations is used In relations (3), summation over one and two contracted index is denoted by "·" and ":", the outer tensor product is denoted by "⊗", k cor ≡ . k + k · ω-ω · k is the rate of change (corotational derivative) of the Kirchhoff stress tensor, in which the principle of frame-indifference is satisfied, ω is the spin tensor of the rigid movable coordinate system (RMCS) related to the crystallite lattice, o is the RMCS orientation tensor, п is the fourth order tensor of crystallite elastic properties (its components are Crystals 2020, 10, 822 6 of 18 constant in the crystallite RMCS basis), z in is the inelastic strain rate component at the mesolevel; c are the cumulative shear, shear rate, effective tangential and critical tangential stress of the slip system k, b (k) , n (k) is the direction vector and normal to the plane of the slip system k, H(·) is the Heaviside function (hereinafter), and m is a dimensionless parameter. A doubled number of slip systems is used (only positive shears along them are considered). Note that in the constitutive equation, the transposed relative velocity gradient (l-ω) is used as a measure of the total strain rate minus the RMCS rotation [54]. Parameter . γ 0 is assumed to be related to the strain rate intensity . γ 0 = l u = (l-ω) : (l-ω) T (a similar relationship of the nominal shear rate due to IDS with macroscopic strain rate modulus is introduced in [57]).
It was shown in [53,54] that for small elastic deformations typical of metals and alloys and in the absence of GBS, the mesolevel model (3), describing the elastoviscoplastic deformation, yields results close to the data obtained by other alternative mesolevel models, including those with the most commonly used formulation of constitutive relations in the unloaded configuration [28,55]. However, use of the rate formulation (3) in the current configuration is preferable for the model development in the case of an increase in the number of considered mechanisms (in our case, we take into account GBS), since it is possible to accept the strain rate component's additive hypothesis.
The RMCS, in which the crystallite properties are assumed to be constant, is coupled with the symmetry elements of the crystallites, which remain unchanged during the IDS; the definition of the spin in the absence of GBS is discussed in detail in [51,54]. It was shown in [58] that the use of a full on grain rotation is demonstrated, e.g., in [59]. When GBS is activated along the grain boundaries, it becomes possible to rotate it additionally due to the action of shear stresses along the boundaries; this effect is described by the component added to the specified spin [26] here χ is the model parameter characterizing the viscous boundary resistance to rotations, µ τ is the moment due to the tangential forces acted on all grain boundary facets, z in gb grain is the mean strain rate due to the GBS ( grain denotes the averaging over the total number of boundary facets of the current grain), and the norm of the 2nd rank tensor T is defined as T = T : T T . An important characteristic of the mesoscale model is that it includes kinetic equations for critical shear stresses τ (k) c along the slip systems. Relations for describing the change τ (k) c have been proposed in many works, for instance, in [57,60,61]. The purpose of our research is to describe GBS and, as an example illustrating the application of this model, we consider the deformation in the SP regime. According to the above scenario, the grain structure in the SP regime can be considered as a recrystallized structure. In this case, the increase in τ (k) c due to the intragranular barriers can be neglected (in the general case, one can use the equations given in the articles already mentioned, or in [26], or in any other publications). On the other hand, since consideration is given here to a fine-grained structure, then it is reasonable that there occurs an increase in the resistance of dislocations to motion due to the boundary defects (orientation mismatch dislocations) and generation of near-boundary dislocation pile-ups during the misorientation between the neighboring grains [62,63].
In order to take into account the boundary effect on IDS, the shear rate relation used in direct crystal plasticity models is usually modified so that it can be applied to the near boundary region. In [64][65][66], the well-known relation [60] was modified to describe the evolution of the dislocation density over the slip systems, taking into account the pile-ups of dislocations near the boundary. In [67], for the near boundary region, the application of Arrhenius-type relation for IDS shear rates made it possible to analyze the additional activation energy proportional to the energy spent on orientation Crystals 2020, 10, 822 7 of 18 mismatch dislocations' generation; the shear rates in the near boundary area obtained in [68] turned out to be dependent on the mutual orientation of the slip systems and the shear stresses acting on them.
In accordance with the basic principles of constructing statistical models, the model proposed in this work effectively takes into account the rising of resistance to IDS due to orientation mismatch dislocations and dislocation piles-up near the boundary by increasing of the critical shear stresses τ (k) c for the grain (as a homogeneous element in a statistical sample). In the absence of GBS at the grain boundary, we apply the following approximate relation where N is the (outer) external normal to the current grain boundary, R (N) is the estimate (of density) of the rate of dislocation inflow into the outer boundary from it, R (N)neib is the estimate (of density) of the rate of passage of dislocations through the boundary into the neighboring grain. The index "neib" denotes the neighboring grain and ψ is the model parameter (Pa).
After the specified GBS occurs on the current grain boundary, a linear combination of the characteristics of the neighboring grains on this boundary in (5) is used, which is similar to the determination of shear stresses for the intergranular boundary carried out by (9). According to (5), in the case of the incompatibility of IDS in the neighboring grains, part of the lattice dislocations forms orientation mismatch dislocations in the boundary, which is taken into account through an increase in the resistance to IDS (in sight, the model can be expanded through the explicit introduction of variables for the densities of various types of defects).
With account for the possible activity of GBS, (5) takes the form where R GBS(N) is the estimate of the rate of increase in the grain boundary dislocation density, on which orientation mismatch dislocations are dissociated during GBS, β is the dimensionless parameter, is the total value of the shear rates along the grain facet with the normal N against its size (four possible directions of positive shear rates are considered), α the dimensionless parameter, and τ c0 is the minimum possible critical stress for a grain with a low defect density.
gb is the area of the boundary facet, v is the grain volume, v gb0 is the characteristic rate of grain boundary dislocation displacement along the facet, τ cgb are the tangential and critical stresses for the contact boundaries, n is a dimensionless parameter; the relations for determining the above quantities are considered below. Use of (6) makes it possible to explicitly take into account the fact that, during the GBS, some part of orientation mismatch dislocations dissociates in grain boundary dislocations, which results in a weakening of the grain boundary hardening or causes its softening (depending on the ratio of IDS and GBS rates).
According to the above relations, GBS does not occur for a fine-grained material with a structure prepared for SP at a low homological temperature, and therefore a significant grain boundary hardening may happen due to IDS, which is mainly caused by the boundary effects. At a temperature required for SP, GBS will be active and, according to (6), the critical IDS stresses will decrease. Thus, the model is evidently capable of taking into account the effect of GBS on IDS.
During the GBS, relative displacements of crystallites (grains, subgrains) occur along their common boundary due to the motion of grain boundary dislocations. Note that, in the absence of an explicit dislocation description of (as in the mesoscale model, describing intragranular deformation), the averaged characteristics (shears) are used. To determine the kinematic characteristics of GBS for each grain boundary, a basis of vectors n gb , b (1) gb иb (2) gb is introduced into a three-dimensional space: n gb is a normal to the facet plane, b (1) gb and b (2) gb are the mutually orthogonal unit vectors of displacement directions in the facet plane. The dyad of sliding direction and normal determines the GBS slip system. When constructing the model under study, a transition from the local displacements along the contact boundaries to the deformation characteristics is carried out for the macrolevel representative volume. The movement of grain boundary dislocations in the slip plane is considered, and the viscoplastic relation is used to set their average displacement rate. Thus, we write [69] Crystals 2020, 10, x FOR PEER REVIEW 8 of 18 According to the above relations, GBS does not occur for a fine-grained material with a structure prepared for SP at a low homological temperature, and therefore a significant grain boundary hardening may happen due to IDS, which is mainly caused by the boundary effects. At a temperature required for SP, GBS will be active and, according to (6), the critical IDS stresses will decrease. Thus, the model is evidently capable of taking into account the effect of GBS on IDS.
During the GBS, relative displacements of crystallites (grains, subgrains) occur along their common boundary due to the motion of grain boundary dislocations. Note that, in the absence of an explicit dislocation description of (as in the mesoscale model, describing intragranular deformation), the averaged characteristics (shears) are used. To determine the kinematic characteristics of GBS for each grain boundary, a basis of vectors gb n , When constructing the model under study, a transition from the local displacements along the contact boundaries to the deformation characteristics is carried out for the macrolevel representative volume. The movement of grain boundary dislocations in the slip plane is considered, and the viscoplastic relation is used to set their average displacement rate. Thus, we write [69] ( ) where is the characteristic displacement rate of a set of grain boundary dislocations along the entire boundary facet, gb is the examined facet boundary area, V is the macrolevel representative volume size, n is the parameter, and 4K gb is the general number of base displacement directions on all boundary facets in the representative volume at the macrolevel. Generally speaking, rigorous determination of v gb0 requires the introduction of another scale-structural level (microscale) in the model so that one can get an opportunity to explicitly consider the characteristics of the varying defect structure (change ρ (i) gb with account for the interaction with lattice dislocations). Such an in-depth description may be useful for further development of the model. Analysis of the results obtained indicates that at this stage of research we shall confine ourselves to considering a three-level model, which makes it possible to achieve the objectives formulated. Therefore, in order to describe the characteristic displacement rate of grain-boundary dislocations moved along the boundary facet v gb0 , we have taken an approximate relationship between the strain rate intensity l u = (l − ω) : (l − ω) T and the grain size r: v gb0 = l u r (similar relationship between the nominal shear rate due to IDS and the macroscopic strain rate modulus was introduced in [57]). In (7), use is made of the variables which characterize the tangential stresses τ  The inelastic component of the relative velocity gradient due to GBS Z in gb is determined by the shear rates (7) as where K gb is the number of intergranular boundaries considered in the model, for which internal variables are specified, and the corresponding grains are assigned (during GBS, these neighboring grains change by others). For illustration, we consider here a simplest test case, assuming that cubical grains and boundaries are identical. Suppose that, at the current time, deformation develops only at the expanse of GBS, i.e., due to shears along the boundaries indicated in blue in Figure 2. This situation may take place when the tangential stresses on the corresponding GBS slip systems are equal to the critical stresses, i.e., shearing occurs in the grain layers. According to (7), the shear rate in a representative volume generated by the shear along one boundary is calculated as Crystals 2020, 10, x FOR PEER REVIEW 9 of 18 where gb K is the number of intergranular boundaries considered in the model, for which internal variables are specified, and the corresponding grains are assigned (during GBS, these neighboring grains change by others). For illustration, we consider here a simplest test case, assuming that cubical grains and boundaries are identical. Suppose that, at the current time, deformation develops only at the expanse of GBS, i.e., due to shears along the boundaries indicated in blue in Figure 2. This situation may take place when the tangential stresses on the corresponding GBS slip systems are equal to the critical stresses, i.e., shearing occurs in the grain layers. According to (7), the shear rate in a representative volume generated by the shear along one boundary is calculated as where r is the grain size. We get the boundary area . This indicates that the described kinematic effect associated with a simple shear will develop at the expanse of GBS. The tangential stresses for GBS slip systems at the intercrystallite boundary, which, in the context of the accepted model representation, is assumed to be associated with one of the adjacent grains (the neighboring grains are replaced by others in this case), are determined as where gb σ is the characteristic stress for the boundary under study (index is omitted), max γ is the shear value at which the grains constituting the boundary are changed (loss of contact with the original neighboring grain A), σ , A σ , B σ are the stresses in the (base) grain coupled with the where r is the grain size. We get the boundary area S (i) gb = r 2 , the grain volume v = r 3 and the size of a representative macrovolume V = Nr 3 and write Crystals 2020, 10, x FOR PEER REVIEW 9 of 18 where gb K is the number of intergranular boundaries considered in the model, for which internal variables are specified, and the corresponding grains are assigned (during GBS, these neighboring grains change by others). For illustration, we consider here a simplest test case, assuming that cubical grains and boundaries are identical. Suppose that, at the current time, deformation develops only at the expanse of GBS, i.e., due to shears along the boundaries indicated in blue in Figure 2. This situation may take place when the tangential stresses on the corresponding GBS slip systems are equal to the critical stresses, i.e., shearing occurs in the grain layers. According to (7), the shear rate in a representative volume generated by the shear along one boundary is calculated as where r is the grain size. We get the boundary area In the case under consideration, a set of N GBS slip systems develops actively, and for each of these modes the shear indicated above will develop as well. Thus, and, if we assume that for active boundaries and the spin of RMCS for grains = ω 0, then we have in . This indicates that the described kinematic effect associated with a simple shear will develop at the expanse of GBS. The tangential stresses for GBS slip systems at the intercrystallite boundary, which, in the context of the accepted model representation, is assumed to be associated with one of the adjacent grains (the neighboring grains are replaced by others in this case), are determined as where gb σ is the characteristic stress for the boundary under study (index is omitted), max γ is the shear value at which the grains constituting the boundary are changed (loss of contact with the original neighboring grain A), σ , A σ , B σ are the stresses in the (base) grain coupled with the In the case under consideration, a set of N GBS slip systems develops actively, and for each of these modes the shear indicated above will develop as well. Thus, gb /τ (i) cgb n k 1 ⊗ k 2 and, if we assume that for active boundaries τ cgb , and the spin of RMCS for grains ω = 0, then we have Z in gb = l u k 1 ⊗ k 2 . This indicates that the described kinematic effect associated with a simple shear will develop at the expanse of GBS.
The tangential stresses for GBS slip systems at the intercrystallite boundary, which, in the context of the accepted model representation, is assumed to be associated with one of the adjacent grains (the neighboring grains are replaced by others in this case), are determined as where σ gb is the characteristic stress for the boundary under study (index is omitted), γ max is the shear value at which the grains constituting the boundary are changed (loss of contact with the original neighboring grain A), σ, σ A , σ B are the stresses in the (base) grain coupled with the boundary, grain A and grain B which replaces grain A on the boundary, γ gb is the shear which occurred till the current moment of time along all the four base directions in the boundary, γ gb = γ gb are the cumulative shears that occurred in the boundary till the current moment of time (in four different directions, b gb ). At γ gb = 0, the characteristic stress is defined by the average stress for the initially adjacent base grain and grain A, and at γ gb = γ max by averaging stresses for the base grain and grain B.
In determining the cumulative shear, the shear rate direction along the boundary (vector is prescribed, and it corresponds to the onset of the GBS along the facet under study; the spatial arrangement of crystallites is taken into account. The value γ max , characterizing the moment at which the grains are changed for current contact boundary, can be obtained as follows where l gb , S gb denote the boundary size in the shear direction and the boundary area, and V is the representative macrovolume size. After the total shear reaches this limiting value, the shear γ (i) gb for the intercrystallite boundary examined here is reset to 0. Formula (10) is simplified for the homogeneous grain structure: γ max = l gb S gb /V ≈ v/V ≈ 1/N, where N is the number of grains in a representative volume. The important element of the constitutive model is that it involves kinetic equations for the grain boundary critical shear stresses τ (i) cgb , i = 1, 2, . . . , 4K gb (to shorten the writing of formulas, the index of the number of the GBS slip system is omitted). In the general case, these relations should take into account an increase (decrease) in the grain structure resistance to GBS, the emergence of hard-phase particles to the boundaries during DRX and the mechanical smoothing of the boundaries; a version of these relationships is given in [26]. In this work, attention has been focused on the SP regime at the final stage of the SP deformation test, when the indicated factors can be neglected (it is assumed that the corresponding processes have been completed). Therefore, the GBS hardening is expressed as The constituent F characterizes a decrease in the critical stresses caused by an increase in the boundary energy due to the influx of intragranular dislocations. The term G characterizes the processes which are associated with grain-boundary diffusion and which lead to a decrease in critical stresses (diffusion-induced adjustment of the boundaries). The critical stresses cannot be smaller than a certain minimum level τ cgb * .
The rate of decrease in critical stresses caused by an increase in the boundary energy due to the influx of intragranular dislocations (this effect is considered, for instance, in [70,71]) is described by the following equation where j is the index denoting the neighboring crystallite adjacent to the boundary; as in (7), consideration is given to the initial neighboring grain A and grain B that replaces it, µ j is the weight coefficient equal to (γ max − γ gb )/(2γ max ) for grain A and γ gb /(2γ max ) for grain B (thus, in (12), the same linear combination as in (9) is used); . γ (k)( j) is the shear rate along the slip system k in the crystallite j adjacent to the facet under study, K j is the number of slip systems in the grain j, b (k)( j) is the unit vector of the slip direction k in the crystallite j, N (k)( j) is the outer normal to boundary facet plane for the crystallite j and g is the model parameter characterizing the relationship between the boundary energy and the critical shear stress of GBS. The last term describes a decrease in the boundary energy due to GBS, which causes the softening value to decrease, ν is the model parameter characterizing the corresponding effect, four possible directions of positive shear rates are considered; the value of F cannot be positive. The constituent for the rate of change in critical stresses caused by the processes associated with grain boundary diffusion is defined as [26] where gb is the normal stress in the shear direction along the boundary b (i) gb (in the area orthogonal to b (i) gb ), calculated by analyzing the stresses σ gb characteristic of the boundary under study (7), σ, c are the model parameters, U d is the activation energy for grain boundary diffusion, θ is the temperature, κ is the Boltzmann constant. Term (13) characterizes a decrease in the critical shear stress of GBS due to diffusion-induced boundary smoothing, facilitating the dissociation processes of the orientational mismatch dislocations into grain boundary dislocations and lattice dislocations, dislocation climb.

Results and Discussion
The model described in this paper was applied to simulate the uniaxial tensile tests in which a representative volume of polycrystalline aluminum alloy 1420 (Al-5.5%Mg-2.2%Li-0.12%Zr) was investigated.
Since in most experiments the rate of crosshead displacement is assumed to be constant [25], in our numerical experiment the decreasing strain rate corresponding to this case is set in the tension direction (along the OX 3 -axis of the laboratory coordinate system), i.e., D 33 = .
where D 33 is the strain rate tensor component, H 33 is the logarithmic Hencky strain measure component, D 0 is the initial strain rate, and t is the time. The rest of the components of the (transposed) velocity gradient, which determines the kinematic effect used in the model, are found from the condition of uniaxial stress state [72]. The experiments were performed at the constant temperature of 523 K, which corresponds to a homological temperature of 0.56 (it is assumed that the equilibrium exchange rate with the environment having this temperature is significantly higher than the rate of possible dissipation-related temperature change during inelastic deformation). Deformation was considered only at the third stage of the SP deformation test [16] in the deformation interval from H 33 = 1.1 to H 33 = 1.45. Therefore, in modeling the final stage of this experiment, we set D 0 = 3.33 · 10 −4 s −1 (for experiment with initial strain rate 0.001s −1 , these data were used for identification procedure) and D 0 = 3.33 · 10 −3 s −1 (for experiment with initial strain rate 0.01s −1 , these data were used for verification). It is worth mentioning here that at this stage of research we intend to create the structure of a crystal plasticity model capable of describing GBS. For simplicity of the structure representation, the model does not take into account grain growth; only SP regime (changes in grain sizes are ignored) is modeled. Thus, a comparison with the experimental curve is performed at the third stage of the experiment, where grain growth is neglected in contrast to the early stages.
Within the framework of the proposed statistical model, the representative volume of the alloy 1420 was described through sampling from 343 f.c.c. crystallites, each of which had 26 neighboring boundary facets. The model parameters used in the numerical experiments performed were the same: constants for an observer in a rigid moving coordinate system associated with a lattice, elastic crystallite modulus p 1111 = 106.8GPa, p 1122 = 60.4GPa, p 1212 = 28.3GPa; viscoelastic law parameters for IDS m = 10; spin parameter χ = 1.5(MN · m) −1 ; parameters of the hardening law due to the boundary effect ψ = 1.7N, β = 1, α = 0.1; parameters of the submodel describing GBS n = 5, g = 40MPa, c exp − U d kθ = 0.03MPa/s, σ = 35MPa, q = 3. The state of the material is determined by the values of the internal variables, which depend on the history of impacts. In our model, such variables are the critical shear stresses for IDS and GBS and the crystallite orientation tensors. In the numerical experiments, the deformation is considered only at the final stage of the SP deformation test. To take into account deformation at previous stages, we set different values for different strain rates: τ In both cases, a blurred stretch texture was taken as the initial texture. Figure 3 gives the dependences of the Cauchy stress tensor component Σ 33 at the macroscale on the component H 33 of the logarithmic strain measure and the corresponding experimental data [16] (the dependence of stresses on the linear strain measure component from [16] is transformed into the dependence on H 33 ) for two different initial strain rates. It is seen that the simulation results are in a good agreement with the experimental results. Some deviations in stresses can be attributed to the fact that the grains are intensively replaced due to GBS (Figure 4). It follows from Figure 4 that, as the deformation increases, the number of replaced neighboring grains increases in a monotone way. Figure 5 shows the inelastic strain rate values obtained due to GBS (Z in gb ) 33 . These values correspond to the current strain rate D 33 and are determined at equal intervals of the strain H 33 .

33
, corresponding to the current strain rate D 33 and determined at equal intervals of the strain H 33 for the initial strain rate: Note that the lines in Figure 5 are somewhat arbitrary: Z in gb changes at each step, the diagrams contain the values (Z in gb ) 33 averaged on some sequent steps made at equal intervals of the strain H 33 . These diagrams characterize a contribution of GBS to the strain rate at the corresponding moments. This contribution turns out to be great, which provides evidence that GBS plays an essential role in superplastic deformation. At lower strain rates (the time taken for achieving the prescribed deformation increases), the shear rates vary gradually, and a slightly larger contribution of GBS takes place (Figures 4  and 5). The latter can be explained by the more significant effects of diffusion mechanisms, and the model developed is able to take this fact into account.
It is known that in the SP regime the GBS effect results in blurring the texture formed at the initial stage of the SP deformation test [73,74]. Figure 6 presents the pole figures constructed for the initial state (H 33 = 1.1) and the state occurred after reaching the strain H 33 = 1.45 (the figures illustrate the simulations at the initial strain rate D 0 = 3.33 · 10 −4 s −1 ; the results obtained for the higher initial strain rate are in good qualitative agreement with the abovementioned data). The results obtained demonstrate that the model can be used to implement the deformation scenario described above at the third stage in the SP regime with a gradual decrease in the flow stress ( Figure 3). This regime is characterized by the predominance of the GBS mechanism (Figures 4 and 5); the equiaxial fine-grained structure is stable due to the simultaneous implementation of accommodative IDS and grain rotations, when the shape change of grains occurs at each moment of time, but, because of unceasing rotations, the resulting grain deformation is close to zero. GBS leads to blurring the texture formed after the first stage of testing (Figure 6).
At this stage of our investigation, the most important aim was to develop a general structure of the model, which should include all the key physical mechanisms of deformation and take into account their interactions. Within the framework of this structure, specific relationships can be modified on the basis of more detailed information from solid-state physics and through the use of the experimental data. The findings about the action and interaction of mechanisms, changes in the material structure, along with the data on the evolution of the stress-strain state at the macroscale level, are useful for parameter identification. Actually, a simple uniaxial loading at the macrolevel is a three-dimensional loading at the mesoscale, in which all the deformation modes are manifested themselves, and corresponding information makes it possible to use the results of simple tests to develop a model suitable for studying complex three-dimensional loads that can be observed in technological processes.
Note that the issue about the end of grain growth remains debatable. The scenario given above implies that the process of DRX is completed before the onset of the SP regime. This permits describing the SP regime (the final stage of the SP test) with the proposed constitutive model, which does not take into account grain growth. At the same time, the structure of the model, which involves the relations describing an increase in the initially (in the state before SP deformation testing) recrystallized grains due to the absorption of the material of nonrecrystallized grains, will allow one to consider other options, in particular, those including the continuation of the DRX until the end of the deformation process, to describe changes in the stress-strain state of the material during the entire SP test. The advanced model holds promise for simulating technological processes, in which the history of impacts (in particular, grain growth) may differ significantly in different parts of the products.

Conclusions
Even in the case of a simple uniaxial tensile test with access to the SP regime, a complex deformation process scenario takes place. According to this scenario, different physical mechanisms act and interact with each other, their significance changes and the material structure evolves greatly. Practically the same situation may arise in the technological processes which are based on the use of SP. This generates a need to develop mathematical models for describing changes in the material structure and the physical and mechanical properties of the material depending on its state. Recent studies show that a multilevel approach involving internal variables, physical theories of viscoplasticity, and an explicit description of the material structure and physical deformation mechanisms, holds the most promise. One of the main directions in this class model development is to expand the number of significant mechanisms to be taken into account.
In the investigation, we have modified the statistical crystal plasticity model for GBS describing. This model serves as the basis for an improved constitutive model that accounts for SP deformation and transitions to it. In contrast to classical Taylor-type statistical models, this model provides a consideration of the boundaries between crystallites as separate objects, but connected at each moment of time with boundaries (grain facets). For these objects, parameters characterizing their properties are introduced and kinetic equations are formulated. A set of grain boundaries forms a separate structural level (the level of GBS description), and the model is modified into a three-level one. The kinetic equations proposed for GBS include accounting for a decrease in the critical stresses due to an increase in the boundary energy caused by the influx of intragranular dislocations and due to the processes associated with grain-boundary diffusion (diffusion-induced adjustment of the boundaries). The model accounts for the fact that, due to GBS, some orientational boundary mismatch dislocations dissociate into grain-boundary dislocations, which leads, depending on the ratio of IDS and GBS rates, to the weakening of hardening or to softening for IDS. Besides, when GBS is active on the grain boundaries, the model includes an additional spin term due the action of shear stresses on the boundaries. The results of the numerical experiments are found to be in satisfactory agreement with the experimental data, the model reproduces the scenario of the deformation process at the final stage of the SP deformation uniaxial tensile tests.
The authors believe that the constitutive model proposed in this study and analysis of the numerical results has revealed some important and interesting aspects related to the development of crystal plasticity for modeling deformation with active GBS. On the basis of the proposed structure of a multilevel statistical constitutive model, an improved model is currently being developed by means of an explicit account of discontinuous DRX (grain growth) and its effect on other deformation mechanisms. The use of this advanced model will allow an accurate description of the SP regime, a transition to this regime and exit from it, which ensures correct modeling of technological processes (the history of impacts may differ significantly in different parts of the articles manufactured during these processes).

Funding:
The authors gratefully acknowledge the financial support from the Russian Science Foundation Grant no. 17-19-01292.

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