{'fraction_non_alphanumeric': 0.03782190076933708, 'fraction_numerical': 0.013451881881420739, 'mean_word_length': 4.428357089272318, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Thermodynamic dislocation theory: torsion of bars

The thermodynamic dislocation theory developed for non-uniform plastic deformations is used here in an analysis of a bar subjected to torsion. Employing a small set of physics-based parameters, which we expect to be approximately independent of strain rate and temperature, we are able to simulate the torque-twist curve for a bar made of single crystal copper that agrees with the experimental one.


I. INTRODUCTION
The thermodynamic dislocation theory (TDT), proposed initially by Langer, Bouchbinder, and Lookman [1] and developed further in [2][3][4][5][6][7], deals with the uniform plastic deformations of crystals driven by a constant strain rate. During these uniform plastic deformations the crystal may have only redundant dislocations whose resultant Burgers vector vanishes. As shown in [8][9][10], the extension of TDT to non-uniform plastic deformations should account for excess dislocations due to the incompatibility of the plastic distortion [11]. There are various examples of non-uniform plastic deformations in material science and engineering, the most typical of which being the torsion of bars [12] and the bending of beams [13]. The purpose of this paper is to explore use of TDT for non-uniform plastic deformations [8][9][10] in modeling bars made of single crystal copper and subjected to torsion. Our challenge is to simulate the torque-twist curve exhibiting the hardening behavior and the size effect. We also want to compare this torque-twist curve with the experimental curve provided by Horstemeyer et al. [14]. To make this comparison possible we will need to identify from the experimental data obtained in [14] a list of material parameters for single crystal copper under torsion. For this purpose, we will use the large scale least-squares analysis described in [6,7,15].
The thermodynamic dislocation theory is based on two unconventional ideas. The first of these is that, under nonequilibrium conditions, the atomically slow configurational degrees of freedom of dislocated crystals are characterized by an effective disorder temperature that differs from the ordinary kinetic-vibrational temperature. Both of these temperatures are thermodynamically well defined variables whose equations of motion determine the irreversible behaviors of these systems. The second principal idea is that entanglement of dislocations is the overwhelmingly dominant cause of resistance to deformation in crystals. These two ideas have led to successfully predictive theories of strain hardening [1,2], steady-state stresses over exceedingly wide ranges of strain rates [1], thermal softening during deformation [6], yielding transitions between elastic and plastic responses [3,5], shear banding instabilities [4,7], and size and Bauschinger ef- We start in Sec. II with a brief annotated summary of the equations of motion to be used here. Our focus is on the physical significance of the various parameters that occur in them. We discuss which of these parameters are expected to be material-specific constants, independent of temperature and strain rate, and thus to be key ingredients of the theory. In Sec. III we discretize the obtained system of governing equations and develop the numerical method for its solution. The parameter identification based on the large scale least squares analysis and the results of the numerical simulations are presented in Sec. IV. We conclude in Sec.V with some remarks about the significance of these calculations.

II. EQUATIONS OF MOTION
Suppose a single crystal bar with a circular cross section, of radius R and length L, is subjected to torsion (see drawing of the bar with its cross-section in Fig. 1). For this particular geometry of the bar and under the condition R L it is natural to assume that the warping of the bar vanishes, while the circumferential displacement is u ϕ = ωrz, with ω being the twist angle per unit length. Thus, the total shear strain of the bar γ = 2 ϕz = ωr and the shear strain rateγ =ωr turn out to be non-uniform as they are linear functions of radius r. Now, let this system be driven at a constant twist rateω ≡ 0 /t 0 , where t 0 is a characteristic microscopic time scale. Since the system experiences a steady state torsional deformation, we can replace the time t by the total twist angle (per unit length) ω so that t 0 ∂/∂t → 0 ∂/∂ω. The equation of motion for the flow stress becomes with µ being the shear modulus. This equation is derived from Eq. (II.1) in [10] by replacing γ = rω and multiplying both sides by r. Note that for uniform plastic deformations involving only redundant dislocations q(γ)/t 0 equals the plastic shear rateβ, with β being the uniform plastic distortion. However, if β is non-uniform, it is not necessarily so. The state variables that describe this system are the elastic strain γ − β, the areal densities of redundant dislocations ρ r and excess dislocations ρ g ≡ |β ,r + β/r|/b (where b is the length of the Burgers vector), and the effective disorder temperature χ (cf. [3,16]). All four quantities, γ − β, ρ r , ρ g , and χ, are functions of r and ω.
The central, dislocation-specific ingredient of this analysis is the thermally activated depinning formula for q as a function of a flow stress τ Y and a total dislocation density ρ = ρ r + ρ g : This is an Orowan relation of the form q = ρ b v t 0 in which the speed of the dislocations v is given by the distance between them multiplied by the rate at which they are depinned from each other. That rate is approximated here by the activation terms f P (τ Y , ρ) and −f P (−τ Y , ρ), in which the energy barrier e P = k B T P (implicit in the scaling of θ = T /T P ) is reduced by the stress dependent factor e −τ Y /τ T (ρ) , where τ T (ρ) = µ T b √ ρ is the Taylor stress with µ T being proportional to µ (see Section III). Note that antisymmetry is required in Eq. (2.2), especially when dealing with the load reversal, both to preserve reflection symmetry, and to satisfy the second-law requirement that the energy dissipation rate, τ Y q/( 0 r), is non-negative. The pinning energy e P is large, of the order of electron volts, so that θ is very small. As a result, q(τ Y , ρ) is an extremely rapidly varying function of τ Y and θ. This strongly nonlinear behavior is the key to understanding yielding transitions and shear banding as well as many other important features of crystal plasticity. For example, the extremely slow variation of the steady-state flow stress as a function of strain rate discussed in [1] is the converse of the extremely rapid variation of q as a function of τ Y in Eq.(2.2).
The equation of motion for the total dislocation density ρ = ρ r + ρ g describes energy flow. It says that some fraction of the power delivered to the system by external driving is converted into the energy of dislocations, and that that energy is dissipated according to a detailedbalance analysis involving the effective temperature χ.
In terms of the twist angle ω this equation reads: with ρ ss (χ) = (1/a 2 )e −e D /χ being the steady-state value of ρ at given χ, e D a characteristic formation energy for dislocations, and a denoting the average spacing between dislocations in the limit of infinite χ (a is a length of the order of tens of atomic spacings). The coefficient K ρ is an energy conversion factor that, according to arguments presented in [1] and [4], should be independent of both strain rate and temperature. The other quantity that appears in the prefactor in Eq.
The equation of motion for the effective temperature χ is a statement of the first law of thermodynamics for the configurational subsystem: Here, χ 0 is the steady-state value of χ for strain rates appreciably smaller than inverse atomic relaxation times, i.e. much smaller than t −1 0 . The dimensionless factor K is inversely proportional to the effective specific heat c eff . Unlike K ρ , there is no reason to believe that K is a rateindependent constant. In [5], K for copper was found to decrease from 17 to 12 when the strain rate increased by a factor of 10 6 . Since the maximum strain rate (reached at the outer radius of the bar) for the small twist rate in our torsion test is small, we assume that K is a constant.
The equation for the plastic distortion β reads This equation is the balance of microforces acting on excess dislocations. Here, the first term τ = µ(γ − β) = µ(ωr − β) is the applied shear stress, the second term the back-stress due to the interaction of excess dislocations, and the last one the flow stress. This balance of microforces can be derived from the variational equation for irreversible processes [8,9] yielding with ψ m being the free energy density of excess dislocations. Note that the applied shear stress is equal to the flow stress for the uniform plastic deformations.
Berdichevsky [17] has found ψ m for the locally periodic arrangement of excess screw dislocations in a bar under torsion. However, as shown by us in [9], his expression must be extrapolated to the extremely small or large dislocation densities to guarantee the existence of solution within TDT. Using the extrapolated energy proposed in [9] we find that τ B is given by The second condition means that the density of excess dislocations must vanish at the free boundary.

III. DISCRETIZATION AND METHOD OF SOLUTION
For the purpose of numerical integration of the system of equations (2.1)-(2.8) let us introduce the following variables and quantities The variabler changes from zero toR = R/b. The variable φ has the meaning of the total twist angle measured in degree (in [14] φ changes from zero to φ * = 73.35 • ). The calculation of the torque as function of φ is convenient for the later comparison with the torque-twist curve from [14]. Then we rewrite Eq. (2.2) in the form We setμ T = (b/a)µ T = µs and assume that s is independent of temperature and strain rate. Theñ We The dimensionless steady-state quantities arẽ Usingq instead of q as the dimensionless measure of plastic strain rate means that we are effectively rescaling t 0 by a factor b/a. For purposes of this analysis, we assume that (a/b)t 0 = 10 −12 s. In terms of the introduced quantities the governing equations read (3.11) with ξ = |β ,r + β/r|. To solve this system of partial differential equations subject to initial and boundary conditions numerically, we discretize the equations in the interval (0 <r <R) by dividing it into n sub-intervals of equal length ∆r =R/n. The first and second spatial where β i = β(r i ). For the end-pointr =R we introduce β n+1 at a fictitious pointr n+1 = (n + 1)∆r and find it from the discretized condition of vanishing density of excess dislocations β n+1 − β n ∆r + β n /R = 0. (3.14) Then it is possible again to discretize the first and second derivative of β atr =R and write the finite difference equation for β at that point. In this way, we reduce the four partial differential equations to a system of 4n ordinary differential-algebraic equations that will be solved by Matlab-ode15s.
After finding the solution we can compute the torque as function of the twist angle according to (3.15)

IV. PARAMETER IDENTIFICATION AND NUMERICAL SIMULATIONS
The experimental torque-twist curve for the single crystal copper bar provided in [14] along with our theoretical results based on the preceding equations of motion are shown in Fig. 2. In this figure, the circles represent the experimental data for sample 1 in [14] while the solid curve is our theoretical simulation. The experimental data for sample 2 in that paper appear less reliable, especially at large twist angles, and are not analyzed here. In order to compute the theoretical torque-twist curve, we need values for seven system-specific parameters and two initial conditions. The seven basic parameters are the following: the activation temperature T P , the stress ratio s, the steady-state scaled effective temperatureχ 0 , the two dimensionless conversion factors K ρ and K, the two coefficients k 0 , and k 1 defining the functionτ B in Eq. (3.11). We also need initial values of the scaled dislocation densityρ i and the effective disorder temperatureχ i ; all of which are determined by the sample preparation. The other parameters required for numerical simulations but known from the experiment are: the ambient temperature T = 298K, the shear modulus µ = 48GPa, the length L = 17.6mm and radius R = 6.35mm of the bar, the length of Burgers' vector b = 2.55Å, the twist rateφ = 0.25 • /s, and consequently, 0 = 0.2479 × 10 −12 /m. We take a = 10b. In earlier papers dealing with the uniform deformations [1][2][3][4], it was possible to begin evaluating the parameters by observing steady-state stresses σ ss at just a few strain rates q 0 and ambient temperatures T 0 = T Pθ0 . Knowing σ ss , T 0 and q 0 for three stress-strain curves, one could solve equation which is the inverse of Eq. (2.2) for T P , s, andχ 0 , and check for consistency by looking at other steady-state situations. With that information, it was relatively easy to evaluate K ρ and K by directly fitting the full stressstrain curves. This strategy does not work here because the stress state of twisted bars is non-uniform. We may still have local steady-state stresses as function of the radius r, but it is impossible to extract this information from the experimental torque-twist curve. Furthermore, the similar parameters for copper found in [1][2][3][4] cannot be used here, since we are dealing with screw dislocations having the energy barrier T P and other characteristics different from those identified in the above references.
To counter these difficulties, we have resorted to the large-scale least-squares analyses that we have used in [6,7,15]. That is, we have solved the system of ordinary differential-algebraic equations (DAE) numerically, provided a set of material parameters is known. Based on this numerical solution we then computed the sum of the squares of the differences between our theoretical torque-twist curve and a large set of selected experimental points, and minimized this sum in the space of the unknown parameters. The DAE were solved numerically using the Matlab-ode15s, while the finding of least squares was realized with the Matlab-globalsearch. To keep the calculation time manageable and simultaneously ensure the accuracy, we have chosen n = 1000 and the φ-step equal to φ * /7335. We have found that the torque-twist curve for sample 1 taken from [14] can be fit with just a single set of system parameters. These are: T P = 28911 K, s = 0.0156, χ 0 = 0.243, K ρ = 49.2, K = 379, k 0 = 4.31 × 10 −7 , k 1 = 1.79 × 10 8 ,ρ i = 6 × 10 −5 , andχ i = 0.178. The agreement between theory and experiment seems to us to be well within the bounds of experimental uncertainties. Even the initial yielding transition appears to be described accurately by this theory. There is only one visible discrepancy: at large twist angles (φ > 69 • ) the torques are slightly below those predicted by the theory. Nothing about this result leads us to believe that there are relevant physical ingredients missing in the theory. The results of numerical simulations for other quantities are shown in Figs. 3-7. We plot in Fig. 3 the shear stress distribution τ = µ(rφη − β) at three different twist angles φ = 10 • (black), φ = 30 • (red/dark gray), and φ = 50 • (yellow/light gray). Contrary to the similar distribution obtained by the phenomenological theory of ideal plasticity, the stress in the plastic zone does not remain constant, but rises with increasing r and reaches a maximum at r = R. This exhibits the isotropic hardening behavior due to the entanglement of dislocations. Fig. 4 shows the evolution of the plastic distortion β(r) at the above three different twist angles. It can be seen that the plastic distortion is an increasing function of r except very near the free boundary r = R. Since the latter attracts excess dislocations, β(r) should decrease in this region to ensure equilibrium. However, due to the strong external stress field, the influence of this attraction can only be felt in a thin layer near the free boundary. In our approximate finite difference solution the decrease of β(r) occurs between the end-point and the fictitious point and cannot be seen in Fig. 4. Figs. 5 and 6 present the densities of total and excess dislocations, respectively, at the above three different twist angles. Under the applied shear stress, the excess dislocations of the positive sign move to the center of the bar and pile up there. At large twist angles the distribution of excess dislocations over radius r remains almost constant except near the center and the free boundary (cf. [18][19][20]).
To understand the mechanism of formation of excess dislocations, we plot in Fig. 7 distribution ρ g (r) at a small twist angle φ = 0.1 • . Since the flow stress at this twist angle exceeds the Taylor stress, redundant dislocations in the form of dislocation dipoles begin to dissolve according to the kinetics of thermally activated dislocation depinning [1][2][3][4][5]. Under the applied shear stress, positive dislocations then move towards the center and nega- At small angles of twist, the applied shear stress near the center is still small and cannot move dislocations. Therefore, excess dislocations occupy an outer ring, as can be seen in Fig. 7. As the angle of twist increases, the shear stress increases as well, and when it becomes large enough, it can drive these excess dislocations to the center and they pile up there. Thus, we can say that the dissolution of dipoles near the free boundary results in excess dislocations of positive sign. They then move to the center and pile up there, increasing kinematic hardening.
It is interesting to examine the influence of the size of the sample on the torque-twist curve. Fig. 8 shows the three normalized torque T /R 3 (measured in Pa) versus normalized twist ωR curves for three bars with different radii R = 25 micron (black), R = 50 micron (red/dark gray), and R = 100 micron (yellow/light gray). We choose the maximum strain rateωR = 10 −3 /s, while all other parameters are left unchanged. We see that the size strongly influences the slope of the hardening curve, since the accumulated excess dislocations pile up against the center leading to a stronger kinematic hardening for the smaller sample than for the larger one (smaller is stronger). The yielding transition, on the other hand, is almost independent of the radius. This can be explained by the fact that at the onset of yielding transition practi- cally no excess dislocations occur, so that the kinematic hardening is not yet noticeable.
Another important question is how strongly the twist rate affects the torque-twist curve. Fig. 9 shows the three torque-twist curves for three samples loaded at three different twist ratesφ = 0.25 • /s (black),φ = 2.5 • /s (red/dark gray), andφ = 25 • /s (yellow/light gray). The radius of the samples is R = 6.35mm, while all other parameters remain unchanged. We see that the twist rate mainly affects isotropic hardening: the higher the twist rate, the higher the slope of the torque-twist curve. The kinematic hardening is not affected by the change of the twist rate. The reason for this is that the kinematic hardening due to the excess dislocations is much less sensitive to the change in strain (twist) rate.

V. CONCLUDING REMARKS
Overall, these results seem to us to be quite satisfactory. Note that we now use thermodynamic dislocation theory for non-uniform deformations not just to test its validity but also as a tool for discovering properties of structural materials. For example, we could find the mechanism of forming excess dislocations based on the dissolution of dislocation dipoles near the free boundary of the bar and predict their distribution. One of the main reasons for the success of this theory -as has been emphasized here and in earlier papers -is the extreme sensitivity of the plastic strain rate to small changes in the temperature or the stress. Another reason for its success is the inclusion of the excess dislocations in the theory, which leads to size-dependent kinematic hardening. Here, in our opinion, the incompatible plastic distortion is the natural variable that keeps the memory of excess dislocations. It cannot enter the free energy, but the curl of this quantity should enter the free energy causing the back stress. In this way the theory differs substantially from the phenomenological plasticity that introduces the back stress along with an assumed constitutive equation to fit the stress strain curves with kinematic hardening. On the contrary, our theory allows us to find the back stress from the first principle calculation of the free energy of dislocated crystals.
The results obtained show the principal applicability of TDT to non-uniform plastic deformations. As far as the size effect is concerned, we could not find reliable experimental data for single crystal copper under torsion at different bar radii, in contrast to polycrystalline copper under torsion [12]. However, the proposed theory may serve as a useful guide for the future experimental investigation of the torsion of single crystal bars in several directions: (i) the torque-twist curves at load reversals and the analog of the Bauschinger effect, (ii) the size effect, (iii) the sensitivity of the torque-twist curves to the twist rate and temperature, et cetera. The identification of material parameters for polycrystalline copper under torsion and the comparison with experiments in [12] will be addressed in our forthcoming paper.