Mechanisms-Based Transitional Viscoplasticity

Abstract

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 unified constitutive model. The model is first tested under simple loading conditions and, later, is validated in a numerical analysis of the plate impact problem, where a copper flyer 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.

1. 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 relevant material information, such as sensitivities to strain rate, temperature, size effect, etc. In a typical scenario, the models adapt the strain rate-insensitive solving scheme, where radial return seems to be the favorite procedure. In consequence, plastic strain is merely a byproduct of the requirement that stress must adhere to the yield surface. An interesting technique suitable for generating random dislocation arrangements is proposed by Gurrutxanga-Lerma [23]. The main concern of the study is focused on how a given population of dislocation structures influences the dominant features of the stress field. The stochastic analysis provides justification for Taylor’s equation and the Hall–Petch relation. The study, however, is not concerned with how the dislocation structures are formed in the first place. Another approach to crystal plasticity is proposed by Langer [24]. His thermodynamics dislocation theory (TDT) is based on the observation that dislocation experience complex chaotic motions and the motion is responsible for the field fluctuations. One should expect the frequency of the fluctuations to be many orders of magnitude lower in comparison with the frequency of thermal excitations. According to Langer, the non-deterministic nature of dislocation dynamics makes the behaviors suitable for statistical and thermodynamics treatment. Brown [25] offered an interesting idea that plastic flow is responsible for the formation of dislocation structures in the state of self-organized criticality. In this interpretation, the dislocation structures evolve until reaching the state of maximum structural adaptiveness such that every site of the structure has a statistically equal chance to act as a nucleation site of further deformation. At the initial stages of plastic deformation, flow is initiated at points of weakness. However, as the process proceeds, the initiation sites become inactive and are replaced by newly created ones in the evolving system. Consequently, there are no fixed sites in time and space. The state of self-organized criticality is dynamic, the fluctuations are inherently erratic and carry a measurable amount of energy, where the energy is fed back into the transient structures. The fluctuations are micro-mechanical and, while allowing to explore the neighboring states, drive the system back into the state of self-organized criticality. The scale-free patterning of self-organized dislocations is further tested in ice [26]. Again, the results confirm the existence of scale-free dislocation assemblies in single crystals. The power-law distribution of dislocation avalanches is insensitive to temperature. Unlike single crystals, polycrystals display a strong deviation from the power-law distribution. Zaiser and Nikitas [27] offered yet another clarification. In their view, the dislocation avalanches exhibit lamellar geometry and, therefore, a very small volume of the material is actively involved in the deformation process. Avalanches produce small macroscopic plastic strain, thus, suggesting an extrinsic limit to the avalanche size. The explanation supports the idea of self-organized criticality.

Here, the proposed concept is aimed to address several questions raised by Brown, Langer, Zaiser and Nikitas. In this study, the mechanism of plastic flow is introduced first, followed by a discussion on the thermal activation of plastic processes. The following section is focused on flow complexities. Slip bursts are shown to activate microscale dynamic excitations. At high strain rates, the excitations are responsible for the dynamic upturn of stress and magnify plasticity-induced heating. Lastly, I construct an energy-based Hall–Petch relation. In contrast to the accepted stress-based interpretation, the proposed Hall–Petch relation is directly linked to the flow constraints. I argue that the stress-based and kinematics-based interpretations represent the upper and lower bound estimates of the grain size effect. For clarity, I made an effort to write each section in a semi-autonomous manner.

All the mechanisms are incorporated into a unified constitutive model. The model reproduces metal responses from quasi-static loading to extreme dynamics and at temperatures from cryogenic environments to melting points. The average size of crystals is shown to affect the responses in a manner that is consistent with experimental observations. The model is calibrated for OFHC copper and has been tested against available experimental data. Lastly, the model is implemented to a code of deformable discrete element method and is validated in a plate impact analysis.

2. Mechanism of Plastic Flow

In metals, the macroscopic plastic flow ${\dot{H}}^p$ results from slip events along many $\theta$ planes defined by two orthogonal unit vectors, $n^{\theta }$ and $s^{\theta }$. 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

$$\left\{\begin{array}{c}{\dot{n}}^{\theta }\\ {\dot{s}}^{\theta }\end{array}\right\}=\left[\begin{array}{cc}0& \dot{\theta }\\ -\dot{\theta }& 0\end{array}\right]\left\{\begin{array}{c}{n}^{\theta }\\ s^{\theta }\end{array}\right\}.$$

(1)

The dominant slip is defined at $\theta =0$. At each angle, slip events occur with frequency $f_{\theta }\left(\theta \right)$, where the angle is taken between $\pm \pi /4$. The distribution $f_{\theta }\left(\theta \right)=\frac{2 \mathsf{\Gamma }\left(1+1/2r\right)}{\sqrt{\pi } \mathsf{\Gamma }\left(1/2+1/2r\right)}{\cos }^{1/r}2\theta$ is constructed in such a manner that ${{\displaystyle \int }}_{-\pi /4}^{\pi /4}{f}_{\theta }\left(\theta \right)d\theta =1$ and $\mathsf{\Gamma }$ 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 $\pm \pi /4$.

The rate of macroscopic plastic strain${\dot{H}}^p$ results not only from local slip events, but also reflects the process of slip reorganizations (Equation (1)), hence

$${\dot{H}}^p=\frac{1}{2}{{\displaystyle \int }}_{-\pi /4}^{\pi /4}\frac{\partial \left[\left(n^{\theta }{s}^{\theta }+s^{\theta }{n}^{\theta }\right) e^{\theta } \right]}{\partial t} d\theta +{\dot{H}}^d.$$

(2)

The volumetric strain rate ${\dot{H}}^d$ quantifies slip incompatibilities in tension. Plastic strain $e^{\theta }$ is assigned to each $\theta$-slip plane. The $\theta$-strains $e^{\theta }$ are weighted ${\dot{e}}^{\theta }=f_{\theta }\left(\theta \right) {\dot{e}}^p$, where ${\dot{e}}^p$ is the rate of macroscopic strain. Since the reorientations are local events, therefore $\dot{\theta }=A_T\left(T\right) \left(\partial \theta /\partial e^{\theta }\right) {\dot{e}}^{\theta }$. 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 $\left(\theta =0\right)$ do not rotate [28]. The slip reorganizations are maximized at $\theta =\pm \pi /4$. The relation $e^{\theta } \left(\partial \theta /\partial e^{\theta }\right)=2\sin 2\theta$ properly reflects the trend. The dominant slip plane is defined in terms of two orthogonal unit vectors $n^{\sigma }$ and $s^{\sigma }$ and the flow tensor is a dyadic product of the two $N=\left(n^{\sigma }{s}^{\sigma }+s^{\sigma }{n}^{\sigma }\right)$. Since plastic flow is stimulated by stress $\sigma$, therefore, the tensor must be expressed in terms of the current stress $N=N^{\sigma }\left(\sigma \right)$. The representation $N^{\sigma }={\alpha }_0 \mathbf{1}+{\alpha }_1 S+{\alpha }_2 S^2$ of the dyadic product is constructed on the basis of stress deviator $S=\sigma -\mathbf{1} tr\sigma /3$ and the three variables ${\alpha }_0$, ${\alpha }_1$ and ${\alpha }_2$ are functions of stress invariants. The representation must reproduce invariants of the original flow tensor $\left(n^{\sigma }{s}^{\sigma }+s^{\sigma }{n}^{\sigma }\right)$, thus $tr{N}^{\sigma }=0$, $tr{\left(N^{\sigma }\right)}^2=2$ and $N^{\sigma }={\left(N^{\sigma }\right)}^3$. There are three solutions, where the representations are constructed on the planes of principal stresses $\left\{{\sigma }_1, {\sigma }_2\right\}$, $\left\{{\sigma }_2, {\sigma }_3\right\}$ and $\left\{{\sigma }_1, {\sigma }_3\right\}$ and where ${\sigma }_1>{\sigma }_2>{\sigma }_3$. I select the solution which reproduces the maximum shear mechanism [29], hence

$$N^{\sigma }=\frac{2cos\frac{\phi }{3}}{\sqrt{3} cos\phi } \left(\mathbf{1}-\frac{3}{2 J_2} S^2\right)+\frac{cos\frac{2\phi }{3}}{\sqrt{J_2} cos\phi } S.$$

(3)

The angle $\phi ={\sin }^{-1}\left(J_3\sqrt{27/4}/J_2^3 \right)$ and two invariants $J_2$ and $J_3$ of stress deviator $S$ complete the relation. It turns out that$N^{\sigma }:\sigma /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 $\theta$-vectors$n^{\theta }$ and$s^{\theta }$ are co-rotational with the dominant slip plane$n^{\theta }=n^{\sigma }cos\theta +s^{\sigma }sin\theta$ and$s^{\theta }=s^{\sigma }cos\theta -n^{\sigma }sin\theta$. As a result, the plastic flow (Equation (2)) becomes

$${\dot{H}}^p=\frac{1}{2} N^{\sigma } M{\dot{e}}^p+{\dot{H}}^d,$$

(4)

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=\frac{{\mathsf{\Gamma }}^2\left(1+\frac{1}{2r}\right)}{\mathsf{\Gamma }\left(\frac{1}{2}+\frac{1}{2r}\right) \mathsf{\Gamma }\left[\frac{1}{2}\left(3+\frac{1}{r}\right)\right]}$ is closely approximated by$M_s\cong \frac{1+r}{1+\pi r/2}$, where the error of approximation does not exceed one percent. The second term$M_r=\frac{4 A_T\left(T\right)}{\pi \left(1+2r\right)}$ 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\left(e^{pe}\right)$. The exponent$r=r_0\left(1+e^{pe}/e_r^p\right)$ reflects the slip organization, where$r_0$ and$e_r^p$ 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${\dot{H}}^d=\frac{q}{2}{\left(N^{\sigma }\right)}^2 M{\dot{e}}^{pe}$. The rate of damage is controlled by the parameter$q$.

3. 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${\xi }_a=F_a/k_{b}T$, where$k_b$ is Boltzmann constant and$T$ is temperature. Thermal activation is stochastic with respect to${\xi }_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 ${\xi }_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 {\xi }_a^{k_a-1} exp\left[-{\xi }_a^{k_a}\right]$ seems to fit the description. The material is expected to increase the number of activation sites with frequency$f_T\left(T\right)$ 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={{\displaystyle \int }}_{T_c/T_m}^{T_c/T}{f}_T\left(s\right)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_a^{1/k_a} k_b{T}_c$. Activation entropy is estimated at melting point$F_a\left(T_m\right)=0$ and we obtain$S_a=E_a/T_m$. The activation factor becomes

$$A_T\left(T\right)=1-exp\left[-g_a{\left(\frac{T_c}{T}-\frac{T_c}{T_m}\right)}^{k_a}\right],$$

(5)

where$g_a$ determines the magnitude of activation energy. At melting temperature, when$A_T\left(T_m\right)\to 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.

4. 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).

Initially, the particle resides in position$\left\{X\right\}$ and stress is${\sigma }^X$. Shortly later$\left(\delta t\right)$, the particle is forced to change its trajectory and arrives at position$\left\{x\right\}$, where stress is$\sigma$. Should the original path be available, the particle would move to position$\left\{z\right\}$ with velocity${\upsilon }^z$ and the stress would be${\sigma }^z$. In all the scenarios, equations of motion$\mathbf{\nabla }·\sigma =\rho \dot{\upsilon }$ must be preserved. The particle acceleration and mass density are denoted by$\dot{\upsilon }$ and$\rho$. The rerouting alters tractions projected onto the surface$\partial V_0$ that is normal to the direction of the particle velocity$m$. The difference in tractions$\left({\sigma }^z·m\ne \sigma ·m\right)$ between positions $\left\{z\right\}$ and $\left\{x\right\}$ triggers stress perturbations$\delta \sigma =\frac{l_r}{V_0} {\int }_{\partial V_0}\left(\sigma -{\sigma }^z\right)·\left(m m\right)dS$. The perturbations are stretched over distance$l_r$. Next, the volume$V_0$ is reduced to a material point and stress perturbations become$\delta \sigma =l_r \dot{\psi }$, where$\psi =\rho \left(\delta \upsilon m+m \delta \upsilon \right)/2$ is the momentum tensor. The path rearrangements affect particle acceleration$\delta \dot{\upsilon }=\left(\dot{\upsilon }-{\dot{\upsilon }}^z\right)$, while the overstress$\left({\sigma }^z-\sigma \right)$ is partly stored in dislocation structures (Figure 1). The uncompensated stress$\delta S^a$ in Figure 1 activates dynamic excitations. Given sufficient time, the overstress is fully relaxed.

5. Dynamic Overstress

The stress perturbations $\delta \sigma$ slow down plastic flow${\dot{H}}^t={\dot{H}}^e+\left({\dot{H}}^p-\delta \sigma /\kappa \right)$, where${\dot{H}}^t$, ${\dot{H}}^e$ and${\dot{H}}^p$ are the rates of total, elastic and drag-free plastic strains. The term $\delta \sigma /\kappa$ represents the drag on dislocations and the drag coefficient is$\kappa$. Stress is partly stored in the lattice$\dot{\sigma }=C·{\dot{H}}^e$, where$C$ is the elastic matrix. During an active plastic process, stress$C·{\dot{H}}^t$ exceeds the elastic limit and its uncompensated part triggers dynamic excitations$\dot{\sigma }+\kappa \dot{\psi }/\rho l_r=C·{\dot{H}}^t$. The excitations give rise to the viscous overstress $\kappa \dot{\psi }/\rho l_r$. The overstress can be expressed in terms of stress perturbations $\kappa \dot{\psi }/\rho l_r=R_k^{-2}\delta \sigma /\tau$, where the resistance to plastic flow$R_k=l_r/\tau {\upsilon }_s$ carries information on the size of the perturbed domain$l_r$ and specifies the relaxation time $\tau =\kappa /\mu$ required for the equilibration of the dislocation structures. Here,${\upsilon }_s=\sqrt{\mu /\rho }$ is shear velocity and $\mu$ is shear modulus. As a result, we have$\dot{\sigma }+R_k^{-2}\delta \sigma /\tau =\mathit{C}·{\dot{H}}^t$. Once again, the resistance to plastic flow$R_k={\upsilon }_r/{\upsilon }_s$ arises in the actively perturbed domain$l_r$, such that${\upsilon }_r=l_r/\tau$. The excitations are redistributed into the surroundings with velocity slower than shear velocity${\upsilon }_s$, thus the resistance$R_k$ cannot be greater than one. In summary, the material description is given in terms of three relations

$$\begin{array}{c}{\dot{\mathit{H}}}^t={\dot{\mathit{H}}}^e+\left({\dot{\mathit{H}}}^p-\delta \sigma /\kappa \right)\\ \dot{\sigma }+R_k^{-2}\delta \sigma /\tau =C·{\dot{\mathit{H}}}^t\\ \dot{\sigma }=C·{\dot{\mathit{H}}}^e\end{array}.$$

(6)

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

$$\begin{array}{c}{\dot{\mathit{H}}}^t={\dot{\mathit{H}}}^e+{\dot{\mathit{H}}}^{pe}\\ \dot{\sigma }+C·{\dot{\mathit{H}}}^{pe}=C·{\dot{\mathit{H}}}^t\\ \delta \sigma =\tau R_k^2 C·{\dot{\mathit{H}}}^{ep}\end{array}.$$

(7)

The rate of effective (true) plastic strain${\dot{H}}^{pe}=P^{-1}·{\dot{H}}^p$ is scaled by the fourth-order drag tensor$P=I+R_k^2 C/\mu$. 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/\left(1+2 R_k^2\right)$ and the rate of plastic strain${\dot{H}}^{pe}={\dot{H}}^p/\left(1+2R_k^2\right)$ take much simpler forms. Comparison of Equations (6₂) and (7₂) $\left(C·{\dot{H}}^{pe}=R_k^{-2}\delta \sigma /\tau \right)$ 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${\dot{S}}^a+S^a/2 \tau =C·{\dot{H}}^{pe}$, where$S^a$ is the active overstress (Figure 2). Rerouting of plastic flow is enabled by the excess of energy ${\dot{W}}_l=S^a:{\dot{S}}^a/2\mu$, where$tr{S}^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_k^0$ 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$R_k=R_k^0+\left(1-R_k^0\right)\left(1-e^{-{\left(W_l/G_l\right)}^{n_r}}\right)$ 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_k^0$ 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.

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${\left(N^{\sigma }\right)}^2=\left(n^{\sigma } n^{\sigma }+s^{\sigma } s^{\sigma }\right)$ and is${\dot{H}}^d=\frac{q}{2}{\left(N^{\sigma }\right)}^2 M{\dot{e}}^{pe}$. Subsequently, the volumetric change is$tr{\dot{H}}^d=qM{\dot{e}}^{pe}$. 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_d^{pe}/e_d}$ are controlled by ductility$e_d$, where plastic strain$e_d^{pe}$ is accumulated during tensile loading only. The damage is responsible for the degradation of the apparent bulk and shear moduli. The scaling parameter is$\left(1-e^{-H_{kk}^d/H_{kk}^0}\right)$. Here,$H_{kk}^0$ specifies the strain at which the damage process is macroscopically evident. Consequently, in a damage-free metal, the scaling factor $\left(1-e^{-H_{kk}^d/H_{kk}^0}\right)$ is equal to zero. The factor is equal to one in a completely fractured material.

6. 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$\dot{Q}=\left(\beta {\sigma }_{eq}\right) {\dot{e}}^{pe}$, where the Taylor–Quinney coefficient$\beta$ 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$\delta \sigma =l_r \dot{\psi }$ are expressed in terms of the plastic strain rate$\delta \sigma =\tau R_k^2 (C·{\dot{H}}^{pe})$. Thus, the change in momentum is$\dot{\psi }=R_k C·{\dot{H}}^{pe}/{\upsilon }_s$. In Maxwell’s process$({\dot{S}}^a+S^a/2 \tau =C·{\dot{H}}^{pe})$, the part$\left(S^a/2 \tau \right)$ in$C·{\dot{H}}^{pe}$ quantifies the already dissipated excitations

$${\dot{\psi }}_d=R_k S^a/2{\upsilon }_s\tau .$$

(8)

The energy of the excitations is absorbed in dislocation structures. Still, there a fraction${\xi }_{ef}$ of${\dot{\psi }}_d$ which gives rise to phonon vibrations${\dot{\mathit{\psi }}}_Q={\xi }_{ef}{\dot{\mathit{\psi }}}_d$. The term${\mathit{S}}^a/2\tau$ in Equation (8) is replaced by $(C·{\dot{H}}^{pe}-{\dot{S}}^a)$, where the phonon-scale overstress${\dot{S}}^a$ is a vanishing quantity. Therefore, the heat generating excitations${\dot{\psi }}_Q\cong {\xi }_{ef} R_k(C·{\dot{H}}^{pe})/{\upsilon }_s$ are activated by the overstress$C·{\dot{H}}^{pe}$ and are affected by the flow constraints$R_k$. As plastic deformation advances, heat sources continue to evolve${\psi }_Q={\psi }_p+{\xi }_{ef} R_{k}C·H^{pe}/{\upsilon }_s$. In this expression, momentum${\psi }_p$ is linked to shear stress capable of triggering path rearrangements.

Heat is produced when excitations $\dot{\psi }$ act on the existing sources of heat ${\psi }_Q$, such that$\dot{Q}={\psi }_Q:\dot{\psi }/\rho$. Consequently, plasticity-induced heat is generated with rate

$$\dot{Q}=2R_k\left[R_k {\xi }_{ef} C·H^{pe}+{\upsilon }_s {\psi }_p\right]:{\dot{H}}^{pe}.$$

(9)

The relation is further simplified$\dot{Q}=2R_k\left[R_k {\xi }_{ef}\mu e^{ep}+{\sigma }_y\right] {\dot{e}}^{pe}$, where equivalent plastic strain$e^{pe}$ has already been defined and${\sigma }_y$ is yield stress. The Taylor–Quinney coefficient $\beta =\dot{Q}/{\dot{W}}^p$ becomes$\beta =2R_k\left[R_k{\xi }_{ef}\mu e^{pe}+{\sigma }_y\right]/{\sigma }_{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${\xi }_{ef}$ 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.

7. Hall–Petch Relation

In early 1950s, Hall [10] and Petch [11] introduced a relation, which correlates yield stress${\sigma }_y$ with an average size of grains$d$

$${\sigma }_y={\sigma }_y^p+k_y d^{-1/2}.$$

(10)

In this construction, ${\sigma }_y^p$ 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${\sigma }_y={\sigma }_{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$\left(\mu b^{1/2}\right)$. 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 ${\dot{W}}_{{\dot{e}}^p}^p=\left(1+2R_k^2\right) {\sigma }_y^p M {\dot{e}}^p$. Here, the constraints $R_k$ magnify the effective yield stress${\sigma }_y={\sigma }_y^p \left(1+2R_k^2\right)$. When assuming that the Hall–Petch relation prevails, one can obtain$R_k^2=k_y d^{-1/2}/2{\sigma }_y^p$. Next, I construct a kinematics-based interpretation. At the prescribed stress${\sigma }_y$, the resistance$R_k$ is responsible for slowing down the plastic flow. For this reason, we write ${\dot{W}}_{{\sigma }_y}^p={\sigma }_y M {\dot{e}}^p/\left(1+2R_k^2\right)$. The two interpretations represent the upper and lower bound estimates of the size effect.

7.1. 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${\sigma }_y^p$. 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${\sigma }_y^p$, the barriers are partly climbed$g_{gb}^y$. 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}=\left(V_s/S_{gb}\right){{\displaystyle \int }}_{{\sigma }_y^p}^{{\sigma }_y}w\left({\sigma }_s\right) {\sigma }_{s}d{\sigma }_s/2\mu$. The average shear stress${\sigma }_s$ continues to rise and becomes the macroscopic yield stress${\sigma }_y$. At this point, dislocations cross the boundaries in large numbers. The average rate of energy accumulation slows down. The correction function$w=\left({\sigma }_y-{\sigma }_y^p\right)/2{\sigma }_s$ captures the slowing down effect. When energy$E_{gr}$ overcomes the energy barrier$(g_{gb}-g_{gb}^a)\left(S_{gb}/V_s\right)$, the macroscopic plastic flow is initiated

$$(g_{gb}-g_{gb}^a) \left(S_{gb}/V_s\right)=\frac{{\left({\sigma }_y-{\sigma }_y^p\right)}^2}{4 \mu }.$$

(11)

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

$${\sigma }_y={\sigma }_y^p+k_{yS} {\left(\frac{S_{gb}}{V_s}\right)}^{1/2}.$$

(12)

In this expression, the Hall–Petch mechanism is controlled by grain boundary area per unit volume. The parameter$k_{yS}=2\sqrt{\mu \left(g_{gb}-g_{gb}^a\right)}$ is closely related to constant$k_y$ in Equation (10), such that$k_{yS}\cong 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_{gb}^a)=\mu ({\alpha }_{yS}^2 b)$ suggests that$k_{yS}=2 \mu {\alpha }_{yS} b^{1/2}$ and ${\alpha }_{yS}$ carries information on the metal crystallographic structure. It turns out that the parameter${\alpha }_{yS}$ is a material constant in fcc and bcc metals. Based on the data collected in [52], the mean value of ${\alpha }_{yS}$ for fcc metals is${\alpha }_{yS}=0.04604\pm 0.00616$. In bcc metals, the values are found to be${\alpha }_{yS}=0.1089\pm 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.

7.2. Kinematics-Based Construction of Hall–Petch Relation

At constant stress, energy barriers slow down plastic flow and, in this manner, affect plastic power ${\dot{W}}_{{\sigma }_y}^p={\sigma }_y M {\dot{e}}^p/\left(1+2R_k^2\right)$. In the energy-based construction, the term$\left(1+2R_k^2\right)$ in${\dot{\mathit{H}}}^{pe}={\dot{\mathit{H}}}^p/\left(1+2R_k^2\right)$ becomes $\left[1+\frac{2 \mu {\alpha }_{yS} b^{1/2}}{{\sigma }_y^p} {\left(S_{gb}/V_s\right)}^{1/2}\right]$. The pre-existing resistance to plastic flow$R_k=R_k^0$ becomes a quantifiable material constant

$$R_k=\sqrt{\frac{k_{yS}}{2 {\sigma }_y^p}} {\left(S_{gb}/V_s\right)}^{1/4}.$$

(13)

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_k^0$ is approaching one, thus the size effect becomes less pronounced. As shown in Section 5, the resistance$R_k=l_r/\left({\upsilon }_s \tau \right)$ is a function of relaxation time$\tau =l_r/\left(R_k{\upsilon }_s\right)$. The length$l_r$ is directly related to the grain boundary area per unit volume and is$l_r={\left(S_{gb}/V_s\right)}^{-1}$. Based on the above, the relaxation time becomes the size-dependent material property

$$\tau =\frac{{\left(S_{gb}/V_s\right)}^{-5/4}}{{\upsilon }_s} \sqrt{\frac{2 {\sigma }_y^p}{k_{yS}}}.$$

(14)

It is worth stating that relaxation time$\tau$ (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.

8. Transitional Viscoplasticity

Plastic power$\mathit{\sigma }:{\dot{H}}^{pe}={\sigma }_{eq}{\dot{e}}^{pe}$ is expressed in terms of equivalent stress${\sigma }_{eq}=\frac{1}{2}{\mathit{N}}^{\mathit{\sigma }}:\mathit{\sigma }$ and strain rate ${\dot{e}}^{pe}$. During active deformation processes, plastic power${\dot{W}}^p$ aids to configurational entropy (suppleness [25]) of the material. The rate of true plastic strain${\dot{\mathit{H}}}^{pe}=\frac{1}{2}\left[{\mathit{N}}^{\mathit{\sigma }}+q{\left({\mathit{N}}^{\mathit{\sigma }}\right)}^2\right] {\dot{e}}^{pe}$ is defined in terms of the effective rate ${\dot{e}}^{pe}=\frac{M {\dot{e}}^p}{1+2R_k^2}$, where the resistance-free strain rate${\dot{e}}^p$ is a function of equivalent stress and the relation is

$${\dot{e}}^p={\dot{\mathsf{\Lambda }}}_p·{\left({\sigma }_{eq}/{\sigma }_y^p\right)}^{n_p}.$$

(15)

The power-law viscoplasticity was first introduced in [58,59]. Strength${\sigma }_y^p={\sigma }_0 A_T\left(T\right)$ is determined at 0.2% strain, where${\sigma }_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${\dot{\mathsf{\Lambda }}}_p={\mathsf{\Lambda }}_0 {\dot{e}}_N^0{\left[{\dot{e}}_N^t/{\dot{e}}_N^0\right]}^{{\omega }_p}$, where ${\mathsf{\Lambda }}_0$ and ${\omega }_p$ are constants and${\dot{e}}_N^0=1/s$. The normalized rate of total strain${\dot{e}}_N^t=\sqrt{{\dot{\mathit{H}}}^t:{\dot{\mathit{H}}}^t/2}$ completes the expression. As shown earlier, the average Schmid factor$M=\frac{1+r}{1+\pi r/2}-\frac{4 A_T}{\pi \left(1+2 r\right)}$ reflects slip plane misorientations, where the exponent$r$ controls the evolution of the process (Figure 1). Advances of plastic deformation broaden the spectrum of active planes and, therefore,$r=r\left(e^{pe}\right)$.

At the initial stage of deformation, the rate of elastic strain dominates$({\dot{\mathit{H}}}^e\to {\dot{\mathit{H}}}^t)$ and, therefore, we have${\dot{e}}^p\propto {\left({\sigma }_{eq}/{\sigma }_y^p\right)}^{n_p}$. At advanced deformation, plastic flow takes over$({\dot{\mathit{H}}}^p\to {\dot{\mathit{H}}}^t)$ and the mechanism evolves${\dot{e}}^p\propto {\left({\sigma }_{eq}/{\sigma }_y^p\right)}^{n_p/\left(1-{\omega }_p\right)}$, accordingly. Thus, there is a smooth transition from diffusional flow ${\dot{e}}^p\propto {\left({\sigma }_{eq}/{\sigma }_y^p\right)}^{n_p}$ to dislocation glide${\dot{e}}^p\propto {\left({\sigma }_{eq}/{\sigma }_y^p\right)}^{n_p/\left(1-{\omega }_p\right)}$, where we should note that$n_p/\left(1-{\omega }_p\right)\gg n_p$.

9. 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. Text-book material constants.

Bulk Modulus	Shear Modulus	Mass Density	Yield Stress, 298 K	Burgers Vector	Melting Point	Specific Heat, 298 K
$B$	$\mu$	$\rho$	${\sigma }_{y0.2\%}$	$b$	$T_m$	$C_p$
$GPa$	$GPa$	$kg/m^3$	$MPa$	$nm$	$K$	$J/\left(kg K\right)$
$138$	$56$	$8930$	$84$	$0.2555$	$1356$	385

. 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${\omega }_p$, stress exponent$n_p$, thermal efficiency of plastic flow${\xi }_{ef}$, crystallographic constant${\alpha }_{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. Material constants with marginal variability.

Strain Rate Exponent	Stress Exponent	Heat Coefficient	Crystallographic Constant	Schmid Factor	Transition Temperature	Overstress Exponent
${\omega }_p$	$n_p$	${\xi }_{ef}$	${\alpha }_{yS}$	$r_0$	$T_c$	$n_r$
$-$	$-$	$-$	$-$	$-$	$K$	$-$
$0.99$	$0.56$	$0.031$	$0.044$	$0.3$	$596.6$	0.5

. In the absence of experimental data, the constants should properly approximate responses of fcc metals.

The material-specific parameters are listed in

Table 3. Material-specific parameters.

Activation Energy Factor	Thermal Activation	Stress Pre-Factor	Critical Energy	Schmid Factor	Ductility	Damage Strain
$g_a$	$k_a$	${\mathsf{\Lambda }}_0$	${\mathrm{G}}_l$	$e_r^p$	$e_d$	$H_{kk}^0$
$-$	$-$	$-$	$MJ/m^3$	$-$	$-$	$-$
$1.13$	$0.087$	$6.8$	$5.0$	$0.016$	$0.148$	$0.0296$

. Among them are the activation energy constants $\left(g_a, k_a\right)$, the strength pre-factor ${\mathsf{\Lambda }}_0$, critical energy$G_l$, hardening strain$e_r^p$ and damage parameters$\left(e_d,H_{kk}^0\right)$. 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 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.

10. 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.

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$\left(1-e^{-H_{kk}^d/H_{kk}^0}\right)$ is determined in terms of plasticity-induced dilatation$H_{kk}^d$. 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.

11. 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.

• 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.
• 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.
• 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.
• Drag on dislocations is activated by dynamic excitations. As shown, the excitations result from the kinematically-necessary readjustments of flow pathways.
• 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.
• 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.
• 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.