Mechanisms-Based Transitional Viscoplasticity

: When metal is subjected to extreme strain rates, the conversation of energy to plastic power, the subsequent heat production and the growth of damages may lag behind the rate of loading. The imbalance alters deformation pathways and activates micro-dynamic excitations. The excitations immobilize dislocation, are responsible for the stress upturn and magnify plasticity-induced heating. The main conclusion of this study is that dynamic strengthening, plasticity-induced heating, grain size strengthening and the processes of microstructural relaxation are inseparable phenomena. Here, the phenomena are discussed in semi-independent sections, and then, are assembled into a uniﬁed constitutive model. The model is ﬁrst tested under simple loading conditions and, later, is validated in a numerical analysis of the plate impact problem, where a copper ﬂyer strikes a copper target with a velocity of 308 m / s. It should be stated that the simulations are performed with the use of the deformable discrete element method, which is designed for monitoring translations and rotations of deformable particles.


Introduction
Metals subjected to extreme dynamics experience rapid microstructural evolution. Dislocations are generated, become entangled and form fine structures [1]. As a result, microscopic plastic flow is frequently interrupted and rerouted. The complex motion of dislocations at quasi-static conditions becomes increasingly sophisticated at high strain rates. When energy is delivered to metals with rates faster than the rate at which the energy is converted to plastic work and damages, then there is uncompensated energy, which is partly stored in the newly created dislocation structures and the rest of it activates micro-dynamic excitations. Among direct consequences of the excitations are dynamic strengthening [2][3][4] and the magnification of plasticity-induced heating [5]. In the Taylor-Quinney interpretation, a significant portion of plastic work is converted to heat [6][7][8][9], while the remaining work is stored in dislocation structures. Dislocations nucleate, move through the material, are piling up and become entangled. Thus, the majority of plastic work aids configurational entropy of the material, while plasticity-induced heating quantifies the efficiency of the reconfigurations. The grain size dependence is an integral part of the picture. Hall and Petch along with Armstrong [10][11][12][13] explored the grain size sensitivities and proposed the generally accepted Hall-Petch relation. The relation successfully survived for nearly sixty years. We are learning that this relation does not work in nanocrystalline metals [14,15]. In fact, when grains are smaller than 15-20 nanometers, the Hall-Petch effect is reversed [16,17]. The relation is not quite applicable to metals composed of very large grains. Thus, there is still much to learn about the phenomena.
The proposed concept is quite different from the celebrated strength models introduced by Zerilli and Armstrong [18], Johnson and Cook [19], Follansbee and Kocks [20], Preston, Tonks and Wallace [21], among others [22]. The models are based on a salient assumption that the yield stress carries all H p results from slip events along many θ planes defined by two orthogonal unit vectors, n θ and s θ . The planes are somewhat misoriented with respect to the plane of the dominant slip ( Figure 1). Crystal-to-crystal misorientations, grain boundaries and other obstacles broaden the distribution of active slip systems. When slip is arrested along one plane, other less favorable pathways are created. At advanced deformation, the rerouting is linked to dislocation interactions, includes grain rotations, causes lattice distortions, and the processes are generally dynamic. A material point in a continuum description contains a large number of grains and, given grain-to-grain misorientations, there is a much larger number of slip planes. While the distribution remains discrete, the continuum approximation seems to be a justifiable assumption. For this reason, one can assume that the θ-planes experience continuous reorientations, such that (1) Crystals 2020, 9, x FOR PEER REVIEW 3 of 18 boundaries and other obstacles broaden the distribution of active slip systems. When slip is arrested along one plane, other less favorable pathways are created. At advanced deformation, the rerouting is linked to dislocation interactions, includes grain rotations, causes lattice distortions, and the processes are generally dynamic. A material point in a continuum description contains a large number of grains and, given grain-to-grain misorientations, there is a much larger number of slip planes. While the distribution remains discrete, the continuum approximation seems to be a justifiable assumption. For this reason, one can assume that the θ-planes experience continuous reorientations, such that The dominant slip is defined at = 0. At each angle, slip events occur with frequency ( ), where the angle is taken between ± 4 ⁄ . The distribution ( ) = and Γ is a gamma function. As shown in Figure 1, the exponent controls the shape of the distribution. When is approaching zero, slip occurs along the dominant direction. On the other hand, when is a large number, slip is distributed indiscriminately in all orientations between ± 4 ⁄ . The rate of macroscopic plastic strain ̇ results not only from local slip events, but also reflects the process of slip reorganizations (Equation (1)), hence The volumetric strain rate ̇ quantifies slip incompatibilities in tension. Plastic strain is assigned to each -slip plane. The -strains are weighted ̇ = ( ) ̇ , where ̇ is the rate of macroscopic strain. Since the reorientations are local events, therefore ̇= ( ) ⁄ ̇ . The coefficient makes the reorganization a thermally activated process. The concept of thermal activation is introduced in the next section. Slip systems aligned closely with the dominant direction ( = 0) do not rotate [28]. The slip reorganizations are maximized at = ± 4 ⁄ . The relation ⁄ = 2 sin 2 properly reflects the trend. The dominant slip plane is defined in terms of two orthogonal unit vectors and and the flow tensor is a dyadic product of the two = ( + ). Since plastic flow is stimulated by stress , therefore, the tensor must be expressed in terms of the current stress = ( ). The representation = + + of the dyadic product is constructed on the basis of stress deviator = − 3 ⁄ and the three variables , and are functions of stress invariants. The representation must reproduce invariants of the original flow tensor ( + ), thus = 0, ( ) = 2 and = ( ) . There are three The dominant slip is defined at θ = 0. At each angle, slip events occur with frequency f θ (θ), where the angle is taken between ±π/4. The distribution f θ (θ) = 2 Γ(1+1/2r) √ π Γ(1/2+1/2r) cos 1/r 2θ is constructed in such a manner that π/4 −π/4 f θ (θ)dθ = 1 and Γ is a gamma function. As shown in Figure 1, the exponent r controls the shape of the distribution. When r is approaching zero, slip occurs along the dominant direction. On the other hand, when r is a large number, slip is distributed indiscriminately in all orientations between ±π/4.
The rate of macroscopic plastic strain . H p results not only from local slip events, but also reflects the process of slip reorganizations (Equation (1)), hence The coefficient A T makes the reorganization a thermally activated process. The concept of thermal activation is introduced in the next section. Slip systems aligned closely with the dominant direction (θ = 0) do not rotate [28]. The slip reorganizations are maximized at θ = ±π/4. The relation e θ ∂θ/∂e θ = 2 sin 2θ properly reflects the trend. The dominant slip plane is defined in terms of two orthogonal unit vectors n σ and s σ and the flow tensor is a dyadic product of the two N = (n σ s σ + s σ n σ ).
Since plastic flow is stimulated by stress σ, therefore, the tensor must be expressed in terms of the current stress N = N σ (σ). The representation N σ = α 0 1 + α 1 S + α 2 S 2 of the dyadic product is constructed on the basis of stress deviator S = σ − 1 trσ/3 and the three variables α 0 , α 1 and α 2 are functions of stress invariants. The representation must reproduce invariants of the original flow tensor (n σ s σ + s σ n σ ), thus trN σ = 0, tr(N σ ) 2 = 2 and N σ = (N σ ) 3 . There are three solutions, where the representations are constructed on the planes of principal stresses {σ 1 , σ 2 }, {σ 2 , σ 3 } and {σ 1 , σ 3 } and where σ 1 > σ 2 > σ 3 . I select the solution which reproduces the maximum shear mechanism [29], hence The angle ϕ = sin −1 J 3 √ 27/4/J 3 2 and two invariants J 2 and J 3 of stress deviator S complete the relation. It turns out that N σ : σ/2 is the maximum shear stress. It is worth mentioning that several tensor representations have been developed and some of them are reported in [30]. In the current (Tresca) description, the θ-vectors n θ and s θ are co-rotational with the dominant slip plane n θ = n σ cosθ + s σ sinθ and s θ = s σ cosθ − n σ sinθ. As a result, the plastic flow (Equation (2)) becomes where the average Schmid factor M = M s − M r specifies the current slip arrangement M s and captures the contribution of the flow reorganizations M r . The term M s = Γ 2 (1+ 1 2r ) is closely approximated by M s 1+r 1+π r/2 , where the error of approximation does not exceed one percent. The second term M r = 4 A T (T) π(1+2r) is takes as derived. Note that the Schmid factor is a function of the exponent r. Advances of plastic deformation broaden the spectrum of the active planes and, therefore, r = r(e pe ). The exponent r = r 0 1 + e pe /e p r reflects the slip organization, where r 0 and e p r are constants. Here, e pe is the magnitude of effective plastic strain.
As stated earlier, ductile damage results from slip incompatibilities [31,32]. In tension, the compatibility is restored by the stress-favorable nucleation and growth of voids. The damage predominantly occurs on the plane of maximum shear, as shown in [24], and is The rate of damage is controlled by the parameter q.

Thermal Activation
Thermal activation (TA) in fcc metals results from the dislocation interactions, where the intersections of dislocations are the controlling mechanism [33]. In bcc metals, TA is determined by Peierls' mechanism, where much stronger temperature dependence of the yield stress is observed at low temperatures. In a continuum-level description, the activation processes can be described in terms of free activation energy F a = E a − S a T, where E a and S a represent activation energy and activation entropy [34]. Note that the entropic contribution S a T weakens the energy barriers F a , where entropy S a is nearly independent of temperature. I introduce a non-dimensional function ξ a = F a /k b T, where k b is Boltzmann constant and T is temperature. Thermal activation is stochastic with respect to ξ a . When dislocation glide is the dominant mechanism, the temperature has a minor influence on the formation of dislocation structures [25]. The distribution of dislocation avalanches is shown to be independent of temperature [26]. For these reasons, temperature in ξ a is just a scaling factor. One can visualize a local landscape of activation energy that contains energetically weakest sites. In the state of self-organized criticality, some of the sites expire and new ones are created, thus it is an evolving landscape. A macroscopic material point contains a very large number of the sources and, therefore, a continuum-type distribution would be an acceptable assumption. The Weibull distribution f T = k a ξ k a −1 a exp −ξ k a a seems to fit the description. The material is expected to increase the number of activation sites with frequency f T (T) and the exponent k a determines the rate of the process. The number of the sites at temperature T includes all the sites which are theoretically available up to this temperature. Therefore, the actual thermal resistance to plastic flow A T is integrated over all frequencies taken in the thermal domain that has not been explored, that is between T c /T m and T c /T, hence A T = T c /T T c /T m f T (s)ds. I introduce the transition temperature T c . At temperatures above T c , the dislocation mobility becomes more probable. Here, activation energy E a is expressed in terms of k b T c such that E a = g 1/k a a k b T c . Activation entropy is estimated at melting point F a (T m ) = 0 and we obtain S a = E a /T m . The activation factor becomes where g a determines the magnitude of activation energy. At melting temperature, when A T (T m ) → 0 , the dislocation activities are replaced by diffusional flow. At cryogenic temperatures, the factor A T is slightly smaller than one and the material exhibits the strongest resistance to plastic deformation.

Rerouting of Plastic Flow and Consequences
At small scales, plastic flow is an erratic process made of irreversible bursts of slip and includes other forms of local instabilities. The sources are activated, deactivated, are spatially fluctuating and undergo continuous reorganizations [35]. As Brown stated [25], the dynamic reconfigurations drive the system toward the state of self-organized criticality. One can visualize the process by monitoring the trajectory of a selected material particle ( Figure 2). Activation entropy is estimated at melting point ( ) = 0 and we obtain = / . The activation factor becomes where determines the magnitude of activation energy. At melting temperature, when ( ) → 0, the dislocation activities are replaced by diffusional flow. At cryogenic temperatures, the factor is slightly smaller than one and the material exhibits the strongest resistance to plastic deformation.

Rerouting of Plastic Flow and Consequences
At small scales, plastic flow is an erratic process made of irreversible bursts of slip and includes other forms of local instabilities. The sources are activated, deactivated, are spatially fluctuating and undergo continuous reorganizations [35]. As Brown stated [25], the dynamic reconfigurations drive the system toward the state of self-organized criticality. One can visualize the process by monitoring the trajectory of a selected material particle ( Figure 2).

Figure 2.
Path rerouting activates dynamic excitations. Note that overstress is only partly stored in newly created dislocation structures.
Initially, the particle resides in position { } and stress is . Shortly later ( ), the particle is forced to change its trajectory and arrives at position { }, where stress is . Should the original path be available, the particle would move to position { } with velocity and the stress would be . In all the scenarios, equations of motion • = ̇ must be preserved. The particle acceleration and mass density are denoted by ̇ and . The rerouting alters tractions projected onto the surface that is normal to the direction of the particle velocity . is the momentum tensor. The path rearrangements affect particle acceleration ̇= ( ̇− ̇ ), while the overstress ( − ) is partly stored in dislocation structures ( Figure 1). The uncompensated stress in Figure 1 activates dynamic excitations. Given sufficient time, the overstress is fully relaxed. Initially, the particle resides in position {X} and stress is σ X . Shortly later (δt), the particle is forced to change its trajectory and arrives at position {x}, where stress is σ. Should the original path be available, the particle would move to position {z} with velocity υ z and the stress would be σ z . In all the scenarios, equations of motion ∇·σ = ρ . υ must be preserved. The particle acceleration and mass density are denoted by . υ and ρ. The rerouting alters tractions projected onto the surface ∂V 0 that is normal to the direction of the particle velocity m. The difference in tractions (σ z ·m σ·m) between positions {z} and {x} triggers stress perturbations δσ = l r

Dynamic Overstress
The perturbations are stretched over distance l r . Next, the volume V 0 is reduced to a material point and stress perturbations become δσ = l r . ψ, where ψ = ρ(δυ m + m δυ)/2 is the momentum tensor. The path rearrangements affect particle acceleration δ structures ( Figure 1). The uncompensated stress δS a in Figure 1 activates dynamic excitations. Given sufficient time, the overstress is fully relaxed.

Dynamic Overstress
The stress perturbations δσ slow down plastic flow .

H
The excitations give rise to the viscous overstress κ . ψ/ρ l r . The overstress can be expressed in terms of stress perturbations κ . ψ/ρ l r = R −2 k δσ/τ, where the resistance to plastic flow R k = l r /τυ s carries information on the size of the perturbed domain l r and specifies the relaxation time τ = κ/µ required for the equilibration of the dislocation structures. Here, υ s = µ/ρ is shear velocity and µ is shear modulus. As a result, we Once again, the resistance to plastic flow R k = υ r /υ s arises in the actively perturbed domain l r , such that υ r = l r /τ. The excitations are redistributed into the surroundings with velocity slower than shear velocity υ s , thus the resistance R k cannot be greater than one. In summary, the material description is given in terms of three relations Stress perturbations slow down the slip process, while the part of stress which exceeds the elastic limit is converted to viscous overstress. The relations are further rearranged The rate of effective (true) plastic strain is scaled by the fourth-order drag tensor P = I + R 2 k C/µ. In textured metals, the drag tensor P channels plastic deformation and damages along crystallographically favorable pathways. In an isotropic metal, the drag tensor P −1 = I/ 1 + 2 R 2 k and the rate of plastic strain k take much simpler forms. Comparison of Equations (6 2 ) and (7 2 ) C· k δσ/τ confirms that the viscous overstress in Equation (6) is acting on the plastic strain. Relaxation of the viscous overstress is best described by Maxwell's process where S a is the active overstress ( Figure 2). Rerouting of plastic flow is enabled by the excess of energy . W l = S a : . S a /2µ, where trS a = 0. During slow processes, a large portion of W l is dissipated. At high strain rates, W l explicitly affects the plastic flow, increases storage of energy, intensifies plasticity-induced heating and influences the damage mechanism.
The initial resistance to plastic flow R 0 k is a consequence of pre-existing defects. At dynamic conditions, the excess of energy W l further amplifies the resistance. Consistently with the description of thermal resistance (Section 3), the relation is constructed on the basis of the weakest link hypothesis. The constants G l and n r specify dynamic constraints. The process occurs along the lowest energy pathways, while the frequency of the rerouting events is controlled by the exponent n r . The dynamic resistance to plastic flow R k is responsible for the stress upturn ( Figure 3). It is shown that the upturn becomes much less pronounced in metals composed of fine grains, where the initial value of R 0 k is closer to one. A full justification of the result is provided in Section 7 (Equation (13)). In Figure 3, the data points marked in red are collected from multiple sources [5,[36][37][38][39][40][41][42][43][44]. Experiments performed on fine-grained tantalum [14] and copper [15] confirm the predictions constructed for microcrystalline copper.

Plasticity-Induced Heating
In 1934, Taylor and Quinney calculated plasticity-induced heating and arrived at the conclusion that about 90% of plastic work is turned into heat. Recent estimates contradict the assessment [45]. Moreover, plasticity-induced heat is found to be a function of plastic strain, strain rate and type of loading. The strongest generation of heat occurs at the initial stage of plastic deformation. Further increase of plastic strain slows down the process, and then, the rate increases again at an even more advanced stage of deformation. In the Taylor and Quinney interpretation, a part of plastic power is converted to hear ̇= ̇ , where the Taylor-Quinney coefficient quantifies the conversion process. I challenge the Taylor-Quinney interpretation and postulate that plastic work contributes to the increase of configurational entropy of the material, while heat is a measure of the process inefficiencies caused by the unavoidable rerouting of the deformation pathways.
The argument is constructed as follows. First, we note that stress perturbations = ̇ are expressed in terms of the plastic strain rate  Ductile damage results from slip incompatibilities. In tension, as stated earlier, kinematical compatibility is restored by the stress-favorable nucleation and growth of voids. In the Tresca description, the growth of damage is defined on the plane of dominant shear (N σ ) 2 = (n σ n σ + s σ s σ ) and is pe . Subsequently, the volumetric change is tr The relation was derived in [24]. The internal friction parameter q scales the damage process. When q = 1, voids nucleate and grow on the plane of maximum tensile stress. In compression, q is equal to zero. Here, the flow incompatibilities in tension q = 1 − e −e pe d /e d are controlled by ductility e d , where plastic strain e pe d is accumulated during tensile loading only. The damage is responsible for the degradation of the apparent bulk and shear moduli. The scaling parameter is 1 − e −H d kk /H 0 kk . Here, H 0 kk specifies the strain at which the damage process is macroscopically evident. Consequently, in a damage-free metal, the scaling The factor is equal to one in a completely fractured material.

Plasticity-Induced Heating
In 1934, Taylor and Quinney calculated plasticity-induced heating and arrived at the conclusion that about 90% of plastic work is turned into heat. Recent estimates contradict the assessment [45]. Moreover, plasticity-induced heat is found to be a function of plastic strain, strain rate and type of loading. The strongest generation of heat occurs at the initial stage of plastic deformation. Further increase of plastic strain slows down the process, and then, the rate increases again at an even more advanced stage of deformation. In the Taylor and Quinney interpretation, a part of plastic power is converted to hear . Q = βσ eq . e pe , where the Taylor-Quinney coefficient β quantifies the conversion process. I challenge the Taylor-Quinney interpretation and postulate that plastic work contributes to the increase of configurational entropy of the material, while heat is a measure of the process inefficiencies caused by the unavoidable rerouting of the deformation pathways.
The energy of the excitations is absorbed in dislocation structures. Still, there a fraction ξ e f of . ψ d which gives rise to phonon vibrations The term S a /2τ in Equation (8)  H pe and are affected by the flow constraints R k . As plastic deformation advances, heat sources continue to evolve ψ Q = ψ p + ξ e f R k C·H pe /υ s . In this expression, momentum ψ p is linked to shear stress capable of triggering path rearrangements.
Heat is produced when excitations . ψ act on the existing sources of heat ψ Q , such that The relation is further simplified e pe , where equivalent plastic strain e pe has already been defined and σ y is yield stress. The Taylor-Quinney coefficient β = .
Q/ . W p becomes β = 2R k R k ξ e f µ e pe + σ y /σ eq . The heat sources are formed discretely both in time and space, where spatial resolution of the sources is estimated to be in a sub-micrometer range. Experimental techniques are yet to be developed for the measurement of such small-scale temperature fluctuations. Large-scale atomistic simulations might provide invaluable insight into the phenomena. It should be stated that the contours in Figure 4a,b are constructed under the assumption that the copper specimen is subjected to uniaxial compression and the temperature rise is depicted at strain 20%. The coefficient ξ e f is calibrated using experimental data presented in [46]. In Figure 4a, the temperature rise is calculated as a function of grain size d and initial temperature T 0 /T m . The samples are deformed at a strain rate of 10 4 /s. The contours indicate that plasticity-induced heating is significantly stronger in microcrystalline copper. The environment T 0 affects the rate of heating. In large-grain copper, high strain rates noticeably magnify plasticity-induced heating, while small grains make the strain rate sensitivity an irrelevant factor (see Equation (13) and Figure 4b). For example, at a strain rate of 2.5 10 5 /s the same temperature rise is predicted in copper with average grains 2.5 µm and 77 µm. It would be very interesting and, perhaps, technologically important to experimentally validate the theoretical predictions. The grain size dependence is discussed in the next section.
In large-grain copper, high strain rates noticeably magnify plasticity-induced heating, while small grains make the strain rate sensitivity an irrelevant factor (see Equation (13) and Figure 4b). For example, at a strain rate of 2.5 10 ⁄ the same temperature rise is predicted in copper with average grains 2.5 m and 77 m. It would be very interesting and, perhaps, technologically important to experimentally validate the theoretical predictions. The grain size dependence is discussed in the next section.

Hall-Petch Relation
In early 1950s, Hall [10] and Petch [11] introduced a relation, which correlates yield stress with an average size of grains = + ⁄ .
(10) In this construction, is the size-independent strength, is a constant and is the average size of grains. In this notation, yield stress corresponds to equivalent stress = . The relation is generally applicable to metals and other polycrystalline materials. A few years later, Armstrong [12] proposed a dislocation pileup interpretation and further explained the mechanism in numerous writing [47]. In this construction, the parameter is considered a Griffith type stress intensity factor, where the grain size strengthening results from the dislocation pileup along slip bands blocked at a grain boundary, thus behaving similarly to a shear crack. Among other descriptions, it is worth mentioning the concept by Ashby [48]. In Ashby's formula, the constant includes the additional term ⁄ . It will be shown that the term is reproduced in the proposed energy-based formula.
Strong experimental evidence exists that dislocation pileups are responsible for the observed heterogeneous distribution of stresses. The supporting evidence was obtained in micro-diffraction (DAXM) experiments further supported by high-resolution electron backscatter diffraction (HR-EBSD) measurements [49]. It is worth noting that nanoindentation tests [50] have shown that the

Hall-Petch Relation
In early 1950s, Hall [10] and Petch [11] introduced a relation, which correlates yield stress σ y with an average size of grains d In this construction, σ p y is the size-independent strength, k y is a constant and d is the average size of grains. In this notation, yield stress corresponds to equivalent stress σ y = σ eq . The relation is generally applicable to metals and other polycrystalline materials. A few years later, Armstrong [12] proposed a dislocation pileup interpretation and further explained the mechanism in numerous writing [47]. In this construction, the parameter k y is considered a Griffith type stress intensity factor, where the grain size strengthening results from the dislocation pileup along slip bands blocked at a grain boundary, thus behaving similarly to a shear crack. Among other descriptions, it is worth mentioning the concept by Ashby [48]. In Ashby's formula, the constant k y includes the additional term µb 1/2 . It will be shown that the term is reproduced in the proposed energy-based formula. Strong experimental evidence exists that dislocation pileups are responsible for the observed heterogeneous distribution of stresses. The supporting evidence was obtained in micro-diffraction (DAXM) experiments further supported by high-resolution electron backscatter diffraction (HR-EBSD) measurements [49]. It is worth noting that nanoindentation tests [50] have shown that the material hardness increases in the proximity of grain boundaries. Thus, the experimental results convey a clear message that dislocation pileups near grain boundaries are indeed responsible for the stress and energy heterogeneities.
The Hall-Petch interpretation is based on the postulate that, at a given rate of plastic strain, the upsurge of plastic power is controlled by the increases of yield stress Here, the constraints R k magnify the effective yield stress σ y = σ p y 1 + 2R 2 k . When assuming that the Hall-Petch relation prevails, one can obtain R 2 k = k y d −1/2 /2σ p y . Next, I construct a kinematics-based interpretation. At the prescribed stress σ y , the resistance R k is responsible for slowing down the plastic flow. For this reason, we write . W p σ y = σ y M . e p / 1 + 2R 2 k . The two interpretations represent the upper and lower bound estimates of the size effect.

Energy-Based Hall-Petch Relation
I propose an energy-based interpretation of the Hall-Petch strengthening. As stated above, the Hall-Petch size effect is mostly caused by grain boundaries, which slow down the plastic flow.
Suppose that polycrystalline metal is subjected to stress σ p y . This stress activates slip inside grains, but the deformation is blocked at grain boundaries. Grain boundaries are energy barriers defined in terms of surface energy g gb . At stress σ p y , the barriers are partly climbed g y gb . As loading advances, blocked dislocations magnify the pileup stress. It should be stated that dislocations are emitted in other areas, such as grain boundary corners. After some time (the delay time), the pileups accumulate sufficient energy for a dislocation to cross the barriers g gb . In essence, it is a stochastic process. Stress concentrations appear in different locations and grow at different rates. Large grains contain a large reservoir of energy which, in turn, feeds the energy concentrations more efficiently. The pileup-stimulated increase of internal energy is taken per unit grain boundary area S gb and becomes E gr = V s /S gb σ y σ p y w(σ s ) σ s dσ s /2µ. The average shear stress σ s continues to rise and becomes the macroscopic yield stress σ y . At this point, dislocations cross the boundaries in large numbers. The average rate of energy accumulation slows down. The correction function w = σ y − σ p y /2σ s captures the slowing down effect. When energy E gr overcomes the energy barrier (g gb − g a gb ) S gb /V s , the macroscopic plastic flow is initiated Grain boundary area per unit volume S gb /V s characterizes granular structure of the polycrystal. This ratio is sensitive to the grain size, shape and number of grains in a sample. The Hall-Petch relation becomes In this expression, the Hall-Petch mechanism is controlled by grain boundary area per unit volume. The parameter k yS = 2 µ g gb − g a gb is closely related to constant k y in Equation (10), such that k yS k y /2. It has been shown [51] that the grain boundary energy is linearly dependent on shear modulus multiplied by a length parameter, where the parameter specifies the grain boundary type. Thus, the relation (g gb − g a gb ) = µ (α 2 yS b) suggests that k yS = 2 µ α yS b 1/2 and α yS carries information on the metal crystallographic structure. It turns out that the parameter α yS is a material constant in fcc and bcc metals. Based on the data collected in [52], the mean value of α yS for fcc metals is α yS = 0.04604 ± 0.00616. In bcc metals, the values are found to be α yS = 0.1089 ± 0.0353. Estimates are also made for hcp metals, but the data scatter is much larger. The hcp deformation twins might be responsible for the scatter. As usual, the magnitude of the Burgers vector is denoted by b.
In Equation (12), the grain boundary area per unit volume is the size variable. The variable can be estimated with the use of stereology methods [53]. In bulk samples composed of nearly equiaxed grains, the average size of grains is about d ∼ 4V s /S gb . However, the validity of the approximation can be questioned in samples with a small number of grains. The grain shape matters too. In one case, yield stress is significantly reduced in thin films of nickel, where the film's thickness is comparable with the size of grains [54]. In the case of mild steel, Armstrong [12] made similar observations. . The pre-existing resistance to plastic flow R k = R 0 k becomes a quantifiable material constant

Kinematics-Based Construction of Hall-Petch Relation
Now, it becomes clear that the resistance to plastic flow is strongly dependent on the size of grains. In microcrystalline copper, the value of R 0 k is approaching one, thus the size effect becomes less pronounced. As shown in Section 5, the resistance R k = l r /(υ s τ) is a function of relaxation time τ = l r /(R k υ s ). The length l r is directly related to the grain boundary area per unit volume and is l r = S gb /V s −1 . Based on the above, the relaxation time becomes the size-dependent material property It is worth stating that relaxation time τ (Equation (14)) matches earlier estimates made for OFHC copper [55,56]. In nanocrystalline copper, relaxation time is in the range of hundred picoseconds. The kinematically-based Hall-Petch relation fits well the existing experimental data in [52] (Figure 5). Similar trends are reported for other metals [52,57] as well.
At constant stress, energy barriers slow down plastic flow and, in this manner, affect plastic power ̇ = ̇ (1 + 2 ) ⁄ . In the energy-based construction, the term (1 + 2 ) in ̇ = The pre-existing resistance to plastic flow = becomes a quantifiable material constant Now, it becomes clear that the resistance to plastic flow is strongly dependent on the size of grains. In microcrystalline copper, the value of is approaching one, thus the size effect becomes less pronounced. As shown in Section 5, the resistance = ( ) ⁄ is a function of relaxation time = ( ) ⁄ . The length is directly related to the grain boundary area per unit volume and is = ⁄ . Based on the above, the relaxation time becomes the size-dependent material property It is worth stating that relaxation time (Equation (14)) matches earlier estimates made for OFHC copper [55,56]. In nanocrystalline copper, relaxation time is in the range of hundred picoseconds. The kinematically-based Hall-Petch relation fits well the existing experimental data in [52] (Figure 5). Similar trends are reported for other metals [52,57] as well.

Transitional Viscoplasticity
Plastic power : ̇ = ̇ is expressed in terms of equivalent stress = : and strain rate ̇ . During active deformation processes, plastic power ̇ aids to configurational entropy (suppleness [25]) of the material. The rate of true plastic strain ̇ = [ is defined in terms of the effective rate ̇ = ̇ , where the resistance-free strain rate ̇ is a function of equivalent stress and the relation is
The power-law viscoplasticity was first introduced in [58,59]. Strength σ p y = σ 0 A T (T) is determined at 0.2% strain, where σ 0 is athermal strength. Exponent n p determines the elastic-plastic transition. What makes the relation uncommon is that the strain rate sensitivity is further tuned by the rate factor  (Figure 1). Advances of plastic deformation broaden the spectrum of active planes and, therefore, r = r(e pe ). , where we should note that n p / 1 − ω p n p .

OFHC Copper
The constitutive model captures copper responses in a broad range of temperature and strain rates, predicts ductile damage, accounts for the size of grains, and includes the description of plasticity-induced heating. There are three groups of parameters. Parameters readily available in the literature are listed in Table 1. Among the constants are elastic properties, mass density, yield stress at 0.2% strain, Burgers vector, melting point and specific heat. Several parameters display a marginal variability. Among the parameters are the quasi-static strain-rate sensitivity ω p , stress exponent n p , thermal efficiency of plastic flow ξ e f , crystallographic constant α yS , the initial distribution of slip planes r 0 , the transition temperature T c and the rate of defect formation n r . The parameters are shown in Table 2. In the absence of experimental data, the constants should properly approximate responses of fcc metals. The material-specific parameters are listed in Table 3. Among them are the activation energy constants (g a , k a ), the strength pre-factor Λ 0 , critical energy G l , hardening strain e p r and damage parameters e d , H 0 kk . In this group, the activation energy constants show low variability and, therefore, there are only five fully tunable constants. One may argue that the total number of parameters is large. It is worth noting that the parameters have well-defined interpretations and most of them are material constants. I emphasize that the constitutive model captures material responses in a very broad range of temperature and strain rates, predicts ductile damage, accounts for the size of grains, and includes plasticity-induced heating. The model capabilities are summarized in the deformation map (Figure 6), where stress at 20% strain is plotted as a function of strain rate and temperature. Two maps are constructed. In Figure 6a, the average grain size is 80 µm. The model capabilities are summarized in the deformation map (Figure 6), where stress at 20% strain is plotted as a function of strain rate and temperature. Two maps are constructed. In Figure 6a, the average grain size is 80 m. The red points depict experimental data collected from several sources [5,[35][36][37][38][39][40][41][42][43][60][61][62] and the blue mesh represents the model predictions. The model correctly predicts copper responses (grains 80 m) at temperatures from cryogenic environments up to melting points and strain rates from quasi-static loading to extreme strain rates. The second plot in Figure 6b is constructed for microcrystalline copper, where the grain size is 1 m. Note that in the second case the stress upturn is nearly gone. This result matches experimental measurements performed on small-grained tantalum [14] and copper [15]. The strain rate sensitivity was also studied in nickel [63], where ultrafine-grained nickel exhibits a diminishing rate-sensitivity.
For completeness, I included selected stress-strain responses (Figure 7). The predictions (solid lines) are calculated at three temperatures and three strain rates. The experimental data points (red, blue and black dots) are collected from [37,61]. In all cases, I included the contribution of plasticityinduced heating, where specific heat varies with temperature. Figure 7. Selected stress-strain responses are constructed for three strain rates and three temperatures. The average size of grains is assumed to be 80 m. The data is collected from [37,61]. The red points depict experimental data collected from several sources [5,[35][36][37][38][39][40][41][42][43][60][61][62] and the blue mesh represents the model predictions. The model correctly predicts copper responses (grains 80 µm) at temperatures from cryogenic environments up to melting points and strain rates from quasi-static loading to extreme strain rates. The second plot in Figure 6b is constructed for microcrystalline copper, where the grain size is 1 µm. Note that in the second case the stress upturn is nearly gone. This result matches experimental measurements performed on small-grained tantalum [14] and copper [15]. The strain rate sensitivity was also studied in nickel [63], where ultrafine-grained nickel exhibits a diminishing rate-sensitivity.
For completeness, I included selected stress-strain responses (Figure 7). The predictions (solid lines) are calculated at three temperatures and three strain rates. The experimental data points (red, blue and black dots) are collected from [37,61]. In all cases, I included the contribution of plasticity-induced heating, where specific heat varies with temperature. The model capabilities are summarized in the deformation map ( Figure 6), where stress at 20% strain is plotted as a function of strain rate and temperature. Two maps are constructed. In Figure 6a, the average grain size is 80 m. The red points depict experimental data collected from several sources [5,[35][36][37][38][39][40][41][42][43][60][61][62] and the blue mesh represents the model predictions. The model correctly predicts copper responses (grains 80 m) at temperatures from cryogenic environments up to melting points and strain rates from quasi-static loading to extreme strain rates. The second plot in Figure 6b is constructed for microcrystalline copper, where the grain size is 1 m. Note that in the second case the stress upturn is nearly gone. This result matches experimental measurements performed on small-grained tantalum [14] and copper [15]. The strain rate sensitivity was also studied in nickel [63], where ultrafine-grained nickel exhibits a diminishing rate-sensitivity.

Activation Energy Factor
For completeness, I included selected stress-strain responses (Figure 7). The predictions (solid lines) are calculated at three temperatures and three strain rates. The experimental data points (red, blue and black dots) are collected from [37,61]. In all cases, I included the contribution of plasticityinduced heating, where specific heat varies with temperature. Figure 7. Selected stress-strain responses are constructed for three strain rates and three temperatures. The average size of grains is assumed to be 80 m. The data is collected from [37,61].

Plate Impact Problem
The constitutive model is implemented to a code of deformable discrete element method (DDEM). In DDEM, each discrete element represents a crystallographic grain. The grains experience large deformation and are subjected to translation and rotation, where the rigid motions are directly measurable quantities. A more detailed description of the method is presented in Appendix A. In the DDEM construction shown in Figure 8, the Cu flyer is discretized with the use of 7200 elements, while the Cu target is composed of 14,400 elements. The thicknesses of the flyer and target are 1.9 mm and 3.9 mm, respectively. The flyer travels with velocity 308 m/s and, after striking the target, shock waves are sent into both of the samples.

Plate Impact Problem
The constitutive model is implemented to a code of deformable discrete element method (DDEM). In DDEM, each discrete element represents a crystallographic grain. The grains experience large deformation and are subjected to translation and rotation, where the rigid motions are directly measurable quantities. A more detailed description of the method is presented in Appendix A. In the DDEM construction shown in Figure 8, the Cu flyer is discretized with the use of 7200 elements, while the Cu target is composed of 14,400 elements. The thicknesses of the flyer and target are 1.9 mm and 3.9 mm, respectively. The flyer travels with velocity 308 m/s and, after striking the target, shock waves are sent into both of the samples. Figure 8. Copper target is struck by a copper flyer with velocity 308 m/s. The entire system is constructed with the use of 7200 triangle particles, where each particle is composed of three elements. Artificial viscosity is not used in this simulation.
A simple Mie Gruneisen equation of state for copper is used. The equation is further modified to account for the contribution of dynamic excitations. It is worth stating that the calculations are made without the use of artificial viscosity. In Figure 9, VISAR measurements [64] are compared with the results of this simulation. Note that the constitutive model correctly predicts the entire impact responses, and this includes the pull-back signals. Four points, A, B, C and D, are selected for plotting contours of damage. The damage is scaled between zero and one, as described in Section 5. The scaling factor refers to the degradation of elastic properties from zero (no damage) to one (open crack). The damage factor 1 − is determined in terms of plasticity-induced dilatation . As discussed earlier, the volume change is triggered by slip incompatibilities in tension. The contact between flyer and target experiences periodic separations. The spallation damage is clearly seen in the middle section of the target. Less severe damage is noticed in other areas of the samples. The damages might be responsible for the observed ringing in the VISAR plot.  A simple Mie Gruneisen equation of state for copper is used. The equation is further modified to account for the contribution of dynamic excitations. It is worth stating that the calculations are made without the use of artificial viscosity. In Figure 9, VISAR measurements [64] are compared with the results of this simulation. Note that the constitutive model correctly predicts the entire impact responses, and this includes the pull-back signals. Four points, A, B, C and D, are selected for plotting contours of damage. The damage is scaled between zero and one, as described in Section 5. The scaling factor refers to the degradation of elastic properties from zero (no damage) to one (open crack). The damage factor 1 − e −H d kk /H 0 kk is determined in terms of plasticity-induced dilatation H d kk . As discussed earlier, the volume change is triggered by slip incompatibilities in tension. The contact between flyer and target experiences periodic separations. The spallation damage is clearly seen in the middle section of the target. Less severe damage is noticed in other areas of the samples. The damages might be responsible for the observed ringing in the VISAR plot.

Plate Impact Problem
The constitutive model is implemented to a code of deformable discrete element method (DDEM). In DDEM, each discrete element represents a crystallographic grain. The grains experience large deformation and are subjected to translation and rotation, where the rigid motions are directly measurable quantities. A more detailed description of the method is presented in Appendix A. In the DDEM construction shown in Figure 8, the Cu flyer is discretized with the use of 7200 elements, while the Cu target is composed of 14,400 elements. The thicknesses of the flyer and target are 1.9 mm and 3.9 mm, respectively. The flyer travels with velocity 308 m/s and, after striking the target, shock waves are sent into both of the samples. Figure 8. Copper target is struck by a copper flyer with velocity 308 m/s. The entire system is constructed with the use of 7200 triangle particles, where each particle is composed of three elements. Artificial viscosity is not used in this simulation.
A simple Mie Gruneisen equation of state for copper is used. The equation is further modified to account for the contribution of dynamic excitations. It is worth stating that the calculations are made without the use of artificial viscosity. In Figure 9, VISAR measurements [64] are compared with the results of this simulation. Note that the constitutive model correctly predicts the entire impact responses, and this includes the pull-back signals. Four points, A, B, C and D, are selected for plotting contours of damage. The damage is scaled between zero and one, as described in Section 5. The scaling factor refers to the degradation of elastic properties from zero (no damage) to one (open crack). The damage factor 1 − is determined in terms of plasticity-induced dilatation . As discussed earlier, the volume change is triggered by slip incompatibilities in tension. The contact between flyer and target experiences periodic separations. The spallation damage is clearly seen in the middle section of the target. Less severe damage is noticed in other areas of the samples. The damages might be responsible for the observed ringing in the VISAR plot.

Conclusions
The paper summarized a five-year effort aimed for the development of a mechanisms-based viscoplasticity concept for metals. The constitutive description comprises of several novel ideas.

1.
The macroscopic plastic flow results from plastic slippages and slip reorganizations. The description is constructed on the basis of the tensor representation concept. It is my conviction that tensor representations derived for generic dyads represent useful tools in the hands of a modeler.

2.
In the proposed model, thermally activated processes are considered stochastic. The concept explains the transition of flow mechanisms from power-law creep to high strain rate dislocation glide. 3.
The proposed description of plasticity-induced heating is based on the hypothesis that plasticity-induced heating quantifies the efficiency of the plastic flow process, while plastic work aids in configurational entropy (suppleness) of the material.

4.
Drag on dislocations is activated by dynamic excitations. As shown, the excitations result from the kinematically-necessary readjustments of flow pathways. 5.
The stress-strain relations are constructed in the framework of transitional viscoplasticity. The power-law relations enable a smooth elastic-plastic transition during loading and unloading processes. 6.
I developed an energy-based Hall-Petch relation, where the commonly known stress-based relation is replaced by its kinematics-based counterpart. The proposed Hall-Petch concept was born out of extensive discussions with Ron Armstrong, who walked me through the sixty years of Hall-Petch interpretations, for which I am grateful. 7.
The model is calibrated for OFHC copper, implemented to a deformable discrete element code (my Ph.D. thesis) and validated in simulations of a plate impact problem. The method itself describes a semi-Cosserat medium, where grain translations and rotations are accounted for.
The most important conclusion of the study is that the flow mechanisms, transitional viscoplasticity, thermal activation, drag on dislocations and size effect are inseparable parts of the constitutive description. The motion of deformable particles includes translations and rotations of mass centers (inodes), while particle-to-particle interactions are defined in e-external nodes. The particle deformation is measured in terms of relative displacements Δ = − ′ between points ′ and ′′ . The relations become Δ = − ( + × ) . In this expression, vectors connect selected external node with the mass center. The deformation matrix links strains and relative displacements = • Δ .
(A1) The rotational deformation is defined with respect to the mass center and is = : Δ , (A2) where = ∑ Δ ⁄ ⁄ . The relative displacements Δ represent rotational deformation subtracted from the rigid core rotations. Each internal node is coupled with −external nodes. The rotation-induced strain is responsible for additional storage of elastic energy During an active deformation process, particles may become separated, and then, the number of interconnected nodes may vary. Forces in mass centers are collected from external nodes. Moments are calculated as described in Equation (A4). In this notation, mass in external nodes is denoted , mass and moment of inertia in internal nodes are and . In DDEM, mass is redistributed between external and internal nodes in such a manner that inertia forces due to rotations and translations are decoupled. Finite elements (particles) can be subjected to large reshaping and repositioning. Thus, the DDEM method is not just a numerical solver, but it carries information on the size and potential shape of grains. The motion of deformable particles includes translations u i and rotations ϕ of mass centers (i-nodes), while particle-to-particle interactions are defined in e-external nodes. The particle deformation is measured in terms of relative displacements ∆u ei = u e − u e between points e and e . The relations become ∆u ei = u e − (u i + r × ϕ). In this expression, vectors r connect selected external node with the mass center. The deformation matrix B links strains H ie and relative displacements The rotational deformation is defined with respect to the mass center and is where H i ϕ = N i ∆u i ϕ /r i /N i . The relative displacements ∆u i ϕ represent rotational deformation subtracted from the rigid core rotations. Each internal node is coupled with N i − external nodes. The rotation-induced strain is responsible for additional storage of elastic energy U ϕ = µ β ϕ H i    During an active deformation process, particles may become separated, and then, the number of interconnected nodes may vary. Forces in mass centers are collected from N ie external nodes. Moments are calculated as described in Equation (A4). In this notation, mass in external nodes is denoted M e , mass and moment of inertia in internal nodes are M i and J M . In DDEM, mass is redistributed between external and internal nodes in such a manner that inertia forces due to rotations and translations are decoupled. Finite elements (particles) can be subjected to large reshaping and repositioning. Thus, the DDEM method is not just a numerical solver, but it carries information on the size and potential shape of grains.