{'fraction_non_alphanumeric': 0.034150557198089605, 'fraction_numerical': 0.0144280191039345, 'mean_word_length': 5.468518300576674, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Microstructural constitutive model for polycrystal viscoplasticity in cold and warm regimes based on continuum dislocation dynamics

Viscoplastic flow of polycrystalline metallic materials is the result of motion and interaction of dislocations, line defects of the crystalline structure. In the microstructural (physics-based) constitutive model presented in this paper, the main underlying microstructural processes influencing viscoplastic deformation and mechanical properties of metals in cold and warm regimes are statistically described by the introduced sets of postulates/axioms for continuum dislocation dynamics (CDD). Three microstructural (internal) state variables (MSVs) are used for statistical quantifications of different types/species of dislocations by the notion of dislocation density. Considering the mobility property of dislocations, they are categorized to mobile and (relatively) immobile dislocations. Mobile dislocations carry the plastic strain (rate), while immobile dislocations contribute to plastic hardening. Moreover, with respect to their arrangement, dislocations are classified to cell and wall dislocations. Cell dislocations are those that exist inside cells/subgrains, and wall dislocations are packed in (and consequently formed) the subgrain walls/boundaries. Therefore, the MSVs incorporated in this model are cell mobile, cell immobile and wall immobile dislocation densities. The evolution of these internal variables is calculated by means of adequate equations that characterize the dislocation processes dominating material behavior during cold and warm monotonic viscoplastic deformation. The constitutive equations are then numerically integrated; and the constitutive parameters are determined/fitted for a widely used ferritic-pearlitic steel (20MnCr5).


Introduction
Nowadays, finite element (FE) simulation of manufacturing processes such as metal forming is an important part of process and product design and development in the industry. Correct and accurate description of material behavior and properties is always the biggest challenge in simulation of industrial manufacturing processes that are based on viscoplastic deformation. Dislocation-density-dependent physics-based constitutive models of metal plasticity while are computationally efficient and history-dependent, can accurately account for varying process parameters such as strain rate and temperature. Since these models are founded on essential phenomena dominating the deformation, they have a wide range of usability and validity. Moreover, they are suitable for manufacturing chain simulations as they can efficiently compute the cumulative effect of the various manufacturing processes by following the microstructure state through the entire manufacturing chain including interpass periods and give a realistic prediction of material behavior and final product properties. The constitutive models are mainly divided into the following main categories (Lin and Chen, 2011;Rusinek et al., 2010):

Empirical constitutive models
Empirical constitutive models provide description of the yield/flow stress based on empirical observations, and consist of some mathematical functions that lack the physical background. In these models, yield stress is usually an explicit function of accumulated plastic strain, strain rate and temperature, which makes empirical models not history dependent. Moreover, the problem with the accumulated plastic strain is that it is a nonmeasurable virtual variable. Furthermore, empirical constitutive models are normally characterized by reduced number of material constants and easy calibration. However, due to their empirical characteristics, they are usually covering limited range of applicability and flexibility and offer low accuracy. Empirical models are determined by fitting parameters of model equations to experimental data without considering the physical processes causing the observed material behavior. Empirical or phenomenological models are also named engineering models as they are more common in engineering applications than the physics-based material models. A number of common empirical constitutive models of metals viscoplasticity are usually incorporated in commercial FE programs (Follansbee and Kocks, 1988;Hockett and Sherby, 1975;Johnson and Cook, 1983;Khan and Liang, 1999;Rusinek and Klepaczko, 2001;Sung et al., 2010).

Microstructural constitutive models
Microstructural or physics-based constitutive models account for microstructural (physical) aspects of the material behavior. These are the models where knowledge about the underlying microstructural processes including dislocation processes, is applied to formulate the thermo-micro-mechanical (TMM) or thermomechanical-microstructural constitutive equations. In addition, since microstructural material models simulate the main microstructural phenomena influencing the overall mechanical response of the material to plastic deformation, they can be used in wide range of deformation parameters (strain rate and temperature) and loading/deformation modes (tension, compression, creep, and relaxation). Additionally, since in industrial metal forming processes, material usually undergoes variety of loading types and parameters, history-dependent microstructural constitutive models are much more suitable and robust for comprehensive simulations of complex industrial metal forming processes. Physics-based models may follow different approaches to describe microstructure evolution/kinetics in polycrystalline metallic materials under plastic deformation:  Discrete dislocation dynamics (DDD) in which slip/glide/motion and interaction of individual dislocations are considered; and thereby, the stress-strain response of the material is a result of direct simulation of a huge assemble of dislocations in a very small representative volume element (RVE). Some good examples can be found in works of van der Giessen and Needleman (1995), Zbib et al. (1998), Devincre et al. (2001), Zbib and La Diaz de Rubia (2002), Arsenlis et al. (2007), Groh et al. (2009), Zhou et al. (2010), Huang et al. (2012) and Chandra et al. (2018). The algorithms based on DDD are extremely costly in terms of computation time, they do not account for size effect, and cannot be readily implemented in standard FE software for industrial application. Nevertheless, DDD simulations are more efficient than those of molecular dynamics (MD) because the RVEs in DDD are much bigger in size than those used in MD which run in atomistic level. DDD simulations provide insights into larger scale behavior (mesoscale). Therefore, models based on MD and DDD are very helpful for studying of dislocation processes and construction of statistical continuum models based on dislocation density (Kubin, 2013;Li et al., 2014;Monavari et al., 2016). Furthermore, DDD and MD have proven to be very useful tools for stochastic modeling of microplasticity experiments such as micro-pillar/column compression, micro-bending and nano-indentation, that together with (in-situ) electron microscopy provide a deep understanding into collective behaviors of dislocations such as dislocation sources, arrangements, configurations and interactions (Csikor et al., 2007;Cui et al., 2017Cui et al., , 2016aCui et al., , 2016bDerlet and Maaß, 2013;El-Awady et al., 2009;Greer et al., 2008;Lee et al., 2016;Miller et al., 2004;Motz et al., 2008;Ng and Ngan, 2008;Oh et al., 2009;Papanikolaou et al., 2018;Parthasarathy et al., 2007;Po et al., 2014;Shiari et al., 2005;Yamakov et al., 2002;Zaiser, 2013;Zhang et al., 2015;Zhu, 2004).  Continuum dislocation dynamics (CDD) describe the microstructure indirectly, so that the effects of the micro level processes are accounted for, in an average way on the macro level. Such type of approach is the subject of this study using the notion of dislocation density which unlike (accumulated) plastic strain, is measurable to some extent by electron microscopy and X-ray techniques. Due to their physical nature, besides plastic/strain/work hardening, constitutive models based on different types of dislocation densities have the potential of predicting many other important processes such as creep, relaxation, dynamic strain aging, static aging, and bake hardening.
Opposite to the DDD approach, constitutive models based on CDD are formulated at the macro level, i.e. the microstructural (internal) state variables (MSVs) are calculated for a mesoscale representative material volume (RMV). In macroscale simulation of plastic deformation, material/integration points are considered to be RMVs. Additionally, an ideal test sample under homogenous uniaxial normal load (tension or compression) which is used for obtaining flow curve is assumed to be a RMV.
With indirect approach of CDD, dislocation density-dependent constitutive models provide a bridge between the micro-level phenomena and macro-level continuum quantities, such as stress and strain rate. Furthermore, simulations performed using these constitutive models are much less costly (in the same range of common empirical constitutive models) and less complicated compared to the algorithms based on discrete dislocation dynamics. Hence, they can be easily implemented in standard FE software and are suitable for industrial applications. This study is limited to the isotropic case meaning that Bauschinger, asymmetrical and anisotropic effects are not included in the presented TMM constitutive model. However, the constitutive relations developed in this context can be applied in the crystal plasticity framework in order to account for anisotropic effects caused by nonuniform dislocation density evolution on each slip system which is negligible in steel alloys. Viscoplastic deformation of crystalline materials with respect to temperature may occur in one of the following regimes/domains:  cold regime: the maximum temperature in cold regime is normally characterized by temperatures above which diffusion controlled dislocation mechanisms such as dislocation climb and pinning become dominant (approximately < 0.3 , where is the absolute temperature; and is the melting absolute temperature) (Galindo-Nava and Rae, 2016);  warm regime: warm viscoplastic deformation of crystalline materials occur above cold but below hot temperature regime (approximately 0.3 < < 0.5 ) (Doherty et al., 1997;Sherby and Burke, 1968); and  hot regime: hot viscoplastic deformations are carried out above warm temperature regime. Hot metal forming processes are characterized by at least one of the hot/extreme microstructural processes such as recrystallization, phase transformation, notable precipitate processes, etc. (roughly 0.5 < < ). Strain rate has different regimes as well, however, independent from the material (Field et al. (2004)):  creep or static: ̇< 10 −4 s −1 (where ̇ is the strain rate);  quasi-static: 10 −4 s −1 ≤̇< 10 −2 s −1 ;  intermediate-rates: 10 −2 s −1 ≤̇≤ 10 s −1 ;  dynamic: 10 s −1 <̇≤ 10 3 s −1 ; and  shock/highly-dynamic: ̇> 10 3 s −1 .
In the present paper, the focus is on monotonic viscoplastic deformation of polycrystalline metallic materials in cold and warm temperature regimes with quasi-static to intermediate strain rate levels. Nevertheless, further studies must be conducted for generalization to hot regime probably in the framework of continuum microdynamics (CMD) which its scope encompasses CDD. Microstructural constitutive models based on CDD and CMD can be coupled with finite element method as microstructural solvers (in addition to the regular thermal and mechanical FE solvers) to simulate not only metal forming processes but also the entire material processing chain including casting, heat treatment, interpass periods, etc., one after the other. However, the microstructural constitutive model based on CMD are characterized with application of extra MSVs in addition to different types of dislocation density. These additional statistical MSVs can be phase volume fractions, void volume fraction, recrystallized volume fraction, twinned volume fraction, transformed volume fraction, precipitate concentrations, etc. In this paper, postulates of CDD are listed. Based on these postulates, a microstructural constitutive model for polycrystal isotropic viscoplasticity in cold and warm regimes is derived. Model's kinetics differential equations are then numerically integrated and subsequently its parameters are determined for a ferritic-pearlitic steel alloy 20MnCr5 which is widely used in industrial forging of automotive components such as bevel gears.

Background
The foundation of CDD was formed in 1930's when the pioneers of the theory, Orowan (1934) and Taylor (1934) introduced the concept of dislocation density and its relationship with plastic strain rate and yield stress. They considered the mean effect rather than individual aspects of dislocations motion and interactions in an attempt to describe macroscopic plastic flow. Johnston and Gilman (1959) were the first to propose an evolutionary equation for dislocation density, which was simply the superposition of a multiplication term and a recovery term. Webster (1966) applied an analogous methodology to creep by assuming that the time rate of change of dislocation density is due to multiplication, immobilization and annihilations processes. Later, widespread adoption of this approach was established by the works of Bergström (1970), Kocks (1976), Mecking and Kocks (1981).
Later, many physics-based constitutive models were proposed with more than one type of dislocation density (multi-MSV models). Ananthakrishna and Sahoo (1981), Bammann and Aifantis (1982), Estrin and Kubin (1986) and Hähner and Zaiser (1999) classified dislocations with respect to their mobility feature. They predicted the flow curve by constitutive models based on two MSVs, mobile/glissile/glide and immobile/sessile dislocation densities. 2-MSV models of mobile and immobile dislocation densities are still being developed and applied in different frameworks (Austin and McDowell, 2011;Hansen et al., 2013;Li et al., 2014). Mughrabi (1983), Nix et al. (1985) and Estrin et al. (1998) differentiated dislocations with regard to their arrangement in the dislocation network; and introduced models with two MSVs, cell and wall dislocation densities. Estrin et al. (1996) and Roters et al. (2000) proposed models with three MSVs, cell mobile, cell immobile, and wall immobile dislocation densities. These models appreciate different dislocations based on their mobility property and their arrangement. Likewise, Blum et al. (2002) approached the creep problem in metals by using a 2-MSV constitutive model that decomposed dislocations with respect to their singularity property, namely singular and dipolar dislocations. There exist models that account for dislocation character (edge and screw dislocations) (Cheong and Busso, 2004) and dislocation polarity (right-hand and left-hand) (Arsenlis and Parks, 2002;Roters, 2011). Ma and Roters (2004) classified dislocations further by allocating extra MSVs to their model for parallel and forest dislocation densities. Moreover, Estrin and Mecking (1992) incorporated effective grain size to the classical single variable Kocks-Mecking model Mecking and Kocks (1981). Sandström and Lagneborg (1975), Busso (1998), Mukherjee et al. (2010) and Babu and Lindgren (2013) all developed multi-MSV microstructural constitutive models to characterize metals behavior under hot deformation by introduction of an additional statistical MSV, recrystallized fraction, with its corresponding kinetics equation. Fan and Yang (2011) and Bok et al. (2014) went further by allocating additional MSVs for each phase fraction to build a microstructural constitutive model for hot sheet metal forming.
Recently, some microstructural constitutive models have been developed for special purposes. Viatkina et al. (2007), Kitayama et al. (2013), Pham et al. (2013), Knezevic et al. (2013) and Zecevic and Knezevic (2015) proposed multi-MSV models for strain path dependent evolution of dislocation structures during cyclic plastic deformation to account for kinematic hardening and Bauschinger effect. Austin andMcDowell (2011), Lloyd et al. (2014), Luscher et al. (2017) and Nguyen et al. (2017a) utilized dislocation density-based constitutive modeling for viscoplastic deformation of metals at dynamic and shock regimes. Patra and McDowell (2012) developed a physical-based constitutive model for inelastic deformation of irradiated bcc ferritic-martensitic steels by introduction of an extra statistical MSV, namely number of interstitial loops that are formed due to irradiation. Bouaziz and Guelton (2001), Allain et al. (2004) and Steinmetz et al. (2013) incorporated twinned volume fraction as an additional MSV to existing dislocation density-based models in order to reflect plastic hardening behavior of twinning-induced plasticity (TWIP) steels. These models were further developed and implemented in crystal plasticity framework by Wong et al. (2016) to account also for transformation-induced plasticity (TRIP) effect by incorporating one more MSV, namely transformed volume fraction. In addition, Kubin et al. (2002) and Ananthakrishna (2007) have published reviews on theoretical approaches for modeling of collective behavior of dislocations which readers are encouraged to read.
In the microstructural state variable approach, the mechanical state at a nonlocal/macroscale material point in a continuum body is characterized in terms of hidden or internal variables that statistically represent the nonlocal microstructural state (NMS) in addition to the observable statistical/nonlocal external variables such as temperature, strain rate and yield stress. To date, many statistical physics-based approaches have been applied for constitutive modeling of metals. Even though extensive work has been conducted in this area, there is not a universal agreement on the number and kind of MSVs to be used (Horstemeyer and Bammann, 2010) as well as the influencing microstructural processes and their associated kinetics equations that determine the values of MSVs. This study is an attempt for unification and completion of the previous works in the field of continuum dislocation dynamics, by gathering and modifying some of the most important postulates of the CDD theory. In the following section, it is argued that the statistical state of microstructure of pollycrystalline materials under monotonic and isotropic viscoplastic deformation in cold and warm regimes is fully defined by three MSVs that are cell mobile, cell immobile and wall immobile dislocation densities. For the first time, evolution of these MSVs with respect to time (or plastic strain) are described considering every statistically notable dislocation process/interaction affecting values of the aforementioned dislocation densities. Without comprehensive decomposition of dislocation types and processes as suggested in this paper, accurately capturing the mechanical response of complex metal alloys such as steels particularly in warm regime is not achievable.

Postulates and the constitutive model
To construct a constitutive model based on microstructural processes, first a set of axioms or postulates must be established as the basis for reasoning and subsequent derivation of constitutive relations. In the framework of CDD, the following postulates are introduced, although not all of them are independent. Additionally, for derivation of the constitutive equations, the consequences of each postulate in combination with the earlier ones (or the results of earlier postulates) are provided as well.

Fundamental postulates
1) Nonlocal microstructural state (NMS): the mean/nonlocal yield/flow/critical shear stress resolved at slip systems (̃) as the nonlocal mechanical response of material is an implicit function of temperature, shear strain rate, and statistical state of microstructure, given by the following equation (Mecking and Kocks, 1981): where ̃̌ is the TMM constitutive function; ̃̇ is mean plastic shear strain rate at slip systems; is temperature, and is a set containing all MSVs and is referred to as nonlocal microstructural state since it represents the statistical state of microstructure: where is the i-th nonlocal MSV which can be various types of dislocation density, grain size, phase fractions, recrystallized fraction, precipitate concentration and size, etc. The evolution of each MSV which is often expressed as time rate of change of MSV, is a function of thermo-mechanical loading and (Nadgornyi, 1988): ̇≡̇̌( , ,̃̇) =̇̌( , 1 , 2 , 3 , … ,̃̇) ; where ̇̌ is the function that determines evolution of with respect to time. As shown in Fig. 1, nonlocally represents the microstructural state of a material point on a macroscale continuum body. As deformation proceeds, the microstructure state evolves towards a saturation/steady state sat (Estrin and Mecking, 1984). Furthermore, it is emphasized that in this description, the NMS set and its components (MSVs) are treated based on the nonlocal principle/treatment (Eringen, 1983;Gao and Huang, 2001) as they are averaged over a mesoscale RMV. The mesoscale RMV must be a polycrystalline aggregate that represents the bulk material properly at the considered nonlocal continuum material point. Hence, it must consist of sufficient number of constituent single crystal grains in order to capture the size effect (Adams and Olson, 1998). This implies that inside the mesoscale RMV, locally the values of MSVs are not necessarily equal to their integral average over RMV due to their heterogeneous (local) distribution. For instance, some MSVs are highly concentrated at small regions while some are statistically distributed inside the grains of RMV (postulate (3)). It should be noted that accumulated plastic strain is not an MSV but a mechanical internal state variable and due to its virtual characteristic, it cannot be measured directly. Two identical material samples that are plastically (uniaxially) deformed to an equal amount of accumulated plastic strain but with different histories of temperature and strain rate (e.g. in cold or warm regime), if again deformed under an equal temperature and strain rate condition, do not necessarily yield the same stress response.
Furthermore, the mean/nonlocal normal yield/flow stress ( ) (or simply the yield stress) nonlocally applied on a polycrystalline aggregate is related to the mean yield shear stress resolved at its slip systems (̃) by Taylor factor (Kocks, 1970;Taylor, 1938): where ̇ is the mean/nonlocal normal plastic strain rate (or simply the plastic strain rate). It is safe to assume constant ≅ 3 for a random orientation distribution of texture for bcc and fcc polycrystalline aggregates. Nevertheless, Taylor factor evolves as plastic strain accumulate; and it is also dependent on the deformation mode but for sufficiently random-textured polycrystals (weak texture or textureless) with random loading during deformation (isotropic case), these dependencies can be neglected (Kocks and Mecking, 2003).
2) Dislocation mobility: with respect to their mobility property, dislocations are divided into mobile and immobile dislocations. While mobile dislocations carry plastic strain, immobile dislocations contribute to plastic hardening (Estrin and Kubin, 1986;Hunter and Preston, 2015). Upon confronting obstacles, mobile dislocation segments may become fully or partially immobilized. Partial immobilization of a (prior) mobile dislocation segment and consequently its division to mobile and immobile dislocation segments is schematically shown in Fig. 2. As illustrated in Fig. 2, an individual dislocation might be consisted of several mobile and immobile dislocation segments. During movement, the length and mean bow-out radius of a bowed-out mobile segment bounded by adjacent immobile segments is reduced while the lengths of its bounding immobile dislocation segments are increased proportionally, as long as the immobile segments are not remobilized by remobilization mechanisms. As pointed out in postulate (1), according to the non-locality principle, mobile and immobile dislocation densities at each nonlocal material point are defined as follows: where is the considered volume which in nonlocal case, is the volume of mesoscale RMV; ( ) and ( ) are respectively lengths of j-th mobile dislocation segment and immobile dislocation segment of type x; d ( ) and d ( ) are infinitesimal elements of arc length along the dislocation segments ( ) and ( ) , respectively; and and are (nonlocal) mobile dislocation density and immobile dislocation densities (of type x), respectively. x can be c or w that respectively stand for cell and wall (see postulate (3)).

Fig. 2.
Schematic of a mobile dislocation segment that is partially immobilized by immobile dislocations intersecting its slip plane (adopted from Hunter and Preston (2015)).
As shown in Fig. 2, it is postulated that motion of bowed out mobile dislocation segments moving on slip systems is always restricted at their both ends by immobile dislocation segments (either cell or wall immobile dislocations). Moving a bowed-out mobile dislocation segment requires a certain amount of shear stress acting on it which is inversely related to the mean radius of bowed-out segment (Hull and Bacon, 2011). In postulate (8), a statistical relation between bow-out radius and immobile dislocation density is derived. Hence, only the collection of immobile dislocations determines the required mean shear stress acting on bowed-out mobile dislocation segments to move them in order to accommodate plastic straining.
3) Dislocation arrangement: due to heterogeneous distribution of dislocations in crystal grains, dislocations with regard to their arrangement are categorized to cell dislocations that exist inside the cell/subgrain blocks/interiors and wall dislocations that form the cell walls/boundaries (Mughrabi, 1983). Wall dislocations are immobile and closely clustered (dense dislocation walls); and statistically, their Burgers vector is very similar to their surrounding wall dislocations. As such, they cause considerable lattice torsion/warp/twist/bending that is manifested by measurable misorientation angle across wall regions that can be detected and measured by electron backscatter diffraction (EBSD) (Gardner et al., 2010;Jiang et al., 2013) and X-ray methods. These dislocations are often referred to as geometrically necessary dislocations (GNDs) (Ashby, 1970;Nye, 1953) as they accommodate lattice curvature (incompatibility of plastic deformation) that arises by geometry change due to local gradient of plastic strain (Arsenlis and Parks, 1999;Gao and Huang, 2003) through formation of semi-planar geometrically necessary boundaries (GNBs) or dense dislocation walls (DDWs) (Kuhlmann-Wilsdorf and Hansen, 1991). Consequently, GNDs are in fact wall dislocations that contribute the most to plastic hardening due to the long range internal stress produced by them (Kassner et al., 2013;Mughrabi, 2006). As plastic straining proceeds, the degree of misorientation angle between adjacent cells or cell blocks (CBs) increases. Physically, GNBs containing a high local dislocation density with a net Burgers vector are very different than spatially relatively random distributions of cell dislocations (Hughes et al., 2003). Cell dislocations that can be either mobile or immobile do not necessarily adopt any considerable particular semi-stable arrangement unless they become part of walls. Thus, cell dislocations are assumed to be statistically/homogenously distributed inside the CB structure (subgrain); and hence are known as statistically stored dislocations (SSD) (Ashby, 1970). However, stationary cell immobile dislocations form another type of semi-temporary accumulates/clusters/ pile-ups/bundles/tangles/nets named incidental dislocation boundaries (IDBs) (Kuhlmann-Wilsdorf and Hansen, 1991), or forest dislocations with relatively negligible misorientation angle. IDBs form secondary dislocation cells or sub-cells (SCs) inside subgrains. Moreover, the cell volumes bounded by GNBs may form tiny channel-shaped shear/deformation bands such as micro-bands (MBs) and lamellar bands (LBs) (Bay et al., 1992;Bay et al., 1989;Hughes, 1993;Hughes and Hansen, 1993). Fig. 3 schematically illustrates cell-wall substructure inside crystal grains. Therefore, there are three main independent types of dislocation density, cell mobile dislocation density ( ), cell immobile dislocation density ( ) and wall immobile dislocation density ( ): where denotes dislocation density; subscripts , , , , and respectively stand for total, immobile, mobile, cell and wall; and SS and GN represent statistically stored and geometrically necessary, respectively. To sum up, wall immobile dislocations are locally (inside the RMV) highly concentrated/compacted at GNBs. Thus, local density of wall immobile dislocations at GNBs is much higher than density of cell immobile dislocations at IDBs or sub-cell interiors. Mobile dislocations are nearly homogeneously distributed in subcells. Consequently, local density of (cell) mobile dislocations is much lower than immobile dislocations.
Initial cell immobile dislocation density is determined by density of trapped and accumulated cell dislocations. Wall immobile dislocation density of an undeformed polycrystal is determined by the area of grain boundaries, phase and precipitate interfaces and other two-dimensional crystal defects. Even though grain boundaries and interfaces do not fit in classical definition of dislocations, as they are two-dimensional crystal defects like dislocations and because they interact with dislocations (acting as relatively impenetrable obstacles), technically, they can be included in undeformed wall immobile dislocation density. The equivalent dislocation density of these interfaces is relatively small enough that does not lead to a considerable error in the calculated initial yield stress, although they largely influence the initial plastic hardening. On the other hand, their equivalent density is not affected by plastic straining in cold and warm regime (constant during deformation). Since the interfaces are stable and strictly immobile in cold and warm regime, at the earliest stage of deformation they quickly transform to walls due to the misorientation (relative to undeformed configuration) they inflict, and hence are considered as sources of walls. Therefore, in the undeformed material state, (initial) wall immobile dislocation density depends on the effective grain size which includes the influence of all the interfaces (size effect). Furthermore, in some dual phase metals (e.g. ferritic-martensitic DP steels), due to phase transformation associated expansion during quenching, softer phase (ferrite) becomes plastically deformed by the harder phase (martensite) which results in formation of wall immobile dislocations (GNDs) in the softer phase. In such cases, initial wall immobile dislocation density has a relatively high value. 4) Viscoplastic decomposition: the mean yield shear stress at slip systems (̃) has two major contributions that obeys the linear superposition rule (Kumar et al., 1968;Mecking and Kocks, 1981): where ̃ is plastic/athermal/rate-independent/internal/back shear stress; and ̃ is referred to as viscous/rate-dependent/thermal/effective shear stress/drag, mean Peierls-Nabarro stress, or overstress which is the mean viscous lattice resistance to move mobile dislocations in a nearly obstacle/dislocation-free lattice (with relatively very low dislocation density) (Nabarro, 1997(Nabarro, , 1952Peierls, 1940). Viscous shear stress is affected by point defects such as vacancies, alloying elements and solute atoms (interstitial and substitutional). Additionally, plastic shear stress needs to be overcome to move the bowed-out mobile dislocations (Fig. 2). As mentioned in postulate (2), the lower the bow-out radius, the lower the required (extra) plastic shear stress to move it. The viscoplastic decomposition can be depicted by the rheological model shown in Fig. 4, similar to Perzyna-type formulation (Perzyna, 1966), which consists of a parallel set of nonlinear dashpot/damper and nonlinear friction elements that are in series linkage with a linear spring element. Furthermore, given Eq. (4), Eq. (11) may take the following form: 5) Cell-wall decomposition: as emphasized in postulate (2), only immobile dislocations contribute to plastic hardening. The total surface density (total surface area per unit volume) of cells and walls are equal, as dislocation walls (GNBs) encompass the cell volumes (subgrains). In other words, the cell-wall dislocation substructure can be treated as a cellular composite. Hence, it is postulated that the two types of immobile dislocations contribute to plastic hardening in parallel (additive decomposition of plastic stress) in order to satisfy the compatibility for isochoric viscoplastic flow of material in which both cells and walls need to be deformed simultaneously at an equal rate. Thus, mean plastic shear stress is linearly decomposed to its constituent internal stresses corresponding to cell immobile dislocations (immobile SSDs) and wall immobile dislocations (immobile GNDs) (Columbus and Grujicic, 2002;Mughrabi, 1987;Voyiadjis and Al-Rub, 2005): where ̃ and ̃ are respectively contributions of cell and wall immobile dislocations to plastic shear stresses. ̃ and ̃ are axial and torsional plastic stress, respectively. In postulate (8), the relations of ̃ and ̃ with immobile dislocation densities are described. Considering Eq. (4), Eq. (13) is rewritten as follows: 6) Local dislocation density: in saturation state, dislocation substructure and density is steady and statistically remain constant in a representative material volume. One can assume that the exact locations of walls (GNBs) inside RMV at saturation state are known prior to viscoplastic deformation. Hence, is assumed to be the summation of immediate (local) surrounding volumes of GNBs in the saturation state (including mobile dislocations constrained by GNBs constituent wall immobile dislocations). Cell immobile dislocations and the mobile dislocations bounded by them are envisaged to be uniformly distributed in the rest of RMV's volume ( ). Therefore, Notice that , and are virtual volumes. Moreover, (volume of RMV) is constant during deformation. In monotonic viscoplastic deformation, it is known that GNBs are relatively sharp boundaries containing locally dense wall immobile dislocations (postulate (3)). This indicates that ≪ and ≪ . However, in materials under long and high-amplitude cyclic viscoplastic deformation, and become closer but still often < . Thus, in monotonic viscoplastic deformation which is the subject of this study, one can assume that cell and wall immobile dislocations during the entire deformation are homogenously locally distributed over constant volumes of and , respectively. Therefore, According to Mughrabi (1983), constant cell volume fraction ( ) and wall volume fraction ( ) are defined as follows: Combining Eqs (6), (16) and (17) results in (Mughrabi, 1983): From the view point of an observer who observes the macroscale material point shown in Fig. 1, viscoplastic deformation is homogenous. This leads to the assumption of nearly uniform local distribution of mobile dislocations which carry the plastic strain. Consequently, local and nonlocal mobile dislocation densities are almost equal ( loc = ). Combining this with Eq. (18) and the fact that is constant gives: 7) Dislocations characteristic spacing: dislocations intersect each other and form dislocation network with characteristic average spacing of Λ. It can be easily shown that the average spacing among dislocations (also known as mean free path of dislocations) of type (Λ ) is inversely proportional to the square root of respective local dislocation density (Kocks, 1966;Nes, 1997;Seeger, 1955): It should be noted that Λ has a local character. Additionally, it can be also readily shown that the volumetric number of junctions/intersection/nodes of dislocation type or dislocation junction density of type ( ) is proportional to dislocation density loc or (Eq. (20)) while inversely proportional to the average spacing Λ (Gottstein and Argon, 1987): 8) Plastic stress: as pointed out in postulate (2), the average bow-out radius of cell mobile dislocations that are constrained at both ends by immobile dislocations of type = , , (̃) (Fig. 2) is proportional to average spacing of immobile dislocations of type (Nes, 1997): On the other hand, the mean plastic resolved shear stress at slip systems is directly proportional to Burgers length (magnitude of Burgers vector) ( ) and shear modulus ( ) and inversely proportional to ̃ Hull and Bacon, 2011;Nabarro, 1952): where ̃ is the mean plastic resolved shear stress associated with immobile dislocations of type . Given Eqs. (21) and (23), the Taylor relation (Bailey and Hirsch, 1960;Seeger et al., 1957;Taylor, 1934) is derived: where loc is a material constant known as local dislocation interaction strength/coefficient; and ̃l oc is statistical average of loc for different configurations of interacting mobile-immobile dislocations in various slip systems at the considered local point. Given Eq. (18): where ̃ is the nonlocal interaction strength related to local density and geometrical arrangement of immobile dislocations of cell and wall species ( = , ). Variation of ̃ with plastic strain is assumed to be negligible (Kocks and Mecking, 2003). Since most of the plastic hardening is due to wall immobile dislocations, essentially ̃>̃. Further, combining Eqs. (19) and (26) yields: Finally, given Eq. (4), Eq. (26) may be rewritten as: 9) Kinetics superposition: dynamic evolution of dislocation density is described as the rate of change of dislocation density with respect to time or plastic strain which is linear superposition of increase/production (with positive sign) and decrease/elimination (negative sign) terms (Johnston and Gilman, 1959;Kocks, 1976;Webster, 1966): Therefore, where ̃≡̃ and ≡ are partial derivative operator with respect to mean shear plastic strain at slip systems and mean plastic strain, respectively. Static and overall (combination of static and dynamic) evolution of dislocation densities are represented by time derivative relations: where subscripts and respectively represent static and dynamic states. Furthermore, for convenience in calculations and dimensional balancing of equations, the normalized/dimensionless dislocation density of type (̂) is defined as follows: where the hat-sign (̂) indicates the normalization; and 0 is constant reference dislocation density. Therefore, in all the previous relations that are presented so far in this paper (except Eq. (26) and (28)) can be replaced by ̂. Thereby, Λ and ̂ are respectively normalized average spacing among dislocations of type and normalized dislocation junction density of type .
10) Dynamic dislocation processes: there are six main classes of dynamic dislocation processes that are statistically considerable in constitutive modeling of metal viscoplasticity in cold and warm regimes:  Generation/multiplication of mobile dislocations: Motion/displacement of mobile dislocations leads to their elongation.  Annihilation of dislocations: dislocations of all three types, upon contact, can be annihilated by mobile dislocations with opposite Burgers vector slipping on neighboring parallel slip planes.  Accumulation of immobile dislocations: mobile dislocations can be immobilized by immobile dislocations and produce accumulates of immobile dislocations.  Trapping of mobile dislocations: mobile dislocations interacting each other can be immobilized by trapping process. Trapping process has two underlying mechanisms/reactions that result in immobilization of infected parts of mobile dislocations: locking of mobile dislocations; and pinning of mobile dislocations by interstitial solute/impurity atoms. Dynamic strain aging (DSA) effect in part is associated with pinning of mobile dislocations. At macroscopic level, pinning process reveals itself in the existence of enhanced upper initial yield stress which is followed by yield point elongation associated with Lüders bands (Hahn, 1962;Hall, 1970) in stress-strain curves. Pinning phenomena and their macroscopic effects are classified in Fig. 5.  Nucleation of wall dislocations: cell immobile dislocations can become wall dislocations if the accumulates they belong to, while are being immobilized sufficiently strong, reach a critical size and local density.  Remobilization/mobilization of immobile dislocations: immobile dislocations can be mobilized and contribute to plastic straining. 11) Static dislocation processes: there are two major static (at nonlocal ̃̇≈ 0) dislocation processes during viscoplastic deformation of metals in cold and warm regimes (Pham et al., 2015):  Static pinning (Fig. 5) of cell dislocations by diffusion of interstitial solute atoms such as Carbon and Nitrogen to surroundings of dislocation cores and junctions (Cottrell cloud) which are statistically distributed (SSDs). This is due to higher stored elastic energy density of dislocation surroundings compared to defect-free regions of lattice (Cottrell and Bilby, 1949).  Static remobilization of immobile dislocations by means of thermal mechanisms such as junction dissociation and dislocation climb (vacancy diffusion).
It should be noted that the aforementioned dynamic and static dislocation processes often consist of multiple individual underlying mechanisms. In the following sections, evolution of different types of dislocations by the above-listed dislocation processes is quantitatively and statistically described through several postulates.

Dynamic dislocation processes
12) Dynamic generation: motion of bowing mobile dislocation segments that are held at both ends by immobile dislocations generates more mobile dislocations by different mechanisms (Stricker et al., 2018) such Frank-Read mechanism (Frank and Read, 1950), multiple cross-slip mechanism (Koehler, 1952), or mobile junctions mechanism (Stricker and Weygand, 2015;Weygand, 2014). Thus, the rate of dynamic generation rate of mobile dislocations in plastic shear increment (̃g n ) is proportional to cell mobile dislocation density. The generation rate is higher where average bow-out radius of mobile dislocation segments bounded by immobile dislocations (̃) and the average length of mobile dislocation segments (Λ ∝̃) is larger (Nadgornyi, 1988;Steif and Clifton, 1979): Given Eqs. (4), (21), (23) and (34) where ̂g n is normalized dynamic generation rate of cell mobile dislocation density with respect to plastic strain; gn is probability amplitude of dynamic generation of cell mobile dislocations; and gn is material coefficient associated with probability amplitude of dynamic generation of cell mobile dislocations. Frequency of occurrence of a dynamic dislocation process with probability amplitude of is equal to = ̃̇. In addition, dynamic generation of mobile dislocations is an athermal dislocation process meaning that it is not directly dependent on the temperature.
13) Dynamic annihilation: the annihilation process takes place between a mobile dislocation and another dislocation that can be of any type. Thus, there are three sorts of dislocation annihilation: mutual annihilation of cell mobile dislocations, mutual annihilation of cell mobile and cell immobile dislocations, and mutual annihilation of cell mobile and wall immobile dislocations. Hence, the magnitude of dynamic annihilation rate of mobile dislocations and dislocations of type with respect to plastic shear strain increment (̃) is proportional to cell mobile dislocation density and the dislocation density of type (Ananthakrishna and Sahoo, 1981): Given Eqs. (4) and (34): ̃a n = an ; an = an ̂ ; = , , ; ̂a n = an ̂ ̂ ; = , , ; where ̂a n is normalized dynamic annihilation rate of dislocation density of type with respect to plastic strain; an is probability amplitude of dynamic annihilation of dislocation density of type ; And an are material coefficients related to frequency of dynamic annihilation of dislocation density of type . Annihilation events take place between a mobile dislocation and another dislocation with opposite/antiparallel Burgers vector on parallel planes when two dislocations are within a critical distance from each other. Two near-screw dislocations of opposite sign, slipping on two neighboring parallel slip planes can annihilate by cross slip of one of them (Brown, 2002;Essmann and Mughrabi, 1979;Nix et al., 1985;Oren et al., 2017;Pauš et al., 2013;Püschl, 2002;Seeger et al., 1957). Two near-edge dislocations may annihilate as well by spontaneous formation and disintegration of a very narrow unstable dislocation dipole which can be assisted by dislocation climb. (Eisenlohr and Blum, 2005;Monavari and Zaiser, 2018;Vegge and Jacobsen, 2002). Thus, thermally activated dislocation climb and cross-slip mechanisms facilitate dislocation annihilation. Therefore, dynamic annihilation of dislocations are thermal dislocation processes meaning that probability amplitudes associated with dislocation annihilation processes have temperature dependence (increase with increasing temperature). On the other hand, since at least dislocation climb as a contributing mechanism for dynamic annihilation processes is controlled by time-dependent diffusion of vacancies, annihilation processes have strain rate dependence with negative strain rate sensitivity. Besides, dislocation cross-slip can be as well considered a rate-controlled mechanism (Nes, 1997).
Thus, one might suggest that when a mobile dislocation and another dislocation (of type ) are on the same slip system within the critical distance (some authors call it capture radius) for annihilation of dislocations of type ( an ), they will likely annihilate each other, in case of having opposite Burgers vectors (Essmann and Mughrabi, 1979). In light of this, for instance, the dimensionless material parameter associated with frequency of dynamic annihilation among cell mobile dislocations ( an ) can be decomposed as an = añ , where ̃ is the average number of active slip systems (̃≥ 3). Statistically, one can assume equal density of dislocations on all active slip systems (in isotropic case) which gives rise to normalization factor 1 (Roters et al., 2000).
14) Dynamic accumulation: immobilization of mobile dislocations by immobile dislocations produces or increases the size of immobile dislocation accumulates. This dislocation process is known as dislocation accumulation; and its magnitude with respect to plastic shear strain increment (̃a c ) is proportional to cell mobile dislocation density. It is also inversely proportional to average spacing among immobile dislocations of type = , (Λ ) (Kocks, 1976) as well as the average radius of mobile dislocation segments bounded by immobile dislocations of type = , (̃∝ Λ ): Having Eqs. (4), (21) and (34): ̃a c = ac ; ac = ac √̂ ; = , ; ̂a c = ac √̂ ̂ ; = , ; where ̂a c is normalized dynamic accumulation rate of immobile dislocations of type with respect to plastic strain; ac is probability amplitude of dynamic accumulation of immobile dislocations of type ; and ac are material parameters associated with probability of dynamic accumulation of immobile dislocations of type . Dynamic accumulation of immobile dislocations is an athermal process. 15) Dynamic trapping: interaction between mobile dislocations moving in different intersecting slip systems results in the formation of dislocation junctions/intersections/nodes. These dislocation junctions can be lock junctions such as Lomer-Cottrell lock junctions (Cottrell, 1953;Lomer, 1951), Hirth lock junctions (Hirth, 1961), collinear (Madec et al., 2003) and coplanar lock junctions (Thompson, 1953). Lock junctions restrict movement of involved mobile dislocations by immobilizing parts of them near the formed lock junction, leading to latent hardening (Franciosi, 1985). In addition, some junctions formed due to interaction of mobile dislocation pairs are not lock junctions (depending on the configuration of respective slip systems and mobile dislocations). They can as well be mobile/glissile or temporal/temporary junctions which do not lead to immediate immobilization directly. Essentially, dislocation junctions due to their relative high energy density are attractive regions for diffusion and accumulation of interstitial solute atoms. Therefore, some of the mobile and temporal junctions become pinned/arrested/anchored due to diffusion and accumulation of interstitial solutes at their vicinity; which leads to subsequent immobilization of respective mobile dislocations or infected parts of them (Cottrell and Bilby, 1949;Mulford and Kocks, 1979;van den Beukel and Kocks, 1982). Density of potential junctions of interacting mobile dislocations, as mentioned in postulate (7), is proportional to density of mobile dislocations ( ) and inversely proportional to the average spacing among mobile dislocations (Λ ). Thus, frequency of occurrence of mobile dislocation trapping events comprised of two major underlying mechanisms (locking and pinning), must be proportional to junction density of cell mobile dislocations ( ). Henceforth, ̃t r ∝ .
Given Eqs. (4), (22) and (34) where superscript tr denotes the trapping process; ̂t r is normalized dynamic trapping rate of cell mobile dislocations with respect to plastic strain; tr is probability amplitude of dynamic trapping of cell mobile dislocations; and tr is a material coefficient associated with probability amplitude of dynamic trapping of cell mobile dislocations. Since the pinning mechanism of trapping process is a diffusion controlled reaction, dynamic trapping of cell mobile dislocations is a thermal dislocation process that its magnitude increases with increasing temperature. In other words, by increasing temperature, particularly in metals with notable interstitial concentration, the rate of immobilization of dislocation junctions (trapping) becomes higher. Additionally, since pinning occurs due to diffusion of interstitial solute atoms to dislocation cores and junction, tr is a function of concentration of interstitial solutes atoms such as Carbon and Nitrogen. Moreover, one might argue that tr has a negative strain rate sensitivity because the pinning process is associated with time-dependent diffusion of interstitial solute atoms. The DSA effect is observed when the interstitial atoms reorientation time is shorter than the waiting time of temporal and mobile junctions at the locations nearby interstitial solute complex (Kubin et al., 1988;Kubin and Estrin, 1990). Hence, there exist a critical strain rate beyond which the intensity of mobile dislocation pinning abruptly drops. That strain rate belongs to the dynamic regime. 16) Dynamic wall nucleation: accumulates of cell immobile dislocations that are immobile sufficiently strong (meaning that their remobilization requires higher stress than mean yield resolved shear stress) will grow and become denser by immobilizing more mobile dislocations at their vicinity (cell immobile accumulation process). These cell immobile pile-ups that reach a critical size and local density, subsequently become adequately stable to produce sufficient local stress concentration which results in necessary lattice curvature in their vicinity to make them part of walls. This process is called wall nucleation. From this point of view, wall nucleation process converts some cell immobile to wall immobile dislocations. As such, its rate is proportional to cell immobile dislocation density. On the other hand, the accumulation of cell immobile dislocations increases the probability of wall nucleation process. Hence, statistical magnitude of dynamic nucleation rate of wall immobile dislocations with respect to plastic shear strain increment (̃n c ) is proportional to cell immobile dislocation density ( ) and the rate of cell immobile dislocation accumulation (̃a c ): Given Eqs. (4), (21), (34) and (41): ̃n c = nc ; nc = nc √̂ ̂ ; where ̂n c is normalized dynamic nucleation rate of wall immobile dislocations with respect to plastic strain; nc is probability amplitude of dynamic nucleation of wall immobile dislocations; and nc is a material coefficient associated with probability amplitude of dynamic nucleation of wall immobile dislocations. Pinning of cell immobile dislocations at their junctions by solute atoms contributes to stronger immobilization of cell immobile accumulates and consequently enhances the stability of the respective IDBs which are about to transform to GNBs. This is also suggested by existence of the factor ̂3 /2 in Eq. (49) which is proportional to junction density of cell immobile dislocations ( ∝ 3/2 ), considering that the pinning process occurs mainly at dislocation junctions. From another perspective, given Eq. (26), presence of the factor √̂ in Eq. (48) is proportional to the average plastic shear stress concentration at IDBs. With this approach, twin formation can be considered as a wall nucleation mechanism, in case of slipdominated metals. Regular dislocation walls are formed by high local stress concentration that causes local elastic bending of crystal in metals with relatively high stacking fault energy (SFE). Analogously, twins are also boundaries of lattice misorientation (GNBs) that plastically nucleate due to local stress concentration at sufficiently stable dislocation pile-ups (IDBs) (Christian and Mahajan, 1995;Kibey et al., 2007;Venables, 1964) in metals with adequately low SFE. Each stacking fault that nucleates a twin is bounded by partial dislocations such as Shockley and Frank partial dislocations (Hirth and Lothe, 1982;Hull and Bacon, 2011) that are geometrically necessary (Mahajan and Chin, 1973) and hence are treated as wall dislocations. Therefore, nc also reflects frequency of twin nucleation which in turn depends on the mean SFE at the nonlocal material point under consideration. Twin forming dislocations must overcome the intrinsic SFE barrier to nucleate twin's associated stacking fault. Further, most of the energy required for twinning is spent on creating its associated stacking fault. Thus, especially in relatively low SFE metals such as high manganese TWIP steels, nc has a relatively high value. It is emphasized that pinning of cell immobile dislocations facilitates twin nucleation by increasing stress concentration in the cells to reach the critical shear stress for dislocation separation/dissociation/splitting mechanism (Byun, 2003;Koyama et al., 2015). In slip dominated plastic deformation in which twinning occurs, the share of TWIP effect in total accumulated plastic strain is negligible compared to that of slip. However, as already highlighted, twinning has a remarkable influence on plastic hardening by facilitating the formation of dislocation walls.
As mentioned before, pinning is a diffusion-controlled mechanism that its magnitude increases by increasing temperature in metals with adequate interstitial solute content. Moreover, the rate of stacking fault formation due to dislocation separation mechanism depends on SFE which in turn is generally an increasing function of temperature. In view of this, dynamic nucleation of wall immobile dislocations is a thermal dislocation process that its statistical magnitude by increasing temperature might increase or decrease depending on the metal structure and composition.
Further, since in nucleation of wall immobile dislocations, time-dependent diffusion-controlled pinning mechanism is partly involved, one might argue that nc has negative strain rate sensitivity. On the other hand, wall and twin nucleation can be amplified by increasing strain rate provoking positive strain rate sensitivity of nc . Therefore, overall strain rate sensitivity of nc can be negative, zero or positive (Yang et al., 2017). 17) Dynamic remobilization: probability of recurrence of remobilization of cell or wall immobile dislocations depends on their density (Bergström, 1970): where ̃r m is the remobilization rate of immobile dislocations of type with respect to plastic shear strain increment. Given Eqs. (4) and (34): ̃r m = rm ; rm = rm ; = , ; ̂r m = rm ̂ ; = , .
where ̂r m is normalized dynamic remobilization rate of immobile dislocations of type with respect to plastic strain; rm is probability amplitude of dynamic remobilization of immobile dislocations of type ; and rm are material parameters associated with probability amplitude of dynamic remobilization of immobile dislocations of type . Dislocation remobilization process is performed by different underlying mechanisms such as climb, cross-slip, bow-out, junction dissociation/unzipping, etc. (Hunter and Preston, 2015). Owing to the involvement of thermal mechanisms in the remobilization processes, such as dislocation climb and crossslip, they are treated as thermal dislocation processes. At low and medium temperatures, dynamic dislocation remobilization and annihilation processes (together known as dynamic recovery) are mainly governed by the cross-slip mechanism while at high temperatures dislocation climb is the dominant mechanism of dynamic recovery processes (Essmann and Mughrabi, 1979;Galindo-Nava et al., 2012;Kubin et al., 1992;Nix et al., 1985;Püschl, 2002;Rivera-Díaz-del-Castillo and Huang, 2012). Essentially, by definition inherently ̂r m should be much lower than ̂r m (̂r m ≪̂r m ).
Finally, since dislocation climb as an underlying mechanism for dynamic remobilization processes is controlled by time-dependent diffusion of vacancies, remobilization processes might have strain rate dependence with negative strain rate sensitivity. On the other hand, by increasing strain rates, the viscous stress is also increased which assists the remobilization process. This provokes positive strain rate sensitivity of remobilization parameters. Therefore, rm might have negative, zero or positive strain rate sensitivities.
By approaches analogous to the one adopted in postulate (13) (dynamic annihilation), one may also physically interpret the constitutive parameters associated with frequency of other dislocation processes ( ) by defining various mean critical distances, local densities, shear stresses, etc. Nevertheless, the mean critical physical parameters have virtual characters (like probability amplitude), as they are extremely difficult and mostly even impossible to be determined accurately enough by means of independent experimental measurements. Now that every main dynamic dislocation process is statistically and quantitatively characterized, the approach first introduced by Ananthakrishna and Sahoo (1981), and Estrin and Kubin (1986) is applied to construct the overall kinetics equations for dynamic evolution of each dislocation type.
The term ̂a n is considered twice in dynamic decrease of cell mobile dislocations because two mobile dislocations annihilate each other in the process of dynamic annihilation of cell mobile dislocations.
Eq. (64) is the dynamic evolutionary equation of total dislocation density.

21)
Static pinning: nonlocal static pinning of cell (mobile and immobile) dislocations is very similar to dynamic pinning (postulate (15) and (16)) with the difference that it occurs at very low local strain rates in relatively long durations. Time rate of static pinning of cell dislocations (̇p n ) by diffusion of solute atoms at their junctions is proportional to their junction density ( ): Given Eqs. (22) and (34): where pn is a material coefficient associated with frequency of static pinning of cell dislocations of type = , . Since pinning is a diffusion-controlled process, static pinning of cell dislocations is a thermal dislocation process that its magnitude increases with increasing temperature. As mentioned in postulate (15), since pinning process occurs due to diffusion of interstitial solute atoms to dislocation surrounding, pn is a function of interstitial solute content.

22)
Static remobilization: nonlocal static remobilization process consists of the same contributing mechanisms of its dynamic counterpart. Time rate of static remobilization of immobile dislocations of type = , (̇r m ) is proportional to density of them: Given Eqs. (22) and (34): where rm is a material parameter associated with frequency of static remobilization of immobile dislocations of type = , . Moreover, static remobilization is a thermal dislocation process due to thermal character of its underlying mechanisms such as dislocation climb and junction dissociation.
Notice that pinning of already immobilized cell dislocations (cell immobile dislocations), strengthens their immobilization (resistance to remobilization). Consequently, upon application of external stress, just prior to local yielding, some of those pinned cell immobile dislocations will convert to GND due to elastic bending of their surrounding lattice to a sufficient mean misorientation angle.

Temperature and strain rate dependencies
24) Temperature dependence of dislocation processes: as pointed out earlier, each thermal dislocation process comprised of at least one underlying thermally-activated mechanism. Temperature is a statistical variable (representative of mean amplitude of atomic fluctuations) as well, which directly affects probability amplitude of different thermal dislocation processes. It is postulated that the change in probability amplitude of a thermal dislocation process with respect to temperature in constant strain rate is proportional to a power of temperature change: Therefore, where is absolute temperature; 0 is reference absolute temperature; is probability amplitude associated with thermal dislocation process that involves dislocations of type ; and 0 are respectively, current (at current temperature) and reference (at reference temperature and strain rate) material parameters associated with probability amplitude of dislocation process of dislocations of type ; and are respectively temperature sensitivity coefficient and exponent associated with dislocation process of dislocations of type . Reference temperature is assumed to be the minimum temperature in the temperature regime under consideration. Thus, in case of cold and warm regimes, reference temperature is the room temperature. Likewise, the reference strain rate is assumed to be the lowest strain rate in the investigated strain rate regime. The reference strain rate must be lower than the maximum/critical strain rate beyond which the viscoplastic deformation cannot be considered isothermal anymore because of adiabatic heat generation.
Referring to the postulates of dynamic and static dislocation processes (sections 3.2 and 3.3), thermal dislocation processes are dynamic annihilation of dislocations, dynamic trapping of cell mobile dislocation, dynamic nucleation of wall immobile dislocations, dynamic remobilization of immobile dislocation, static pinning of cell mobile dislocations, and static remobilization of immobile dislocations: and warm regimes). The temperature dependence of dislocation mobility (stress dependence of velocity of mobile dislocations) and underlying mechanisms of thermal dislocation processes, e.g. dislocation climb (speed), are often much more complex than that can be described purely by Arhenius relation and the respective activation energy barriers (Amodeo and Ghoniem, 1990;Argon and Moffatt, 1981;Blum et al., 2002;Eisenlohr and Blum, 2005;Gu et al., 2015;Hirth and Lothe, 1982;Yuan et al., 2018). In addition, the temperature dependence of some of the involved thermal mechanisms has not been properly understood yet. Therefore, it is reasonable to assume a phenomenological power-law description for temperature dependence of constitutive parameters associated with probability amplitude of different thermal dislocation processes.

25)
Strain-rate dependence of dislocation processes: the strain rate dependence of probability amplitude of ratedependent dynamic dislocation processes is suggested to be described by power-law relation as: where is strain rate sensitivity associated with dislocation process of dislocations of type . Therefore, at reference temperature: where ̇0 is reference shear strain rate; and ̃̇ is normalized mean plastic shear strain rate. Combination of Eqs. (73) and (76) produce: where ̂ is normalized absolute temperature; ̂ is normalized material coefficient associated with frequency of dislocation process that involves dislocations of type . Moreover, given Eqs. (4) and (76), and owing to Taylor factor ( ) being constant in the isotropic case: where ̇0 is reference strain rate. Furthermore, as pointed out, material coefficients associated with probability amplitudes of dynamic thermal dislocation processes depend on strain rate, although their rate sensitivity in quasi-static and intermediate regimes is often negligible. Therefore, { < 0 ∶ = an, tr ⋛ 0 ∶ = nc, rm = 0 ∶ = gn, ac, pn, rm 26) Temperature dependence of plastic stress: strain rate sensitivity of mean dislocation interaction strengths (̃) are reportedly negligible compared to strain rate sensitivity of viscous stress (Mecking and Kocks, 1981). Nonetheless, mean interaction strengths decrease with increasing temperature for each material in a characteristic way (Kassner, 2015;Kocks and Mecking, 2003). Shear modulus ( ) has a mild temperature dependence (decreasing with increasing temperature) as well which can be expressed by power-law relation (Argon, 2012;Galindo-Nava and Rae, 2016;Ghosh and Olson, 2002). Variation of temperature-dependent factors in Eqs. (26) and (28), ̃ and , with respect to temperature can be described simultaneously similar to temperature dependence relation of frequency of thermal dislocation processes: where 0 and ̃0 are respectively reference (at reference temperature) shear modulus and mean dislocation interaction strength associated with immobile dislocations of type ; ̃ is normalized ̃; and and are temperature sensitivity coefficient and exponent associated with ̃. Fig. 4, the viscous contribution of yield shear stress at slip systems is represented by a nonlinear dashpot in rheological representation of the elasto-viscoplastic (EVP) constitutive model. Stress response of the nonlinear dashpot is:

27) Viscous stress and its temperature dependence: as shown in
where ̃0 is mean viscous shear stress at reference strain rate; and is strain rate sensitivity parameter of viscous stress. At constant reference strain rate, by the assumption of ∆̃∝ (∆ ) , temperature dependence of viscous resistance can be expressed as: where ̃0 0 is the reference viscous shear stress (at reference temperature and strain rate); and and are respectively viscosity's temperature sensitivity coefficient and exponent. Eq. (82) is similar to the relations proposed by Kocks et al. (1975) and Argon (2012) where 00 is the reference viscous stress (at reference temperature and strain rate); and ̂ is normalized viscous stress. The reason for assuming an ambiguous phenomenological description for temperature dependence of viscous stress is that the temperature dependence of mean (mixed/curved) dislocations mobility function is extremely complex and has not been yet properly understood (Argon, 2012;Fleischer, 1962;Gilbert et al., 2011;Gilman, 1965;Kocks et al., 1975;Li, 1967;Tang and Marian, 2014).
28) Rate dependence of strain rate sensitivities: strain rate sensitivity of viscous drag increases with increasing temperature (Cereceda et al., 2016;Khan and Liu, 2012;Kocks, 1976;Rusinek and Rodríguez-Martínez, 2009) and is assumed to have a similar form as aforementioned temperature and strain rate dependence relations: where 0 is reference strain rate sensitivity (at reference temperature and strain rate); ̂ is normalized strain rate sensitivity of viscous stress; and are respectively temperature sensitivity coefficient and exponent associated with strain rate sensitivity of viscous stress; and is strain rate sensitivity parameter associated with strain rate sensitivity of viscous stress. Moreover, strain rate sensitivity coefficients associated with dislocation processes are assumed to have similar form of temperature and strain rate dependencies: where , 0 and ̂ are respectively, current, reference (at reference temperature) and normalized strain rate sensitivities associated with dislocation process of dislocations of type ; and and are temperature sensitivity coefficient and exponent associated with strain sensitivity of dislocation process of dislocations of type , respectively. Further, and are very small compared to and . Thus, only often shock regimes they have considerable impact.

29)
Plastic dissipation/adiabatic heating: volumetric adiabatic heat generation rate due to plastic work (̇) which is a fraction of volumetric plastic power (̇) (Taylor and Quinney, 1934), is obtained as follows: where ̇ is volumetric stored elastic power which is a fraction of volumetric plastic power that is stored in material by dislocations; and is referred to as dissipation/conversion factor, inelastic heat fraction, efficiency of plastic dissipation, or the Taylor-Quinney coefficient. In other words, is the fraction of plastic power that is not stored elastically and consequently is converted to heat. Taylor and Quinney (1934) emphasized that the fraction increases as plastic deformation progresses until saturation state where = 1. In other words, in saturation state, the entire input volumetric plastic power is converted to heat (̇=̇). Moreover, Rosakis et al. (2000) and Zehnder (1991) proposed models for variation of as a function of plastic strain and plastic hardening ( ≡ ) where approaches one by decreasing as plastic strain increases. At the beginning of plastic deformation of a nearly undeformed/annealed polycrystalline metallic material, large portion of the input plastic energy is stored in the crystal structure by generation of dislocations and dislocation structures. However, as the deformation progresses the generation rate of dislocations diminishes whereas the annihilation rate increases (annihilation releases the previously stored energy of dislocations in form of heat) until the saturation state where these two rates are equivalent. Therefore, one can assume: where is a material constant associated with dissipation factor. Furthermore, given Eqs. (12), (14), (28), (34), (78), (83) and (84), plastic/strain hardening ( ) can be calculated as follows: where is plastic hardening associated with dislocations of type . In addition, according to Eqs.
where ̇≡̇ is partial derivative operator with respect to plastic strain rate (̇); and are viscous hardenings respectively associated with viscous and plastic stress; and is viscous hardening associated with plastic stress of type .
For FE implementation, in return mapping (using Newton-Raphson scheme) of EVP constitutive models as well as calculation of viscoplastic consistent tangent operator (implicit FE framework), the viscoplastic tangent modulus ( ) is required. Given Eqs. Details of FE implementation of the above discretized constitutive equations in the framework of hypo-EVP finite deformation based on isotropic associative J2 plasticity is planned to be published soon.

Identification of constitutive parameters
Generally, constitutive/material parameters/constants of constitutive models are directly derived from experimental flow curves (yield stress versus plastic strain) through parameter identification techniques. Here, as an important hypothesis, it is envisaged that having the flow curves of a material in various temperatures, strain rates, etc. in the considered regimes is sufficient for calibration of a properly devised microstructural constitutive model using a robust parameter identification procedure. After all, under viscoplastic deformation, the underlying microstructural processes statistically result in a specific material response that is macroscopically observable through flow curves. However, although due to their physical nature, some of the constitutive parameters can be obtained by independent characterization methods (other than flow curves), it is suggested that the most accurate and effective approach for finding them is the simultaneous parameter identification, given the flow curves as the reference of fitting. Identification of microstructural parameters is carried out in two steps: 1) In the first step, isothermal uniaxial compression/upsetting or tensile tests are conducted in reference strain rate and various temperatures. Then, based on the measured experimental flow curves, most of the constitutive parameters are determined by means of pointwise parameter optimization using RMV-based analysis (RMVA). 2) In the last step, nonisothermal uniaxial tests are performed with different strain rates (in intermediate strain rate regime) and various initial temperatures (in cold and warm regime). Subsequently, the remaining constants which are strain rate sensitivity parameters ( ), temperature sensitivity coefficients ( ) and exponents ( ) associated with strain rate sensitivities, as well as the material parameter associated with dissipation factor ( ) are simultaneously calibrated by parameter optimization using finite element model updating (FEMU) method also known as finite element-based inverse strategy (Steenackers et al., 2007). After conducting uniaxial tests, force-displacement data is processed to true stress-strain, and then to flow curves. Afterwards, noise of flow curve data is reduced (here, "smooth" function in MATLAB library is used for noise filtration) and interpolated for constant plastic strain intervals. For the first step of parameter identification by RMVA, the following mean relative error which is also called objective, fitness, or residual function, is constructed as follows to be optimized/minimized: where is the set containing constitutive constants for fitting; is total number of experiments at different temperatures with constant reference strain rate ̇0; is total number of data points of isothermal experimental flow curve at temperature ; Δ is the constant plastic strain interval between data points; num is yield stress calculated numerically by the constitutive model at temperature which is corresponding to = Δ ; is the set of MSVs of -th step; exp is yield stress in isothermal experimental flow curve at temperature and plastic strain = Δ ; and are weighting factors respectively for isothermal flow curve at temperature , and data point corresponding to = Δ . In default, = 1 and = 1, unless some individual flow curve or ranges of plastic strain, have higher weight/importance ( , > 1) or lower importance ( , < 1) than the rest. For the example provided in the next section, Δ = 0.001 and = = 1. For many cases, finding the global solution/minimum to this optimization problem is difficult and requires complex and robust mathematical optimization methods which follow the steepest descent with the crossover functionality. Many methods are available for searching for the global solution such as gradient-free minimization methods and evolutionary algorithms including genetic algorithm. These methods exist in literature and their corresponding computer codes are developed and embedded in many commercial mathematical software. For the example provided in the next section, minimization of Eq. (106) for calibration of material parameters is carried out by a script written in MATLAB software. The optimization script employs MATLAB's global search class "GlobalSearch" and genetic algorithm "ga" which are coupled with its local minimization solver for constrained nonlinear multivariate functions "fmincon". The local minimization solver optionally utilizes various robust optimization algorithms such as interior point and reflective trust region.
For the second step of fitting based on FEMU, firstly the constitutive model must be implemented and programmed as a material subroutine based on hypo-EVP finite deformation using the constitutive parameters derived in the first step of calibration and with an initial guess for the constitutive parameters yet to be determined in the second step. Then, depending on the type of conducted experimental uniaxial tests, an FE model of either uniaxial compression or tension test must be constructed and linked to the material subroutine. Next, optimization of remaining constitutive parameters is done by iterative running of FE simulation of uniaxial tests with consecutives correction and update of remaining parameters. Running of FE simulation of uniaxial tests with updated parameters can be done by an optimization script that controls the pre-processing and post-processing of FE simulation. However, with a good initial guess for the remaining material constants, the second step of optimization can be done manually as well. For the example provided in the next section, the presented constitutive model has been implemented as implicit and explicit user-defined material subroutines for hypo-EVP finite deformation in ABAQUS FE package based on isotropic associative J2 plasticity.

Results, discussion and validation
In order to determine constitutive parameters, and also for validation of the constitutive model, uniaxial compression tests at various temperatures and strain rates using a deformation-type dilatometer (DIL-805A/D by TA Instruments) were conducted as dilatometer provides better precision than conventional compression test machines. Temperature measurement was done by thermocouple welded on the middle point in height direction at the lateral surface of the compression specimen. Temperature and strain rate were controlled by the dilatometer's integrated PID controller. Moreover, GND density measurements using high resolution EBSD for deformed compression specimens at different plastic strains are carried out. The material used in dilatometer compression is case-hardenable ferritic-pearlitic steel alloy 20MnCr5 (1.7147/1.7149) which is widely used in industrial forging of automotive components such as bevel gear. The chemical composition of 20MnCr5 used in experiments is presented in Table 1. The material constants required by the constitutive model that are independently measurable are also presented in Table 2. As already mentioned, reference temperature and strain rate are selected as the minimum possible temperature and strain rate in the considered domain (Table 2). Since for the steel studied here, the minimum temperature for onset of dynamic recrystallization of ferrite is slightly above 600 °C, the maximum temperature for constitutive modeling in warm regime is selected to be 600 °C. For the first step of parameter identification, experimental dilatometer compression tests at reference strain rate and 13 different temperatures (20, 50, 100, 150, ..., 600 °C) are performed. For each temperature at least two compression tests are carried out. Cylindrical dilatometer compression samples with diameter of 3.5mm and height of 5mm without lubrication pockets are used. However, lubricant is applied on contact surfaces. Subsequently, corresponding flow curves including yield point elongation are derived, smoothed, and interpolated in constant plastic strain intervals. Averages of the resultant flow curves that have the same nominal test parameters (temperature and strain rate) are calculated. These flow curves are then used for calculation of constitutive parameters using RMVA scheme. The constitutive constants that are determined in the first calibration step using RMVA method are presented in Table 3; and their corresponding temperature sensitivity coefficients and exponents are listed in Table 4. Mean relative error of the fitting is ̃= 0.39% which is in the same range of experimental scatter/precision (around 0.43%). constitutive parameters identified in the second step of calibration procedure by FEMU technique using nonisothermal uniaxial compression tests at different strain rates and various initial temperatures are presented in Table 5. Table 5 Constitutive parameters identified in the second step of parameter calibration by FEMU method. Temperature dependencies of material coefficients associated with frequencies of various thermal dynamic dislocation processes are shown in Fig. 6. For discussing the temperature dependence of different dynamic dislocation processes demonstrated in Fig. 6, it is useful to consider an intermediate temperature regime as the transition between cold and warm regimes which in case of the investigated material ranges from 200 ℃ to 400 ℃. As shown in Fig. 6, probability amplitudes associated with dynamic annihilation of cell mobile and wall immobile dislocations, dynamic nucleation of wall immobile dislocations and dynamic remobilization of immobile dislocations are nearly constant or have almost linear temperature dependencies in cold (0 ℃ < < 200 ℃) temperatures. However, they have appreciable nonlinear (exponential Arhenius-type) temperature dependencies in the intermediate (200 ℃ < < 400 ℃) and warm (400 ℃ < < 600 ℃) regimes. This is partly ascribed to different temperature dependencies of various thermal dislocation mechanisms in different temperature regimes that contribute unevenly to the aforementioned dynamic dislocation processes.
As Fig. 6 suggests, frequencies of annihilation of cell mobile and wall immobile dislocations as well as remobilization of immobile dislocations are intensified by transition from cold to intermediate and warm regimes. This can be in part explained by the well-accepted phenomenon that the principal recovery mechanism changes from cross-slip to climb by transition from cold to warm regime because climb is energetically more favorable in warm regime (Essmann and Mughrabi, 1979;Galindo-Nava et al., 2012;Kubin et al., 1992;Nix et al., 1985;Püschl, 2002;Rivera-Díaz-del-Castillo and Huang, 2012). As plotted in Fig. 6 (e), probability amplitude of dynamic remobilization of cell immobile dislocations has the strongest temperature dependence compared to the rest of thermal dynamic dislocation processes. At 600 ℃, it is roughly 230 times higher than room temperature. However, dynamic annihilation of wall immobile dislocations at 600 ℃ is only about 10 times higher than that of room temperature. As mentioned in postulates (13) and (17), dynamic annihilation and remobilization (dynamic recovery) processes share some of the underlying thermal mechanisms. Moreover, as shown in Fig. 6 (b), probability amplitude of dynamic annihilation of cell immobile dislocations is almost linear in the entire temperature domain (cold and warm regimes) while the probability amplitude of the other dynamic annihilation processes have nonlinear temperature dependencies ( Fig. 6 (a) and (b)). The reason is that the thermally-activated and diffusioncontrolled dislocation recovery mechanisms such as dislocation climb operating on a cell immobile dislocation in warm regime, more likely results in its remobilization than annihilation compared to its wall immobile counterpart. However, after being remobilized, this former cell immobile dislocation that is now another cell mobile dislocation can become annihilated by another cell mobile dislocation (dynamic annihilation of cell mobile dislocations). This also can partly explain why the frequency of cell mobile dislocation annihilation is more sensitive to temperature raise. Hence, the climb mechanism has much stronger effect on dynamic annihilation of wall immobile dislocations than dynamic annihilation of cell immobile dislocations. Combining this with the fact that the climb mechanism is controlled by vacancy diffusion which has a high nonlinear temperature dependence in warm regime, will explain the difference between forms of temperature dependencies of probability amplitudes of dynamic annihilation of different types of immobile dislocations. As inferred from Fig. 6 (e), the climb mechanism substantially influences cell immobile dislocations because of their immobility and particular loose local arrangement compared to wall immobile dislocations. In cold regime, as shown in Fig. 6 (e), the climb mechanism is not strong enough in order to notably impact the remobilization rate of cell immobile dislocations. However, as shown in Fig. 6 (b), together with cross-slip, its intensity in cold regime is sufficient to linearly affect the annihilation rate of cell immobile dislocations.
Probability amplitude of remobilization of cell immobile dislocations for the investigated material, as mentioned, grows intensively in warm regime by increasing temperature, while it is almost constant in cold regime. Accordingly, cell immobile dislocations and their associated dislocation substructures (IDBs) have relatively high stability in cold regime. However, their stability diminishes abruptly by transition from cold to warm temperatures. Hence, cell immobile dislocation density must be much lower during plastic deformation at warm temperatures than that of cold temperatures. From this it can be inferred that the studied metal tends to formation of dislocation sub-cells (bounded by IDBs) in cold regime while in warm regime it leans toward formation of dislocation cells (bounded by dislocation walls or GNBs).
Furthermore, temperature dependencies of shear modulus, interaction strengths, viscous stress, and strain rate sensitivity are depicted in Fig. 7.

Isothermal uniaxial tests
Experimental and computational flow curves (after optimization by RMVA) at different temperatures (isothermal condition) for compressive deformation at reference strain rate (̇0 = 0.01 s −1 ) are shown in Fig. 8. For the RMVA fitting, the experimental data are used only up to the accumulated plastic strain of = 0.4, because at larger strains the effect of friction between compression tools and specimen becomes dominant. This induces notable inhomogeneous distribution of MSVs and strain rate (in compression specimen) that is in contradiction with the assumptions of homogenous RMV and uniaxial deformation. Moreover, since experimental measurements carry notable amount of noise, for fitting of constitutive parameters as well as for comparison to corresponding computational curves, they have been smoothed. In each temperature regime (cold, intermediate and warm), one of the experimental flow curves plotted in Fig. 8 has not been used for parameter identification. Those experimental flow curves (isothermal condition at reference strain rate) that are associated with temperatures 150 °C, 350 °C and 550 °C, have been later compared to their computational predictions for validation.
Isothermal flow curves as shown in Fig. 8, are categorized in the aforementioned three temperature groups: cold, intermediate and warm regimes. As shown in Fig. 6 and Fig. 7, in cold temperature regime for the investigated material (20 to 200 °C), temperature dependent constitutive parameters are either almost constant or have approximately linear temperature dependence. However, in intermediate and warm regimes (200 to 600 °C), temperature dependence of most of the thermal material constants abruptly changes. Temperature dependence of flow curves is the result of two competing effects:  Dynamic thermal softening (DTS) due to temperature dependence of viscous stress, shear modulus, interaction strength, dynamic annihilation of dislocations, and dynamic remobilization of immobile dislocations.  Dynamic thermal hardening (DTH) because of DSA (by pinning mechanism), i.e. temperature dependence of dynamic trapping of cell mobile dislocations and dynamic nucleation of wall immobile dislocations. According to Fig. 6, Fig. 7 and Fig. 8, with respect to relative average difference among flow curves at different temperatures, in cold regime, DTS mostly due to temperature dependence of viscous stress, is clearly stronger than DTH. In cold regime by increasing temperature, DTS gradually decreases while DTH grows. This can be seen from the decreasing trend (by increasing temperature) of relative difference of yield stress in flow curves of cold regime. As shown in Fig. 8 (a), in constant plastic strains by increasing temperature the yield stress decreases due to DTS. However, the amount of reduction by increasing temperature also decreases because the competing DTH becomes stronger as the temperature continues to rise. Thus, the transition between cold and intermediate regimes occurs in the temperature at which DTS and DTH almost neutralize each other.
In intermediate regimes, DTH is slightly more dominant than DTS. As shown in Fig. 8 (b), in constant plastic strains by increasing temperature the yield stress increases as well due to stronger DTH in intermediate regime compared to DTS. Nonetheless, at the upper bound of intermediate regime, the amount of increase in yield stress by increasing temperature diminishes because this time the competing DTS is getting stronger as the temperature continues to grow. In warm regime, the dominant effect again becomes DTS which enhances more than DTH by increasing temperature, mainly due to dynamic annihilation and remobilization processes, mean interaction strengths and shear modulus. As shown in Fig. 8 (c), in constant plastic strains by increasing temperature the yield stress decreases with an increasing rate due to very strong DTS in warm regime that cannot be competed by DTH effect. This implies that, in warm regime by increasing temperature, even though the pinning mechanism is amplified but the DTS impact becomes much greater that macroscopically results in an accelerated reduction of yield stress by increasing temperature. Experimental and computational plastic hardening curves (plastic hardening vs. yield stress or plastic strain) of compression tests in constant reference condition ( = 20 ℃ and ̇0 = 0.01 s −1 ) are demonstrated in Fig. 9. As plotted in Fig. 9, the plastic hardening increases abruptly at the beginning of deformation of the investigated material at reference temperature and strain rate. As shown in experimental flow curves depicted in Fig. 8, in cold regime and close to lower bound of intermediate regime, yield point elongation occurs which corresponds to this increasing domain in plastic hardening curves (Fig. 9) at small plastic strains. The proposed constitutive model has captured this effect. After reaching a maximum value, plastic hardening begins to decrease with a high rate. The rate of decrease of plastic hardening gradually decreases as the plastic strain continues to accumulate. As shown in Fig. 9 (a), the model predict saturation ( = 0) at yield stress of sat ≈ 900 MPa which corresponds to (accumulated) plastic strain of sat ≈ 1.1. The lowest plastic strain at which saturation occurs ( sat ) usually has a decreasing trend by increasing temperature. As shown in Fig. 8 (c), for the investigated material in warm regime, the saturation starts at plastic strains of less than 0.4 ( sat < 0.4) which can be measured experimentally. However, for materials similar to the one investigated in the present paper (nearly all steel grades), in cold and intermediate regimes, generally it is not possible to experimentally measure the flow curves until the point of saturation with acceptable accuracy (using uniaxial tests). Therefore, in order to perform FE simulation of metal forming processes in cold and warm regimes with adequate precision, one has to predict/extrapolate the flow curves generally for relatively very large plastic strains (until saturation) for which there are not experimental measurements available. The microstructural constitutive model proposed in this paper can be trusted for such predictions because it can reproduce the experimental flow curves very accurately in the plastic strain range for which the experimental data is available. On the other hand, it has a comprehensive (statistical) physical background that covers all the major dislocation types and processes in cold and warm regimes.

Evolution of dislocation densities at constant strain rate
Numerically calculated evolution of dislocation densities based on RMV versus plastic strain in isothermal condition at various temperatures and reference strain rate (̇0 = 0.01 s −1 ) are plotted in Fig. 10. It is emphasized that by increasing temperature in cold and warm regimes, there is trend of reduction in cell mobile, cell immobile and total dislocation densities. However, similar to flow curves of intermediate temperature regime (Fig. 8 (b)), by increasing temperature wall immobile dislocation density increases, but again in warm regime it decreases. This is reasonable because the plastic stress associated with wall (immobile) dislocations ( ) has the largest contribution to the yield stress (postulates (5) and (8)). As demonstrated in Fig. 10 (b) and (c), in cold regime (0 ℃ < < 200 ℃), generally (nonlocal) cell immobile dislocation density is higher than wall immobile dislocation density ( > ). In other words, as mentioned earlier in this section, this material tends to formation of cell immobile dislocation tangles (IDBs) at cold regime (sub-cell forming material). Nevertheless, by increasing temperature from 200 °C to 300 °C (intermediate regime), there is a sudden drop in cell immobile dislocation density and a (increasing) jump in wall immobile dislocation density. In this material, in the intermediate regime, (nonlocal) wall immobile dislocations have higher density than their cell immobile counterparts ( < ). This implies that, in warm regime, dislocation walls (GNDs) are becoming more frequent in the microstructure (of the investigated material) than cell dislocation pile-ups (IDBs). Temperature-dependent dynamic pinning associated with DSA as an underlying mechanism facilitating wall nucleation process (transforming cell dislocations to wall dislocations), is responsible for this effect. Due to contribution of dynamic pinning, as shown in Fig. 6 (d), at around 200 °C, frequency of wall nucleation process in the investigated material increases exponentially by increasing temperature. This results in transformation of more cell immobile to wall immobile dislocations. As Fig. 10 (d) suggests by increasing temperature, total dislocation density always decreases during plastic deformation at constant temperature and strain rate. At 600 °C which is close to the upper bound of warm regime for the investigated steel and therefore very close to the temperature of the onset of dynamic recrystallization (DRX), total dislocation density during deformation is almost constant. Furthermore, as shown in Fig. 10 (a), at the beginning of viscoplastic deformation, the undeformed/annealed material (low initial dislocation density), generates (cell) mobile dislocations with a remarkably high rate. High generation rate of mobile dislocations subsequently will lead to a maximum cell mobile dislocation density which occurs at relatively low accumulated plastic strains. This maximum is followed by a rapid drop in cell mobile dislocation density, very similar to the plastic hardening behavior ( Fig. 9 (b)). As the plastic deformation proceeds, the rate of production of cell mobile dislocations becomes gradually lower until the saturation state at which the rate of production of cell mobile dislocations becomes equal to the rate of their reduction.

GND density measurements
In order to examine and validate the model's prediction of GND density which as mentioned in postulate (3) is equal to wall immobile dislocation density ( ), several compression specimens are deformed to different plastic strain at constant reference temperature and strain rate. Along the symmetry axis of a deformed compression specimen, there is a gradient of accumulated (equivalent) plastic strain ranging from a near zero value at the contact surface to a maximum of more than mean/macroscopic/true accumulated plastic strain of the deformed specimen at its center point. Hence, along the symmetry axis there is always a point at which equivalent plastic strain is equal to the mean plastic strain of the deformed specimen. For each of those deformed specimens, coordinates of such points are derived through corresponding FE simulations. It must be noted that there are other points with this property in the cross section of deformed compression specimens; however, the strain state at those points are much more complex (triaxial) compared to the points on the symmetry axis which are almost under uniaxial strain state. Moreover, it can be easily shown that those material points with the aforementioned property ( ̅ = ), have also the same average and instantaneous equivalent plastic strain rate as the (constant) prescribed mean/macroscopic/true plastic strain rate of the deforming specimen ( ̅ =̇).
For each specimen, around the vicinity of such points (RMV), an EBSD sample is prepared. Sample preparation for EBSD involved standard mechanical polishing to 0.05 μm, followed by electropolishing in a 5% perchloric acid and 95% acetic acid solution (by volume) with an applied voltage of 35 V. Measurements are performed using a field emission gun scanning electron microscope (FEG-SEM), JOEL JSM 7000F, at 20 KeV beam energy, approximately 30 nA probe current, and 100-300 nm step size. A Hikari EBSD camera by Ametek-EDAX, in combination with the OIM software suite (OIM Data Collection and OIM Analysis v7.3) by EDAX-TSL, is used for data acquisition and analysis. Subsequently, at each point, GND density is calculated from kernel average misorientation (KAM) which is the average angular deviation between a point and its neighbors in a distance twice the step size as long as their misorientation does not exceed 5° (Calcagnotto et al., 2010). After mapping KAM values to GND density, over a representative area with the size of 100 × 100 μm, the average GND density is calculated.
Comparison between computationally predicted wall immobile dislocation density and the measured average GND density using high resolution EBSD in different plastic strains at reference condition ( = 20 ℃ and ̇0 = 0.01 s −1 ) is shown in Fig. 11.   Fig. 11. Comparison between computationally predicted wall immobile dislocation density and measured GND density using high resolution EBSD in different compressive plastic strains at reference condition ( = 20 ℃ and ̇0 = 0.01 s −1 ) along with corresponding EBSD images with GND density distribution (calculated from KAM data) over the representative area.

Multistep uniaxial tests
As mentioned in postulate (1), "Two identical material samples that are plastically deformed to an equal amount of accumulated plastic strain but with different histories of temperature and strain rate, if again deformed under an equal temperature and strain rate condition, do not necessarily yield the same stress response". In order to test this statement and further validate the history dependence of the presented microstructural constitutive model, the multistep compression tests are devised. In these tests, first, previously undeformed compression specimens are plastically deformed to a predefined plastic strain in constant temperature with low reference strain rate to make sure temperature remains almost constant during compression (initial compression step). In the next step (final compression step), they are cooled down and again deformed in constant room temperature at the reference strain rate to a certain plastic strain. This cycle is illustrated in the schematic time-temperature diagram shown in Fig. 12 (a). The initial compression step is performed at three constant temperatures, = 20, 300 and 400 ℃ until the accumulated plastic strain = 0.2. In Fig. 12 (b), flow curves derived from the final compression steps are compared to their corresponding RMV-based computational flow curves predicted by the constitutive model.
All three flow curves shown in Fig. 12 (b) are material responses (yield stress) under plastic deformation at identical constant reference temperature and strain rate with the same initial accumulated plastic strain 0 = 0.2. However, for each of them, the initial plastic strain was accumulated in different temperatures (different histories of plastic strain accumulation). As shown in Fig. 12 (b), initial yield stress is slightly higher while plastic hardening is considerably much lower for = 300 and 400 ℃ compared to those of = 20 ℃. The flow curves related to = 300 and 400 ℃ in the investigated plastic strain domain exhibit yield point elongation (Lüders bands) as they are convex (́≡

Nonisothermal uniaxial tests
Nonisothermal compression tests and their respective FE simulations at various strain rates (in intermediaterate regime) and different initial temperatures (in cold and warm regimes) are conducted for the second step of parameter identification to determine constitutive parameters associated with strain rate sensitivity and dissipation factor. In order to do so, first the presented constitutive model is implemented as a user-defined material subroutine UMAT in ABAQUS/Standard with the material constants identified in the RMV-based fitting step with an initial guess for the parameters yet to be determined in the second step of parameter identification using FEMU technique. Then, thermo-mechanical FE model of compression test with the nominal dimensions of experimental specimens and tools has been created in the implicit ABAQUS/Standard software. Mechanical and physical properties of compression material and tools such as elastic modulus, poisson's ratio, mass density, specific heat capacity, thermal conductivity, and thermal expansion are inputted to the FE model as functions of temperature. Convection and radiation along with thermal contact conductance as a function of pressure of contact interface of compression specimen and tools are also considered. In addition, simple coulomb friction model with friction coefficient of 0.05 corresponding to the experimental condition was used. As mentioned, details of FE implementation of the constitutive model and thermo-mechanical FE simulation of an industrial multistep warm forging of a bevel gear for the same material investigated here are planned to be published later. Fig. 13 shows FE-simulated distribution of statistical TMM variables, equivalent stress ( ̅), temperature ( ), equivalent plastic strain rate ( ̅ ), cell mobile dislocation density ( ), cell immobile dislocation density ( ), and wall immobile dislocation density ( ), in radial cross section of compression specimen during compressive plastic deformation at mean plastic strain of = 0.6 and strain rate of ̇= 1 s −1 with initial temperature of 0 = 20 ℃. Experimental and FE-simulated yield stress and temperature change (Δ = − 0 ) due to plastic work for compression tests at different strain rates and initial room temperature are shown in Fig. 14. Comparison between experimental and FE-simulated yield stress and temperature change due to plastic work (Δ = − 0 ) for compression tests at different strain rates (̇= 0.1, 1, 3, 10 s −1 ) and initial room temperature ( 0 = 20 ℃).
In each experimental compression test corresponding to Fig. 14, temperature is measured via a thermocouple welded on the middle point in height direction at the lateral surface of the compression specimen. Likewise, in FE simulation the temperature is read from the respective point. The experimental data (force and temperature versus displacement) obtained from the nonisothermal compression tests at constant strain rates of 0.1, 1.0 and 10.0 s -1 are supplied to the second step of parameter identification by FEMU procedure. As shown in Fig. 14, for validation of the model's rate-dependent features, the experimental data corresponding to the nonisothermal compression test at constant strain rate of 3.0 s -1 is compared to its counterpart predicted by the FE simulation.
As inferred from Fig. 14, at relatively high strain rates in the investigated rate regimes, since there is not enough time for the generated adiabatic heat to dissipate completely, it remain in the material; and as a result elevates the temperature which in turn leads to spontaneous DTS. Nonisothermal FE simulations of metal forming processes highly depend on the input thermo-physical properties of billet/blank and tools, their pressure-dependent friction and thermal contact conductance, as well as convection and radiation heat transfer with the ambient environment which all of them are temperature-dependent. Hence, large portion of the difference in temperature readings from these FE-simulated nonisothermal compression (virtual) tests compared to their experimental counterparts is attributed to the complexity of the thermomechanical FE simulation of metal forming that is independent from the constitutive model. In addition, the accuracy of experimental temperature measurement with welded thermocouples is highly dependent on the weld quality and position, as well as the thermocouple delay. However, the discrepancy between FE-simulated and experimental temperature increase is much less than expected.

Outlook
As every continuum plasticity model has limitations, despite its vast scope and outstanding accuracy, the presented constitutive model in this paper is not an exception: 1) It has been tested thus far only for cold and warm regimes. Therefore, application or extension of it to hot temperature regime demands further research. 2) It has been validated up to now, merely for the intermediate-rate regime. Its validity for other strain rate regimes as well as static state must be verified. Moreover, its overall strain rate dependence needs more tests. 3) In its current state, the model does not account for anisotropy and Bauschinger effect. In order to generalize it to account for such effects, some new postulates need to be added while probably some of the existing ones require slight modification. Subsequently, for validation, the model must be implemented in the crystal plasticity framework. 4) It does not account for strain (rate) localization due to pressure-dependent evolution (nucleation and growth) of void volume fraction. Incorporation of strain localization in the constitutive model is crucial for its application to simulation of metal forming processes, as forming operations of especially sheet metal usually need to be designed in such a way to avoid or minimize localized shearing, as this non-uniform deformation is believed to induce notable damage and microstructural heterogeneity in the final product, which is generally considered detrimental to its in-service performance. Furthermore, the introduced microstructural constitutive model requires further validation using different alloys with various compositions, microstructures and grain size distributions.

Conclusion
A fully coupled thermo-micro-mechanical constitutive model for isotropic viscoplasticity of polycrystalline metallic materials was proposed based on continuum dislocation dynamics. The constitutive model was designed in order to statistically capture and quantify all statistically-considerable microstructural phenomena influencing dislocations in cold and warm regimes and hence accurately describe macroscopic viscoplastic response of material. The following conclusions are drawn from the present approach:  The mechanical response of a macroscale material point under viscoplastic deformation is an implicit function of its stochastic microstructural state, plastic strain rate, and temperature. The presented constitutive model is history dependent and it does not depend on accumulated plastic strain. Multistep compression tests were conducted for verification of the accuracy of model's history dependence features.  Dynamic and static evolutions of different types of dislocations are functions of dynamic and static pinning of dislocations by impurity atoms. Hence, the proposed constitutive model accounts for the effects associated with pinning such as yield point elongation, dynamic strain aging, static aging and bake hardening. However, validation of static evolution of dislocation densities suggested by four postulates is a subject of another paper.  The constitutive model predicts flow curves at various temperatures and strain rates in cold and warm regimes with an exceptionally high accuracy (mean error of less than 0.4% which is in the same range of experimental scatter). Moreover, the evolution of measured GND density (by EBSD) is very close to that of predicted by the model.  Comparison of FE simulations of nonisothermal compression tests at various strain rates with their corresponding compression experiments verifies that the postulates suggested for model's rate sensitivity and varying dissipation factor in adiabatic heat generation (due to plastic work) are realistic.  Implementation of the presented constitutive model as user-defined material subroutines UMAT and VUMAT respectively, in ABAQUS/Standard and ABAQUS/Explicit based on associative isotropic J2 hypoelastoviscoplasticity, has revealed that the constitutive model is computationally efficient (simulation costs are in the same level as empirical models) in thermo-mechanical FE simulations of cold and warm metal forming processes. Mughrabi, H., 2006. Deformation-induced long-range internal stresses and lattice plane misorientations and the role of geometrically necessary dislocations. Philosophical Magazine 86,